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

Exercise 1 The following stochastic volatility model [FPSS00] generalizes the well known Black-Scholes geometric Brownian motion model [BS73], improving some aspects of option pricing. A simplified version of the model reads

dS(t) = rS(t)dt + e^{Y (t)}S(t)dW (t)  (1)
dY (t) = (−α(1 + Y (t)) + 0.4\sqrt{α}\sqrt{1 − ρ^{2}}) dt + 0.4\sqrt{α}dẐ(t),  (2)

where W and Z are independent Wiener process and

Ẑ(t) ≡ ρ W (t) + \sqrt{1 − ρ^{2}}Z(t)  (3)

Note: The correlation coefficient is ρ = −0.3.

To solve the SDE (1), define the following:

f(x) = log(x)  (4)

Solving the SDE (1) by Itô formula:

d log S_{t} = \frac{1}{S_{t}} dS_{t} − \frac{1}{2S_{t}^{2}} (dS_{t})^{2}
= r dt + e^{Y_{t}}d W_{t} − \frac{1}{2}e^{2Y_{t}}d_{t}
= (r − \frac{1}{2}e^{2Y_{t}}) dt + e^{Y_{t}}d W_{t}  (5)

Integrating (5), obtain the following:

log S_{t} = log S_{0} + ∫_{0}^{t}(r − \frac{1}{2}e^{2Y_{s}}) ds + ∫_{0}^{t} e^{Y_{s}}d W_{s}
S_{t} = S_{0}e^{∫_{0}^{t}(r− \frac{1}{2} e^{2Ys})ds+∫_{o}^{t} e^{Ys}d W_{s}}  (6)

Solving the SDE (2):

dY_{t} = (−α(1 + Y_{t}) + 0.4\sqrt{α}\sqrt{1 − ρ^{2}}) dt + 0.4\sqrt{α}dẐ_{t}

= (−α Y_{t} + \underbrace{0.4\sqrt{α}\sqrt{1 − ρ^{2}} − α}_{=M}) dt + \underbrace{0.4\sqrt{α}}_{σ} dẐ_{t}
= (−αY_{t} + M)dt + σ dẐ_{t}
= −α (Y_{t} − \frac{M}{α} ) dt + σ dẐ_{t} α ≠ 0  (7)

Define the following:

X_{t} = Y_{t} − \frac{M}{α}  (8)

Rewriting (7) by substituting (8), obtain the following:

dX_{t} = −α X_{t} dt + σ dẐ_{t}  (9)

To solve (9) use the following integrating factor:

V_{t} = e^{α t}X_{t}  (10)

Differentiating (10) and subsituting (9) obtain:

dV_{t} = d(e^{α t}X_{t})
= αe^{α t}X_{t}dt + e^{α t} \underbrace{dX_{t}}_{=(9)}

= αe^{α t}X_{t}dt − e^{α t}α X_{t}dt + e^{α t}σ dẐ_{t}
dV_{t} = e^{α t}σ dẐ_{t}  (11)

From (11), as well as using (10), obtain an equation for X_{t}.

\underbrace{V_{t}}_{(10)} = V_{0} + ∫_{0}^{s} e^{α s}σ dẐ_{s}

e^{α t}X_{t} = X_{0} + ∫_{0}^{t} e^{α s}σ dẐ_{s}
X_{t} = X_{0}e^{−α t} + ∫_{0}^{t} e^{−α (t−s)}σ dẐ_{s}  (12)

From (8), and substituting into above, we obtain an equation for Y_{t}.

\underbrace{X_{t}}_{(12)} = X_{0}e^{−α t} + ∫_{0}^{t} e^{−α (t−s)}σ dẐ_{s}

Y_{t} − \frac{M}{α} = (Y_{0} − \frac{M}{α} ) e^{−α t} + ∫_{0}^{t} e^{−α (t−s)}σ dẐ_{s}
Y_{t} = (Y_{0} − \frac{M}{α} ) e^{−α t} + \frac{M}{α} + ∫_{0}^{t} e^{−α (t−s)}σ dẐ_{s}  (13)

Compare the following explicit and implicit method for the computation of the option value below. Explicit Method:

S_{n+1} − S_{n} = rS_{n} Δ t + e^{Y_{n}}S_{n} Δ W_{n},
Y_{n+1} − Y_{n} = (−α(1 + \underbrace{Y_{n}}) + 0.4\sqrt{α}\sqrt{1 − ρ^{2}}) Δ t + 0.4\sqrt{α} Δ Ẑ_{n}
Ẑ_{n} = ρ W_{n} + \sqrt{1 − ρ^{2}}Z_{n}

Inside the code: Explicit method

    S1(n+1) = (1+ dt*r + dW(n)*exp(V(n)))*S1(n);
    V(n+1) = V(n)+(-alpha*V(n)-alpha+tmp*sqrt(1-rho^2))*dt + tmp * dhZ(n);

Implicit Method:

S_{n+1} − S_{n} = rS_{n} Δ t + e^{Y_{n}}S_{n} Δ W_{n},
Y_{n+1} − Y_{n} = (−α(1 + \underbrace{Y_{n+1}}) + 0.4\sqrt{α}\sqrt{1 − ρ^{2}}) Δ t + 0.4\sqrt{α} Δ Ẑ_{n}
Ẑ_{n} = ρ W_{n} + \sqrt{1 − ρ^{2}}Z_{n}

Inside the code: Implicit method

    S(n+1) = (1+ dt*r + dW(n)*exp(Y(n)))*S(n);
    Y(n+1) = (Y(n) + (-alpha+tmp*sqrt(1-rho^2))*dt + tmp * dhZ(n))/(1+alpha*dt);

Recall: The implicit Euler is absolutely stable; the explicit Euler method is not. Computation of the option value. Use α = 200, r = 0.04, T = \frac{3}{4} and S_{0} = K = 100.

V_{T} = e^{−rT} E[max(S_{T} − K, 0)]  (14)

a. Motivate which above numerical methods is best in this case. Provide a numerical approximation, based on uniform time steps, to the option value with accuracy TOL=10^{−2}.

Setup for numerical approximation:

\underline{Pseudocode} Selection of Δt and M

• choose Δt, fix M = 10^{3}

while E[g(X_{Δt} − g(X_{\frac{Δt}{2}} )] + \frac{Cσ_{error}}{\sqrt{M}} > \frac{1}{3}TOL

output: Δt = Δt^{∗}

• choose M, fix Δt^{∗} while σ̂^{2} = E(g(X)^{2}) − E(g(X))^{2}

check: |error_{statistical}| < \frac{C_{α}σ̂}{\sqrt{M}} ≤ \frac{2}{3}TOL

b. In the limit α → ∞ one obtains the results from geometric Brownian motion. Can you verify this by numerical experiments?

Note: In geometric Brownian motion with constant volatity, we can compute the value of the option exactly using the the Black-Scholes formula. The numerical results obtained approach this thoerotical price.

From (13), let the limit α → ∞,

Y_{t} = (Y_{0} − \frac{M}{α} ) e^{−α t} + \frac{M}{α} + ∫_{0}^{s} e^{−α (t−s)}σ dẐ_{s}

lim_{α→∞} \frac{Y_{0}}{e^{α t}} = 0

lim_{α→∞} −\frac{0.4\sqrt{α}\sqrt{1 − ρ^{2}} + α}{α e^{α t}} = 0

lim_{α→∞} \frac{0.4\sqrt{α}\sqrt{1 − ρ^{2}} − α}{α} = −1  (15)

lim_{α→∞} −\frac{0.4\sqrt{α}}{e^{α (t−s)}} = 0

From above, (15), we obtain as α → ∞, Y_{t} → −1.

However, when we calculate the V ar(Y_{t}), we see that similiar behavior as above (15), does not occur:

E [∫_{0}^{t} \sqrt{α} e^{−α(t−s)}dẐ_{s}]^{2} = ∫_{0}^{t} αe^{−2α(t−s)}ds
= \frac{e^{−2α(t−s)}}{2α} α|_{0}^{t}
= \frac{1}{2} − \frac{e^{−2αt}}{2}  (16)

Therefore, from(16), let the limit α → ∞, obtain:

lim_{α→∞} E [∫_{0}^{t} \sqrt{α} e^{−α(t−s)}dẐ_{s}]^{2} = lim_{α→∞} (\frac{1}{2} − \underbrace{\frac{e^{−2αt}}{2}}_{=0})

= \frac{1}{2}  (17)

Using this results, we can see that the volatilty from (1):

e^{Y_{t}} → e^{−1} as α → ∞

From (6), as limit α → ∞

S_{t} = S_{0}e^{∫_{0}^{t}(r− \frac{1}{2} e^{2Ys})ds+∫_{o}^{t} e^{Ys}d W_{s}}

= S_{0} exp \{∫_{0}^{t} (r − \frac{1}{2e^{2}} ) ds + ∫_{0}^{t} \frac{1}{e} dW_{s}\}
= S_{0} exp \{(r − \frac{1}{2e^{2}} ) t + \frac{1}{e}W_{t}\}  (18)

• Figure one: Histogram: Option Prices width 20cm sde_{h}w7_{h}ist.png

• Figure two: Numerical Solutions width 17cm sde_{h}w7_{1}.png

• Figure three: Convergence of option price width 17cm sde_{h}w7_{2}.png

• Figure four: Time discretization error width 20cm timesteps.png