Numerical Methods of SDEs-Dr. Raúl Tempone Homework #3

Exercise 1 Formulate and motivate a forward Euler method for approximation of the Stratonovich SDE

dX_{t} = a(t, X_{t})dt + b(t, X_{t}) ∘ dW_{t}  (1)

First approach: (Incorrect approach!) We approximated b(t, X_{t}) by the midpoint rule:

X_{t_{n+1}} = X_{n} + a(t_{n}, X_{t_{n}})Δt_{n} + b (\frac{1}{2}(t_{n+1} + t_{n}), \frac{1}{2}(X_{t_{n+1}} + X_{t_{n}})) ΔW_{n}  (2)

What is wrong with approach (2)?

Consider the following process as an example to illustrate why formulation (2) was incorrect.

dX(t) = bX_{t}dW (t) with a constant b,  (3)

Using (2) to approximate (3) we obtain:

X_{t_{n+1}} = X_{t_{n}} + b (\frac{1}{2}(X_{t_{n+1}} + X_{t_{n}})) ΔW_{n}
X_{t_{n+1}}(1 − \frac{1}{2}bΔW_{n}) = X_{t_{n}}(1 + \frac{1}{2}bΔW_{n})

X_{t_{n+1}} = X_{t_{n}} [\frac{(1 + \frac{1}{2}bΔW_{n})}{(1 − \frac{1}{2}bΔW_{n})}]  (4)

Consider the computation of E(X_{1}^{2})

E(X_{1}^{2}) = \frac{(X_{0})^{2}}{\sqrt{2πΔt}} ∫_{−∞}^{∞} \frac{(1 + \frac{1}{2}bz)^{2}}{(1 − \frac{1}{2}bz)^{2}} e^{\frac{−z^{2}}{2Δt}} dz  (5)

where we have used the fact that E(f(x)) = ∫ (f(x)dp(x)) where p(x) is the probability density of x.

From (5), consider a similiar integral:

∫_{−∞}^{∞} \frac{(1 + x)^{2}}{(1 − x)^{2}} e^{−x^{2}}dx  (6)

However, the integral in (6) does not converge.

Formulation: Rewrite (1), so that may discretize using the forward Euler method.

X_{t} = X_{0} + ∫_{0}^{t} a(s, X(s))ds + \frac{1}{2} ∫_{0}^{t} b(s, X_{s})b_{x}(s, X_{s})ds + ∫_{0}^{t} b(s, X_{s})dW_{s}  (7)

From (7), considering one interval, [t_{n}, t_{n+1}]

X_{n+1} = X_{n} + ∫_{t_{n}}^{t_{n+1}} a(x, X(s))ds + \frac{1}{2} ∫_{t_{n}}^{t_{n+1}} b(s, X_{s})b_{x}(s, X_{s})ds + ∫_{t_{n}}^{t_{n+1}} b(s, X_{s})dW_{s}  (8)

Discretization of (8) to obtain the following:

X_{t_{n+1}} = X_{t_{n}} + a(t_{n}, X_{t_{n}})Δt_{n} + \frac{1}{2}b(t_{n}, X_{t_{n}})\frac{∂b(t_{n}, X_{t_{n}})}{∂x}Δt_{n} + bΔt_{n}ΔW_{t_{n}}  (9)

Exercise 2 Consider the deterministic differential equation

dZ_{t} = a(Z_{t})dt  (10)
Z(0) = x_{0}  0 ≤ t ≤ T,

and consider a peturbation of (10), the Itô stochastic differential equation

dX_{t} = a(X_{t})dt + bdW_{t}  (11)
X(0) = x_{0}  0 ≤ t ≤ T,

where a is a smooth function and b > 0 is a positive constant. The aim of this exercise is to compare the solution of both equations. Define then the difference:

e_{t} = X_{t} − Z_{t}  (12)

a) Consider a(x) = ax (linear case) and compute E(e_{t}) and Var(e_{t}).

Hint: Use Itô’s formula when necessary. From (10), obtain the following solution

Z_{t} = x_{0}e^{at}  (13)

Applying the intergrating factor, define the following:

Y_{t} = e^{−at}X_{t}  0 ≤ t ≤ T  (14)
Y_{0} = x_{0}

Itô’s formula. Let f be a smooth function, assume that X_{t} is the solution of an SDE, then g_{t} = f(X_{t}, t) satisfies the SDE

df(X_{t}, t) = \frac{∂f}{∂t} dt + \frac{∂f}{∂x} dX_{t} + \frac{1}{2} \frac{∂^{2}f}{∂x^{2}} σ^{2}(X_{t}, t)dt,  (15)

Applying (15) to (14), to obtain the following:

dY_{t} = −ae^{−at}X_{t}dt + e^{−at} \underbrace{dX_{t}}_{see (11)} + \underbrace{higher order terms}_{=0}

= −ae^{−at}X_{t}dt + ae^{−at}X_{t}dt + be^{−at}dW_{t}
dY_{t} = be^{−at}dW_{t}  (16)

integrating (16), will obtain:

Y_{t} = Y_{0} + ∫_{0}^{t} be^{−as}dW_{s}

e^{−at}X_{t} = X_{0} + ∫_{0}^{t} be^{−as}dW_{s}
X_{t} = e^{at}X_{0} + ∫_{0}^{t} be^{a(t−s)}dW_{s}  (17)

Therefore, now that we have obtained an equation for X_{t}, see (17), revisit the definition for the difference, (12).

e_{t} = \underbrace{X_{t}}_{see (17)} − \underbrace{Z_{t}}_{see (13)}

e_{t} = ∫_{0}^{t} be^{a(t−s)}dW_{s}  (18)

Solution: Now derived an equation for e_{t}, (18), compute E(e_{t}) and Var(e_{t}).

• Compute E(e_{t}). Will use, Theorem 2.15 (ii), using the basic properities of Itô’s integrals.

Recall: E [∫_{0}^{t} f(s, ·)dW_{s}] = 0

E(e_{t}) = E (∫_{0}^{t} be^{a(t−s)}dW_{s})
= 0 by Thm 2.15  (19)

• Compute V ar(e_{t}).

V ar(e_{t}) = E[e_{t}^{2}] − [\underbrace{E[e_{t}])}_{see (19)}]^{2}

= E[e_{t}^{2}]
= E [(∫_{0}^{t} be^{a(t−s)}dW_{s})^{2}]  (20)

Will use, Itô isometry in order to simplify (20).

Recall: E [∫_{0}^{t} f(W_{s}, s)dW_{s}]^{2} = ∫_{0}^{t} Ef^{2}(W_{s}, s)ds

From above, (20) will become,

V ar(e_{t}) = E [(∫_{0}^{t} be^{a(t−s)}dW_{s})^{2}]

= E (∫_{0}^{t} b^{2}e^{2a(t−s)}ds)
= \frac{b^{2}}{2a} (e^{2at} − 1)

b. Assume now that |a(x) − a(y)| ≤ C_{a}|x − y| with a positive constant C_{a}. Find bounds for the expectation, E(|e_{t}|^{2}) and use it to bound the variance Var(e_{t}). Discuss what happens as b → 0.

• Find bounds for the expectation, E(|e_{t}|^{2}).

e_{t} = X_{t} − Z_{t}
= ∫_{0}^{t}(a(X_{s}) − a(Z_{s}))ds + ∫_{0}^{t} bdW_{s}  (21)

Taking the absolute value, squaring both sides, and taking the expectation we will obtain the left hand solution. To obtain a upper bound, we will use various methods involving the triangle inequality, Cauchy-Schwarz inequality, Fubini’s theorem, and Grönwall’s lemma.

|e_{t}| = |∫_{0}^{t}(a(X_{s}) − a(Z_{s}))ds + ∫_{0}^{t} b dW_{s}|
≤ ∫_{0}^{t} C_{a}|X_{s} − Z_{s}|ds + b |W_{t}| (b > 0, by assumption)
= ∫_{0}^{t} C_{a}|e_{s}|ds + b |W_{t}|
|e_{t}| ≤ ∫_{0}^{t} C_{a}|e_{s}|ds + b |W_{t}|  (22)

From (22), we will square both sides, and take the expectation to obtain:

E|e_{t}|^{2} ≤ E [(∫_{0}^{t} C_{a}|e_{s}|ds + b |W_{t}|)^{2}]

≤ 2E [(∫_{0}^{t} C_{a}|e_{s}|ds)^{2}] + 2E [(b |W_{t}|)^{2}] Recall: (a + b)^{2} ≤ 2(a^{2} + b^{2})

= 2E [(∫_{0}^{t} C_{a}|e_{s}|ds)^{2}] + 2E [b^{2} W_{t}^{2}]

≤ 2E [(∫_{0}^{t} C_{a}^{2}ds) (∫_{0}^{t} |e_{s}|^{2}ds)] + 2b^{2}t Recall: Cauchy-Schwarz inequality
≤ 2C_{a}^{2}tE [∫_{0}^{t} |e_{s}|^{2}ds] + 2b^{2}t
E|e_{t}|^{2} ≤ 2C_{a}^{2}t [∫_{0}^{t} E|e_{s}|^{2}ds] + 2b^{2}t Recall: Fubini’s theorem  (23)

To simplify the integral in (23), we will use Grömwall’s lemma.

Lemma. Assume that there exist positive constants A and K such that the function f : ℝ → ℝ satisfies:

f(t) ≤ K ∫_{0}^{t} f(s)ds + A

then

f(t) ≤ Ae^{Kt}  (24)

Also, to use Grönwall’s lemma, define the following:

f(t) = E|e_{t}|^{2}  (25)

From (25), and substituting into (23), we will obtain:

f(t) ≤ 2C_{a}^{2}t ∫_{0}^{t} f(s)ds + 2b^{2}t,
≤ 2C_{a}^{2}T ∫_{0}^{t} f(s)ds + 2b^{2}T  (26)

Let K = 2C_{a}^{2}T and A = 2b^{2}T, (note: A > 0) Then by Grönwall’s lemma, obtain a bound on E(|e_{t}|^{2}).

f(t) ≤ 2b^{2}T e^{2C_{a}^{2}T t}  (27)

• Using answer from part 2a), (27), obtain a bound on the variance Var(e_{t}).

Var(e_{t}) = E(e_{t}^{2}) − \underbrace{[E(e_{t})}_{see (19)}]^{2}

≤ E(e_{t}^{2}) (Note : [E(e_{t})]^{2} ≥ 0)
Var(e_{t}) ≤ 2b^{2}T e^{2C_{a}^{2}T t}  (28)

• Discuss what happens as b → 0.

From (28), as b → 0, will lead to:

0 ≤ Var(e_{t}) ≤ 0 ⟹ Var(e_{t}) = 0  (29)

Therefore, the stochastic term disappears as b → 0 in the SDE of X_{t}. Thus, matching the deterministic function, Z_{t} c) Implement a uniform time step forward Euler discretization of the above equations taking

a(x) = cos(x), b = 0.1, and T = 6. Plot the sample estimator for Var(e_{t}) vs. time. Compare it with the bound obtained in part (b). Use M = 10^{3} sample paths and different number of time steps: N = 10, 20, 40.

• Each graph below was graphed using Matlab, please see appendix for code.

• For each figure, the deterministic function, trajectories, and the bounds were plotted.

• Figure one: Number of time steps: N = 10.

width 20cm n10.png

• Figure two: Number of time steps: N = 20. width 20cm n20-2.png

• Figure three: Number of time steps: N = 40. width 20cm n40.png

APPENDIX: Problem 2(c) Mathlab code

clf randn(round(clock),1); T = 6; N = 40; Delta_t = T/N; b=0.1; X0=0;
C1=2*(b^2); C2=8; v_error_mean=zeros(N,1); v_error_square=zeros(N,1);
v_variance=zeros(N,1); v_xdet=zeros(N,1); v_xdet2=zeros(N,1); M = 1e+3;
for s = 1:M
X=X0;
Z=X0;
for j = 1:N
Winc1 = sqrt(Delta_t)*randn;
X=X+cos(X)*Delta_t + b*Winc1;
Z=Z+cos(Z)*Delta_t;
v_xdet(j,1)=v_xdet(j,1)+X;
v_xdet2(j,1)=v_xdet2(j,1)+X^2;
e=abs(X-Z);
v_error_mean(j,1)=v_error_mean(j,1)+e;
e2=(X-Z)^2;
v_error_square(j,1)=v_error_square(j,1)+e2;
end end X=X0; X1=X0; Z=X0; for j=1:N
v_xdet(j,1)=v_xdet(j,1)/M;
v_xdet2(j,1)=v_xdet2(j,1)/M;
v_variance_det(j,1)=v_xdet2(j,1)-((v_xdet(j,1))^2);
Winc1 = sqrt(Delta_t)*randn;
X=X+cos(X)*Delta_t + b*Winc1;
v_notmean(j,1) = X;
Winc2 = sqrt(Delta_t)*randn;
X1=X1+cos(X1)*Delta_t + b*Winc2;
v_notmean2(j,1) = X1;
v_error_mean(j,1)=v_error_mean(j,1)/M;
v_error_square(j,1)=v_error_square(j,1)/M;
v_variance(j,1)=v_error_square(j,1)-((v_error_mean(j,1))^2);
Z=Z+cos(Z)*Delta_t;
v_deterministic(j,1)=Z;
v_time=linspace(0,T,N);
v_minus_2d(j,1)=v_xdet(j,1)-2*sqrt(v_variance_det(j,1));
v_plus_2d(j,1)=v_xdet(j,1)+2*sqrt(v_variance_det(j,1)); end
title('N = 10'); subplot(2,1,1); plot(v_time, v_variance,'r');
legend('Variance(e(t))'); xlabel('Time (in seconds)'); subplot(2,1,2);
plot(v_time, v_deterministic, 'k', v_time, v_xdet,'b', v_time, v_notmean,
'm+--', v_time, v_notmean2,'gx--'); hold on; plot(v_time, v_minus_2d,'r',
v_time, v_plus_2d,'r','LineWidth',3); legend('Z','Avg(X)','Trajectory
1','Trajectory 2','Avg(X) +/- 2 SD');
xlabel('Time (in seconds)');