Växjö University

(1)

Växjö University

School of Mathematics and System Engineering Reports from MSI - Rapporter från MSI

Some recent simulation techniques for diffusion bridges

Yadigar Sekerci

May 2009

MSI

Växjö University SE-351 95 VÄXJÖ

Report 09018 ISSN 1650-2647

ISRN VXU/MSI/MA/C/--09018/--SE

(2)

(3)

Yadigar Sekerci

Some recent simulation techniques for diffusion bridges

Bachelor’s thesis Mathematics

2009

Växjö University

(4)

(5)

Abstract

We apply some recent numerical solutions to diffusion bridges written in Iacus (2008).

One is an approximate scheme from Bladt and Sørensen (2007), another one, from Beskos et al (2006), is an algorithm which is exact: no numerical error at given grid points!

Key-words: Brownian bridge, diffusion bridge, Brownian motion, stochastic differential equation, simulation, numerical methods, Euler-Maruyama’s method, acceptence- rejection method.

Sammanfattning

Vi tillämpar några nya numeriska lösningar till diffusionsbryggor beskrivna i Iacus (2008).

En algoritm, från Bladt Sørensen (2007), är enkel och approximativ, den andra, från Be- skos et al (2006), är exakt: inga numeriska fel vid givna gridpunkter!

Nyckelord: Brownsk rörelse, diffusionsbro, Brownsk rörelse, stokastisk differentialekva- tion, simulering, numeriska metoder, diffusionsbro, Euler-Maruyama’s method, acceptans- förkastelsemetoden

Acknowledgments

I thank Roger Pettersson

1 Introduction

Much of this work is inspired by nice examples of [9] where directions in the field of simulation of stochastic differential equations are presented.

2 Simulation of the Brownian motion

A Brownian motion, also called Wiener process, is a continuous process {Xt}_t≥0 starting at zero with independent stationary increments that are normal distributed, more exactly, X_t+h− X_t are normal distributed with zero mean and standard deviation√

h. One heuristic representation of the Brownian motion is that it satisfies the differential equation

dX = dW, X₀= 0 (2.1)

where dW is√

dtN(0,1), a normal random variable with standard deviation√

dt! In fact the representation (2.1) is a particular case of stochastic differential equations (SDEs) often written in the form

dX= µ(X )dt + σ (X )dW, X₀= a for which there are stringent definitions in its integral form.

It is simple to simulate the Brownian motion during a time interval [0, T ] at gridpoints 0 = t₀< . . . < t_n= T . The algorithm, with the miscellanous notation X_i= X_t_i= X (t_i) is simply

X_i+1= X_i+p

∆_iN(0, 1),

where ∆_i= ∆ti = t_i+1− t_i. It means that for each time step you simulate a normal random variable with mean zero and standard deviation√

∆_i.

3 Simulation of the Brownian bridge

The Brownian bridge is a process {X_t} that heuristically is just a Brownian motion with fixed end points X₀= a, X_T = b. It satisfies the SDE

dX= dW, X₀= a, X_T = b. (3.1)

Note that X_t is not adapted, it means X_t not only depends on the information (the σ - algebra) of the Brownian motion W up till time t. Xt depends on the future values. That means that even though (3.1) looks almost trivial it is actually not a standard SDE. How- ever, it is well known that a solution to (3.1) is given by the closed formula,

X_t= a(1 − t

T) + bt

T + (W_t− t

TW_T), 0 ≤ t ≤ T, (3.2) see any standard book, for example [2]. That makes is particularly simple to calculate the Brownain bridge: First you simulate the Brownian at grid points 0 = t₀< t₁< . . . < t_n= T . Then you simulate the Brownian bridge at those grid points:

X_i= a(1 −t_i T) + bt

T + (W_i− t

TW_n), 0 ≤ t_i≤ T, (3.3)

(8)

Figure 3.1: The blue path (solid) is an exact simulation of the Brownian bridge at the grid points by the algorithm (3.3) starting and ending at zero. The green path (shadowed) is the Euler approximation (3.5) of (3.4). The number of gridpoints is 1000. The error seems to increase. For this simulation the error seems to be around the middle. When I simulated I found that the error is typically largest around the middle. It would be intersting to know a more exact description of the error.

where we again used the miscellanous notation X_i= X_t_i and W_i= W_t_i. Now what is interesting, is that the Brownian bridge satisfies the SDE

dX= b− X_t

T− t dt+ dW_t, 0 ≤ t ≤ T, X₀= a. (3.4) Note that there is no end-point condition X_T = b. The SDE (3.4) looks like an SDE of standard type with solution X adapted to the Brownian motion, i.e. X_tdoes not depend on future values of W . However, the drift coefficient (b − X_t)/(T − t) is dangerously large as t approaches T . Now, an interesting question is, how bad is an Euler simulation of (3.4).

In this bachelor thesis we don’t reply that answer. Note that is not a trivial matter due to the non-boundedness of the drift coefficient for t near T . The Euler simulation is (again with miscellanous notation)

X_i+1= X_i+b− X_i

T− t_i∆i+p

∆iN(0, 1). (3.5)

4 Simulation of stochastic differential equations

A generalization of the Brownian motion is a solution to the stochastic differential equation (SDE)

dX_t = µ(X_t)dt + σ (X_t)dW_t, X₀= a. (4.1) Nowadays there is a huge literature on such equations. One interpretation for tourists in this field is simply that it is the limit of the following naïve Euler-Maruyama algorithm:

(9)

X₀= a and

X_i+1= X_i+ µ(Xi)∆i+ σ (Xi)p

∆iN(0, 1), i = 1, 2, . . . .

As before X_iis X_t_iwhere t₀= 0 < t₁< . . ., ∆i= t_i+1−t_i. In the next chapter on simulation of diffusion bridges, Euler-Maruyama simulation of SDEs will be used in an approximative algorithm for diffusion bridges. There are other more sophisticated methods, see e.g [9] and [5].

Throughout we assume there exists a unique solution of (4.1). Classical assumptions that give existence and uniqueness of a solution to (4.1) is linear growth and a Lipschitz continuity of b and σ¹.

5 Simulation of diffusion bridges

A generalization of (2.1) is

dX_t = µ(Xt)dt + σ (Xt)dW_t, 0 < t < T, X(0) = a, X(T ) = b (5.1) Given µ and σ it is often convenient to call a solution to (5.1) a (0, a, T, b)-bridge.

It is not a trivial matter how to simulate (5.1). In the literature it seems that one im- portant reason why diffusion bridge simulation is needed is that it is used for likelihood inference of discretely observed diffusion processes. That means that simulation of diffusion bridges is used for estimating parameters of SDEs observed at some grid points.

In [3] it is suggested a particular technique based on Markov Chain Monte Carlo (MCMC) methods. In [7] there is a technique emaneting from the observation that limit solutions of certain stochastic partial differential equations are solutions to diffusion bridges! In [4] it is suggested a technique which is exact at any given time points the time interval [0, T ]! In this thesis we apply the approximate technique from [8] and the technique from [4].

5.1 The approximate algorithm

In [8] there is a simple algorithm to a simulate an approximate solution to (5.1), it means an approximate diffusion bridge. For the algorithm to work, time-reversibility of the diffusion process {X_{t∈[0,T ]}} is assumed, that means {X_t− X₀: t ∈ [0, T ]} and the process {X_T_−t− X_T : t ∈ [0, T ]}, have the same statistical properties: the process behaves similarly if the time is running backwards.

We impose the following condition, which ensures that the diffusion X is time-reversible, see [8], [5]. For that condition we let (l, r), −∞ ≤ l < r ≤ ∞ be the state space of {X_{t∈[0,T ]}}, that means l is the smallest possible value for {X_{t∈[0,T ]}} which may be −∞ and r is the largest possible value for {X_{t∈[0,T ]}} which may be ∞.

Condition5.1. Let {X_{t∈[0,T ]}} be the solution to the SDE (4.1). Let

s(x) = exp(−2 Z _x

x#

µ (y) σ²(y)dy) where x# denotes an arbitrary point in the state space (l, r). Let

m(x) = 1

σ²(x)s(x).

1f Lipschitz continuous means | f (x) − f (y)| ≤ L|x − y| for all x, y.

(10)

Figure 5.1: The figures describes the approximate algorithm for the diffusion bridge dX =

−0.5Xdt + dW , X₀= 0, X_T = 1. In the upper figure the blue one (solid line) starts at 0, the green one (dotted) starts at 1 (!) and the red one (dashed) is the flipped version of the green one, it ends at 1. The bridge starts with the blue (solid) one and if it crosses the red (dashed) one the bridge continues with the red one which finishes at 1.

Assume

Z r l

m(x)dx < ∞ and

Z _x

x#

s(x)dx = Z x#

l

s(x)dx = ∞. (5.2)

The approximative algorithm from [8] for simulation of diffusion bridges is as follows.

Let W¹and W²be two independent Wiener processes, and let X¹ and X² as the solutions to

dX_t¹= µ(X_t¹)dt + σ (X_t¹)dW_t¹, X₀¹= a and

dX_t²= µ(X_t²)dt + σ (X_t²)dW_t², X₀²= b.

The idea is to use the time-reversibility of X² and realize a (0, a, T, b) -bridge by simulating the process X¹ from a forward in time and X² from b backward in time! If the samples paths of the two process intersect, they can be combined into a realization of the bridge.

Fact 5.1 ([8]). Let τ = inf (0 ≤ t ≤ 1X_t¹= X_1−t² ), where the infimum over empty set is taken to be infinite. Define

Z_t= X_t¹, if 0 ≤ t ≤ τ X_1−t² if τ < t ≤ 1

Assume Condition 5.1. Then, the distribution of {Z_t: 0 ≤ t ≤ T } conditional on the event {τ ≤ T } is equal to the conditional distribution of {X_t : 0 ≤ t ≤ T } given X₀= a and X_T = b; i.e., Z is a (0, a, T, b)-diffusion bridge.

The fact 5.1 points towards a natural algorithm for simulation of a diffusion bridge.The algorithm consists in simulating two independent diffusion process X¹and X² using one numerical method, for example the Euler-Maruyama on the time interval [0, T ] with discretization step ∆ = T /n applying a rejection sampling procedure. Let Y_i∆¹ and Y_i∆², i= 0, 1, 2, ..., n be independent simulations of X¹ and X². If either Y_i∆¹ ≥ Y_(n−i)∆² and Y_(i+1)∆¹ ≤ Y_{(n−(i+1))∆}² or Y_i∆¹ ≤ Y_(n−i)∆² and Y_(i+1)∆¹ ≥ Y_{(n−(i+1))∆}² , a crossing has been realized. Hence, let ν = min{i ∈ (1, 2, ..., n) | Y_i∆¹ ≤ Y_(n−i)∆² } if Y₀¹ ≥ Y_T² and v = min{i ∈ (1, 2, ..., n) | Y_i∆¹ ≥ Y_(n−i)∆² } if Y₀¹≤ Y_T², and define

B_i∆= Y_i∆¹, if i = 0, 1, . . . , ν − 1 Y_(n−i)∆² if i = ν, ..., n

Then B is simulation of a (0, a, T, b)-diffusion bridge. If no crossing happened,start again by simulating Y_i∆¹ and Y_i∆² and iterate until a crossing of the two trajectories is realized.

(11)

Figure 5.2: 20 simulated OU-bridges with the algorithm exact at the grid points.

Figure 5.3: 20 simulated OU-bridges with the approximate algorithm

5.2 Ornstein-Uhlenbeck bridge

For the Ornstein-Uhlenbeck bridge (OU-bridge)

dX_t= −θ Xtdt+ σ dWt, X₀= a, X_T = b (5.3) there is a simple algorithm that is exact at the grid points, [8]. According to [8], from the transition densities of the Ornstein-Uhlenbeck process we can calculate the transition densities of the OU-bridge. Thus we could in principle simulate the OU-bridge by sampling transitions from these densities. The following alternative method is, however, numerically more stable. For this particular equation there exists an algorithm which at the grid points t₁,t₂· · · is exactly an OU-bridge. That means that at the grid points the simulated process has the same distribution as the OU-bridge. The algorithm from the lemma below is from [8]:

Lemma 5.1. For a partition 0 = t₀< t₁ < . . . < t_n= T , define X_i = X_t_i and Y_i= Y_t_i by X₀= a,

X_i+1= e^{−θ (t}ⁱ⁺¹^−tⁱ⁾X_t_i+Y_i, i= 1, . . . n where

Y_i+1= (σ²

2θ(1 − e^{−2θ (t}ⁱ⁺¹^−tⁱ)))^1/2N(0, 1), i= 1, . . . n

where N(0,1) represents independ N(0,1) random variables. Let Z_ibe defined by

Z_i= X_i+ (b − X_n)e^{θ t}ⁱ− e^−θtⁱ

e^{θ t}ⁿ− e^−θtⁿ, , i= 1, . . . n.

Then Z_t_i= Z_ihas the same distribution as a Ornstein-Uhlenbeck bridge at the grid points 0 = t₀< t₁< . . . < t_n= T , with Z₀= a, Z_T = b.

Lemma 5.1 makes it possible to compare the approximate algorithm with the OU- algorithm exact at the grid points. See figures 5.1–5.3 .

Figure 5.4: Kernel density estimate of the Ornstein-Uhlenbeck bridge at midpoint time t = 1/2 with the algorithm exact at the grid points and the kernel density estimate of the OU-bridge with the approximate algorithm. For the approximate algorithm the number of time steps is 500 and the number of simulations is 100²= 10 000. For the exact algorithm only the mid-point time is simulated (because the random draw is exact) and the number of simulations is also 10 000, due to slowness of the exact algorithm.

(12)

5.3 One exact algorithm

Let us first describe an algorithm for SDEs (4.1). We later describe a modified algorithm for diffusion bridges (5.1).

In [4] is suggested a general technique for simulation of SDE (4.1) for which the simulated process has the same distribution as the process given by (4.1). It is not, as for example, the Euler-Maruyama method only an approximative algorithm.

The conditions we apply here is the following condition.

Condition5.2. (i) The derivative µ_x of µ exists

(ii) There exists k₁and k₂such that k₁≤ ¹₂µ²(x) +¹₂µ_x(x) ≤ k₂for any x∈ R

Condition 5.2 can be relaxed but then the algorithm is more technical. In the algorithm below we use the function φ (x) =¹₂µ²(x) +¹₂µ_x(x) − k₁. Note that 0 ≤ φ (x) ≤ M = k₂−k₁ The exact algorithm for simulating (4.1) is briefly as follows. For a partition {0 = t₀ < t₁< · · · < t_n = T }, given that X_t_i is already simulated, X_t_i+1 is first simulated in a certain way, typically by a acceptance-rejection rule (not a variant of the Euler-Maruyama method). To simulate the process in [ti,t_i+1], a Poisson random number (an integer) is simulated. If that Poisson random integer, N say, is positive, N couples of uniformly random variables {(τ₁, v₁), . . . , (τN, v_N)}, where τ1, . . . , τN are uniform random numbers in [ti,t_i+1] and (v₁, . . . , v_N) are uniform random in the set [0, M]. Then simulate a ’standard’

Brownian bridge in [t_i,t_i+1] beginning at X_t_i and ending at X_t_i+1. Consider the values y₁, . . . y_N of the Brownian bridge at the time points τ₁, . . . , τN. The Brownian bridge is accepted as a simulation of {Xt : t ∈ [ti,t_i+1]} if φ (y1) ≤ v₁, . . . , φ (yN) ≤ v_N.

If we only want values of {X_t} at the grid points 0 = t₁< t₂< . . . < t_n= T , the connecting Brownian bridges are not needed. An alternative would instead be to only simulate the end point XT exactly by an acceptance-rejection rule, and then simulate the Brown- ian bridge in the whole interval [0, T ] and, if the Brownian bridge is accepted, the values of the Brownian bridge at the grid points is obtained as simulated values of the process {X_t : t ∈ [0, T ]}. However, the rejection probability of the XT-value and the Brownian bridge is in that situation probably high which should imply an inefficient algorithm.

For diffusion bridges of type (5.1) the exact algorithm is simpler. Given X₀= a, X_T does not need to be simulated by a time consuming acceptance-rejection rule: XT is just given the value b. Instead the skeleton time points τ₁, . . . , τk can be uniformly picked in the whole time interval [0, T ].

It is applied to diffusion bridges of the type

dX_t = µ(Xt)dt + dW_t, X₀= a, X_T = b, (5.4) i.e. σ ≡ 1.

Note that (5.4) is more general than one may think at the first sight. The equation (5.1) with possibly non-constant σ can, under non-degeneracity conditions of σ (for example σ (x) > ε > 0 for any x) be transformed into one with diffusion coefficient one by applying the Lamperti transform,

Y_t= F(X_t) = Z _X_t

z

1 σ (u)du.

Here z is any arbitrary value in the state space of {X_t}. Indeed, the process {Y_t} solves the SDE

dY_t = µY(Y_t)dt + dW_t, where

µY(Y_t) = µ (F⁻¹(y)) σ (F⁻¹(y) −1

2σx(F⁻¹(y))

(13)

Figure 5.5: 20 simulated bridges of dX = sin(X )dt + dW , X₀= X₁= 0 with the algorithm exact at the grid points.

Figure 5.6: 20 simulated bridges of dX = sin(X )dt + dW , X₀= X₁= 0 with the approximate algorithm

which can also be written as

dY_t= (µ (Xt) σ (X_t)−1

2σ_x(X_t))dt + dW_t.

The ’exact’ algorithm to simulate the diffusion bridge (5.4), from [4], see also [9]:

1. Let X_T = b.

2. Simulate a non-negative integer τ = N from the Poisson distribution with intensity λ = T M.

3. Draw (if N ≥ 1) k pairs of uniformly distributed random numbers {(τ_i, v_i) : i = 1, . . . , N} on [0, T ] × [0, M],i=1,2,...,N.

4. Generate a Brownian bridge X starting at x at time 0 and ending in b at time T at time instants t_i; i.e., generate X_τ_i= y_i, i = 1, ..., N.

5. Compute the indicator function

I=

N

∏

i=1

1_{{φ (y}_i_)≤ν_i_}.

6. If I = 1, the trajectory {Y₀= a,Y_t₁ = y₁, ...,Y_t_k= y_k,Y_T = b} is accepted as an exact draw of a solution X to (5.4) at the points 0 < t₁< . . . < t_k < T . Otherwise, restart from step 1.

Exact draws from X given by (5.4) at pre-determined grid-points 0 = t₁< t₂< . . . <

t_n= T is obtained by noting the values at these points of the connecting Brownian bridges between the points τ₁< . . . < τ_k.

Figure 5.7: Kernel density estimate of the bridge dX = sin(X )dt + dW , X₀= X₁= 1 at midpoint time t = 1/2 with the algorithm exact at the grid points and the kernel density estimate of the bridge with the approximate algorithm. For the exact algorithm the number of time steps is only 3 (we only need the midpoint!) and the number of simulations is 500. For the approximate algorithm the number of time steps is 501 and the number of simulations only 500. In this situation I don’t know a closed formula for the density of X_1/2, the midpoint value.

(14)

Figure 5.8: Simulation of a numerical solution to (5.6),[6]. The solution to (5.6) converges to a diffusion bridge.

5.4 PDE-MCMC-algorithm

From [7], for X given by the diffusion bridge SDE

dX(x) = µ(X (x))dx + σ dW (x), X(0) = X (1) = 0 (5.5) (here the time variable is x; the reason will be clear from the arguments below) can be sampled by MCMC-sampling of u given by stochastic partial differential equation







∂ u

∂ t = ¹

σ²(^∂²^u

∂ x²− µ(u)µ⁰(u) −¹₂µ⁰⁰(u)) +√ 2^∂²^W

∂ t∂ x

u(0,t) = u(1,t) = 0, t∈ [0, T ] u(x, 0) = u₀(x), x∈ [0, 1].

(5.6)

where ^∂²^W

∂ t∂ x is a time-space Gaussian white noise.

MCMC (Markov Chain Monte Carlo) is a particular class of algorithms for sampling from probability distributions based on constructing a Markov chain that has the desired distribution as its equilibrium distribution. The state of the chain after a large number of steps is then used as a sample from the desired distribution. The quality of the sample improves as a function of the number of steps. What is used here is that (for a fixed x∈ [0, 1], as t → ∞, u(x,t) converges to a process X(x) which solves (5.5). What is needed is to simulate u(t, x) for large t in a clever way. Normally when simulations are needed the time interval is bounded, but here the focus is on for large times. For that the, for me complicated MCMC-algorithm, is needed.

References

[1] S: Iacus, Simulation and inference for stochastic differential equations, with R examples, Springer, (2008).

[2] , I Karatzas, I Shreve, Brownian motion and stochastic calculus, Springer (1991).

[3] G. Roberts, O. Stramer, On inference for partial observations from partial observed nonlinear diffusion model using Metropolis-Hastings algorithms, Biometrika 88:3 (2001), 603-621.

[4] A.Beskos, O.Papaspiliopoulos and G.O. Roberts, Retrospective exact simulation of diffusion sample paths with applications, Bernoulli 12:6 (2006), 1077-1098

[5] J. Kent, Time-reversible diffusions, Adv Appl Probab 10 (1978), 819-835

[6] R. Pettersson, Projection scheme for a reflected stochastic heat equation with additove noise, Foundations of Probability and Physics-3, AIP Conference Proceedings, ed: A.

Yu Khrennikov, (2005), 158-167.

[7] A.M. Stuart, J. Voss, and P. Wiberg Conditional path sampling of SDEs and the Langevin MCMC method. Comm. Math. Sci. 2 (2004), 685-697.

(15)

[8] M.Bladt and M.Sorensen, Simple simulation of diffusion bridges with application to likelihood inference for diffusions,Centre for Analytical Finance, Working Paper, University of Copenhagen. February (2007)

[9] S. M. Iacus, Simulation and interference for Stochastic Differential Equations, with R Examples, Springer Series in Statistics, (2008)

A Appendix: Definition of random variables etc

A random variable X is a function Ω 3 ω 7→ X (ω) measurable in a certain sense so that probabilities can be defined for sets of the type {X (ω) ≤ x}. The ω represents a random choice. A real-valued random variable simply takes its values in R. A continuous stochastic proces X = {X_t: t ≥ 0} is a random variable with values in the set of continuous functions C(R). The map then is Ω = C(R) 3 Xt(ω) = ω(t), i.e. the random selection is then a continuous function.

A Brownian motion is a continuous stochastic process {Wt : t ≥ 0}, W (0) = 0, with independent increments W_t+h−W_t that are√

hN(0,1) for h > 0.

An Ornstein-Uhlenbeck process is a process {X_t≥0} satisfying dX_t= −θ Xtdt+ σ dWt, X₀= a i.e. as the OU-bridge (5.3) but without restrictions on X_T.

B Appendix: matlab codes

B.1 Euler-Maruyama’s method for ’standard’ SDEs Matlab code for Euler simulation of SDE of the type

dX= µ(X )dt + σ (X )dW, X₀= x.

In the code below mu(x) = sin(x), σ (x) = 1.

x0=1; % initial value

n=1000; % number of partition points t=linspace (0,1,n); %gridpoints

dt=1/(n-1); %partition step

N=500; % N number of simulations x=x0*ones (1,N);

mu=@(x) sin(x);

sigma=@(x)1;

dW=sqrt (dt)*randn (n,N);

for i=1:n-1;

x (i+1,:)=x (i,:)+mu(x (i,:))*dt+sigma (x (i,:)).*dW (i,:);

end

B.2 Brownian bridge

Matlab code for Brownian bridge X_t= a(1 − t

T) + bt

T + (W_t− t

TW_T), 0 ≤ t ≤ T,

(16)

n=10000; % number of partition points T=1;

t=linspace(0,T,n)’;

dt=1/(n-1);

N=5;

dW=sqrt (dt)*randn (n,N);

W=cumsum(dW);

a=0;

b=0;

X=a*(1-t/T)+b*t/T+(W-t/T.*W(end,:)’);

B.3 Approximate diffusion bridge

The matlab code for the approximate diffusion bridge that appeared in [8].

n=501; %number of partition points N=500; %number of simulations

ok=0;

dt=t(2)-t(1);

bridges=zeros(n,N);

for j=1:N ok=0;

while ok==0;

x=a;y=b; z=a;

dWx=sqrt (dt)*randn (n,1);

dWy=sqrt (dt)*randn (n,1);

for i=1:n-1;

x (i+1)=x (i)+mu(x (i))*dt+sigma (x (i)).*dWx (i);

y (i+1)=y (i)+mu(y (i))*dt+sigma (y (i)).*dWy (i);

end

z=fliplr(y);

dummy=(x(1:end-1)>=z(1:end-1)).*(x(2:end)<=z(2:end))...

|(x(1:end-1)<=z(1:end-1)).*(x(2:end)>=z(2:end));

if max(dummy)==1 ok=1; %they meet

index=min(find(dummy==1));

end end

bridges(:,j)=[x(1:max(index-1,1)) z(max(index,2):end)];

end

B.4 Code for simulation of Ornstein-Uhlenbeck bridge exact at grid points Matlab code for simulation of

dX = −θ (Xt)dt + σ (dWt), X₀= a, X_T = b exact at the grid points, [8].

n=100; % number of partition points N=n^2; %number of simulations

(17)

T=1;

t=linspace(0,T,n)’; %grid points dt=1/(n-1);

a=0;

b=1;

theta=.5;

sig=1;

Y=(sig/sqrt(2*theta)*sqrt(1-exp(-2*theta*diff(t)))*ones(1,N)).*...

randn(n-1,N);

X=a*ones(1,N);

for i=1:n-1

X(i+1,:)=exp(-theta*dt)*X(i,:)+Y(i,:);

end

Z=X+(exp(theta*t)-exp(-theta*t))/(exp(theta*t(n))-...

exp(-theta*t(n))*(b-X(n,:));

Algorithm for more general diffusion bridges exact at grid points Matlab code for

dX_t = µ(Xt)dt + dW_t, X₀= a, X_T = b, i.e. σ ≡ 1, where k₁≤ ¹₂µ²(x) +¹₂µ_x(x) ≤ k₂

N=500;%number of simulations n=3;%number of partition points mu=@(x)sin(x);

mup=@(x)cos(x);

sigma=@(x)1;

a=0;

b=0;

T=1;%final time

prephi=@(x)1/2*mu(x).^2+1/2*mup(x);

mprephi=@(x)-prephi(x);

k1=prephi(fminsearch(prephi,-pi));%-1/2;

k2=prephi(fminsearch(mprephi,pi/2));%5/8;

phi=@(x)prephi(x)-k1;

M=k2-k1;

lambda=T*M;

Ndifb=zeros(n,N);

for j=1:N accept=0;

while accept==0

k=poissrnd(lambda);

if k>0

pretskel=sort(T*rand(k,1));

v=M*rand(k,1);

tskel=[0; pretskel; T];

W=[0; cumsum(sqrt(diff(tskel)).*randn(length(tskel)-1,1))];

skeleton=a+W-tskel/T*(W(end)-b+a);

I=sum((phi(skeleton(2:end-1))<=v));%prod

(18)

if I==k accept=1;

end end end

% [tskel skeleton]

difb=a;

for i=2:length(skeleton)

if max((t>tskel(i-1)&t<tskel(i)))>0 s=t(t>tskel(i-1)&t<tskel(i));

difb=[difb; Bbfn(tskel(i-1),skeleton(i-1),tskel(i),... skeleton(i),s)];

end end

difb=[difb;b];

Ndifb(:,j)=difb;

%plot(t,difb)

%plot(ti,x);

end

(19)

(20)

Växjö universitet

Matematiska och systemtekniska institutionen SE-351 95 Växjö

Tel. +46 (0)470 70 80 00, fax +46 (0)470 840 04 http://www.vxu.se/msi/

Växjö University

Växjö University

Some recent simulation techniques for diffusion bridges

Yadigar Sekerci

Some recent simulation techniques for diffusion bridges

Växjö University

Abstract

Sammanfattning

Acknowledgments

Contents

1 Introduction

2 Simulation of the Brownian motion

3 Simulation of the Brownian bridge

4 Simulation of stochastic differential equations

5 Simulation of diffusion bridges

∏

References

A Appendix: Definition of random variables etc

B Appendix: matlab codes

Växjö universitet