# Stochastic Volatility (Euler's Scheme)
stoch_vol_Euler <- function(S0, V0, t0, T, delta_t, r, kappa, theta, eta, rho) {
time <- seq(t0, T, by=delta_t)
S <- c(S0, rep(NA, T/delta_t))
V <- c(V0, rep(NA, T/delta_t))
W1 <- rnorm(T/delta_t, 0, delta_t)
W2 <- rnorm(T/delta_t, 0, delta_t)
for(n in 1:(length(time)-1)){
V[n+1] = V[n] + kappa * (theta - V[n]) * delta_t + eta * sqrt(V[n]) * (rho * W1[n] + sqrt(1-rho^2)*W2[n])
S[n+1] = S[n] + r * S[n] + S[n]*sqrt(V[n])*W1[n]
}
results <- data.frame(
Price = S,
Volatility = V
)
return(results)
}1 SDE Simulation
In Financial Mathematics when pricing derivatives we often wish to simulate continuous time stochastic price processes to in turn simulate derivative payoffs and hence estimate option prices using SDE Monte-Carlo Methods. In this note we will detail some of the methodologies employed to perform these simulations as well as quantifying their error and work requirements.
Consider the general SDE with differential form \[ d X_{t} = f(t, X_{t})dt + g(t, X_{t})dW_{t}, \] with Lipschitz \(f,g\) and Brownian Motion \(W_{t}\). The integral form is the exact solution \[ X_{t}=X_{t_{0}}+\int_{t_{0}}^{t}{f(s, X_{s})}~d{s}+ \int_{t_{0}}^{t}{g(s,X_{s})}~d{W_{s}}. \] Supposing that no explicit solution exists we pose the problem of how can we discretize and simulate the two integrals on the RHS of the expression?
We introduce two methodologies for simulating this process, first the simple but flawed Euler-Maruyama Method followed by the more robust and commonly used Milstein Method. We define the weak and strong errors of convergence of SDE simulation and compute both error orders for each method. We include several key examples including applications to option pricing, jump processes, multiple correlated Brownian drivers and Brownian Bridges.
1.1 Euler-Maruyama Method
The Euler-Maruyama Method is a SDE Simulation method with dynamics given by stochastic differential equations. The Itô Formula for a scalar-valued function \(U(t, X_{t})\) of the solution \(X_{t}\) of the scalar Itô SDE has the integral representation \[ U(t,X_{t}) = U(t_{0}, X_{t_{0}})+ \int_{t_{0}}^{t}{L^0U(s,X_{s})}~d{s} + \int_{t_{0}}^{t}{L^1U(s, X_{s})}~d{W_{s}}, \] where the differential operators \(L^0\) and \(L^1\) are defined by \[ L^0 = \frac{\partial}{\partial t}+f \frac{\partial}{\partial x}+ \frac{1}{2}g^2 \frac{\partial^2}{\partial x^2}\quad \& \quad L^1 = g \frac{\partial}{\partial x}. \] This differs from the deterministic chain rule by the additional third term in the \(L^0\) operator, which is essentially due to the fact that \(\mathbb{E}[\Delta W^2]=\Delta t\). When \(U(t,x)=x\) the Itô formula above is simply \[ X_{t}=X_{t_{0}}+\int_{t_{0}}^{t}{f(s,X_{s})}~d{s} + \int_{t_{0}}^{t}{g(s,X_{s})}~d{}W_{s}. \] Applying the Itô formula to \(f\) and \(g\) we obtain \[ \begin{align} X_{t} & =X_{0} \\ & \qquad+ \int_{t_{0}}^{t}{\left[ f(t_{0}X_{t_{0}})+\int_{t_{0}}^{s}{L^0f(u,X_{u})}~d{u}+\int_{t_{0}}^{s}{L^1 f(u, X_{u})}~d{W_{u}} \right]}~d{s} \\ & \qquad +\int_{t_{0}}^{t}{\left[ g(t_{0}, X_{t_{0}})+\int_{t_{0}}^{s}{L^0g(u, X_{u})}~d{u}+\int_{t_{0}}^{s}{L^1g(u, X_{u})}~d{W_{u}} \right]}~d{W_{s}} \\ & =X_{t_{0}}+f(t_{0}, X_{t_{0}})\int_{t_{0}}^{t}{}~d{s}+g(t_{0}, X_{t_{0}})\int_{t_{0}}^{t}{}~d{W_{s}}+R_{1}(t,t_{0}), \end{align} \] where the remainder is \[ \begin{align} R_{1}(t,t_{0}) & =\int_{t_{0}}^{s}{\int_{t_{0}}^{s}{L^0f(u, X_{u})}~d{u}}~d{s} + \int_{t_{0}}^{s}{\int_{t_{0}}^{s}{L^1f(u, X_{u})}~d{W_{u}}}~d{s} \\ & + \int_{t_{0}}^{t}{\int_{t_{0}}^{s}{L^0g(u, X_{u})}~d{u}}~d{W_{s}} + \int_{t_{0}}^{t}{\int_{t_{0}}^{s}{L^1 g(u, X_{u})}~d{W_{u}}}~d{W_{s}}. \end{align} \] Replacing \(t_{0}\) by \(t_{n}\), \(t\) by \(t_{n+1}\) and discarding the remainder be obtain the Euler-Maruyama Method which we summarize in the definition below. This is the simplest stochastic Taylor expansion.
The Euler-Maruyama Method discretizes the continuous time dynamics by defining a step-size \(\Delta t=\frac{T}{N}\) for \(N \in \mathbb{Z}\) and computing the approximate solutions at times \(t_{i}=i\Delta t\) using simple left-endpoint Riemann sums in place of the integrals \[ \widehat{X}_{t_{n+1}}=\widehat{X}_{t_{n}}+\Delta tf(t_{n},\widehat{X}_{t_{n}})+\underbrace{\Delta W_{n}}_{=\sqrt{\Delta t}\cdot\xi_{n}}g(t_{n},\widehat{X}_{t_{n}}), \] where \(\Delta W_{n}=W(t_{n+1})-W(t_{n})\sim\mathcal{N}(0,\Delta t)\) (or equivalently with \(\xi_{n}\sim\mathcal{N}(0,1)\) for all \(n\)).
Example: BS Model EM Method and Error Consider the Black-Scholes-Merton Model which consists of two assets, a risk free asset \(B\) and a stock \(S\) with price dynamics given by \[
\begin{align}dB_{t}&=rB_{t}dt \\dS_{t}& =\alpha S_{t}dt+\sigma S_{t}dW_{t},\end{align}
\] on filtered probability space \((\Omega, \mathcal{F}, \mathbb{F}, \mathbb{P})\) where \(r, \alpha\) and \(\sigma\) are deterministic constants and \(W_{t}\) is BM in \(\mathbb{P}\). Under the risk-neutral measure \(\mathbb{Q}\) the dynamics of \(S_t\) are instead given by \[
dS_{t}=rS_{t}dt+\sigma S_{t}dW_{t}^\mathbb{Q},
\] and the true solution under the risk-free measure is \[
S_{T}=S_{0}e^{(r-\frac{\sigma^2}{2})T+\sigma W_{T}}.
\] for some future time \(T\). Suppose we wish to simulate this process on the interval \([0,T]\). We can apply Euler’s algorithm by considering a discretized mesh \(0=t_{0}, t_{1}, \dots, t_{n}=T\) where \(t_{n}=t_{n-1}+\Delta t\) for \(\Delta t=\frac{T}{n}\) and then writing \[
S_{t+\Delta t}=S_{t}+rS_{t}\Delta t + \sigma S_{t}\underbrace{(W_{t+\Delta t}-W_{t})}_{=\sqrt{ \Delta t }X_{t}},
\] where \(X_{t}\stackrel{iid}\sim\mathcal{N}(0,1)\). This is easy to simulate, for example the following R code using a for loop:
## Euler's Algorithm
bs_euler <- function(r, alpha, sigma, S0, t0=0, T, n=1000){
dt = T/n
S <- numeric(n)
S[1] <- S0
for(i in 2:n){
S[i] <- S[i-1] + r*S[i-1]*dt + sigma*S[i-1]*sqrt(dt)*rnorm(1,0,1)
}
return(S)
}We have clearly made some kind of discretization error. To determine how bad this error is we turn back to the true solution and we consider the Taylor expansion \[ \begin{align} S_{t+\Delta t} & =S_{t}e^{\left( r-\frac{1}{2}\sigma^2 \right)\Delta t+\sigma \sqrt{ \Delta t }\xi} \\ & =S_{t}\left\{ \underbrace{1+\left(r-\frac{1}{2}\sigma^2\right)\Delta t+\sigma \sqrt{ \Delta t }\xi}_{\text{Euler's Approximation}}+\frac{1}{2}\left[ \left( r- \frac{1}{2}\sigma^2\right)\Delta t+\sigma \sqrt{ \Delta t }\xi \right]^2 + \dots \right\} \\ & =S_{t}\left\{ 1+\sigma(\Delta t)^{\frac{1}{2}}\xi+\Delta t \left[ \left( r- \frac{1}{2}\sigma^2 \right)+\frac{1}{2}\sigma^2\xi^2 \right]+(\Delta t)^{\frac{3}{2}}\left[ \dots \right]+\dots \right\}, \end{align} \] where we have replaced \(W_{t+\Delta t}-W_{t}\sim \sqrt{ \Delta t }\xi\) where \(\xi\sim\mathcal{N}(0,1)\). From this we see that Euler’s approximation might seem reasonable from the SDE perspective but it is actually rather strange in its use of various orders of \(\Delta t\).
It seems reasonable to suggest that the approximation improves as \(N \to \infty\). To test this presumption, we consider the two standard approaches for measuring the error \(X_{n}\to X(t_{n})\), leading to the concepts of weak error and strong error.
1.2 Milstein Method
The Milstein Method applies a similar methodology to the Euler-Maruyama Method but provides a greater degree of accuracy as it also includes additional terms from the remainder in the stochastic Taylor expansion above. We apply the Itô Formula to the integrand \(L^1g\) in the fourth double integral of the remainder \(R_{1}(t, t_{0})\) we obtain the stochastic Taylor expansion \[ \begin{align} X_{t}=X_{0}+f(t_{0}, X_{t_{0}})\int_{t_{0}}^{t}{}~d{s}+g(t_{0}, X_{t_{0}})\int_{t_{0}}^{t}{}~d{W_{s}}+L^1 & g(t_{0}, X_{t_{0}})\int_{{t_{0}}}^{t}{\int_{t_{0}}^{s}{}~d{W_{u}}}~d{W_{s}}+R_{2}(t, t_{0}), \end{align} \] where the remainder is \[ \begin{align} R_{2}(t, t_{0}) & =\int_{t_{0}}^{s}{\int_{t_{0}}^{s}{L^0f(u, X_{u})}~d{u}}~d{s} + \int_{t_{0}}^{s}{\int_{t_{0}}^{s}{L^1f(u, X(u))}~d{W_{u}}}~d{s}+\int_{t_{0}}^{t}{\int_{t_{0}}^{s}{L^0g(u, X_{u})}~d{u}}~d{W_{s}} \\ & \quad+ \int_{t_{0}}^{t}{\int_{t_{0}}^{s}{\int_{t_{0}}^{u}{L^0L^1g(v, X_{v})}~d{v}}~d{W_{u}}}~d{W_{s}}+\int_{t_{0}}^{t}{\int_{t_{0}}^{s}{\int_{t_{0}}^{u}{L^1L^1g(v, X_{v})}~d{W_{v}}}~d{W_{u}}}~d{W_{s}}. \end{align} \] Making the same replacements as above and once again discarding the remainder we obtain the Milstein Method which we summarize in the definition below.
The Milstein method is \[ X_{n+1}=X_{n}+\Delta t f(t_{n}, X_{n})+\Delta W_{n}g(t_{n}, X_{n})+L^1g(t_{n}, X_{t_{n}})\int_{t_{n}}^{t_{n+1}}{\int_{t_{n}}^{s}{}~d{W_{u}}}~d{W_{s}}, \] where the double integral is given by \[ \int_{t_{n}}^{t_{n+1}}{\int_{t_{n}}^{t}~d{W_{s}}}{}~d{W_{t}}=\frac{1}{2}((\Delta W_{n})^2-\Delta t). \]
Example: BM Model Milstein Algorithm An alternative approach could be to instead use the Milstein Algorithm \[
S_{t+\Delta t}=S_{t}\left\{ 1+\sigma(\Delta t)^{\frac{1}{2}}\xi+\Delta t \left[ r+ \frac{1}{2}\left( \sigma^2(\xi^2-1) \right) \right] \right\},\tag{6}
\] which has improved on Euler’s approximation by grouping terms based on powers of \(\Delta t\) and removed all terms with powers \(>1\). This is again very easy to simulate with a for loop:
## Milstein's Algorithm
S <- S0
for(i in 2:length(t)){
xi <- rnorm(1,0,1)
S[i] <- S[i-1]*(1+sigma*h^(1/2)*xi + h*(r+(1/2)*(sigma^2*(xi^2 - 1))))
}1.3 SDE Simulation Error
1.3.1 Weak Error
Given a function \(\Phi\) from some class (e.g. polynomials with degree \(d\leq k\)) we define the weak error by \[ e_{\Delta t}^\text{weak}:=\sup_{0\leq t_{n}\leq T}\lvert \mathbb{E}[\Phi(X_{n})]-\mathbb{E}[\Phi(X(t_{n}))] \rvert . \] The weak error measures how well the given methodology can approximate the mean of \(\Phi(X(t))\) (weak convergence measures the error of the means). We say that the method converges weakly if, for any \(\Phi\) in this class \[ e_{\Delta t}^\text{weak}\to 0~~\text{as}~~\Delta t \to 0. \] Further, we say that a method has weak order of convergence \(p\) if there exists a constant \(K\) and a step size level \(\Delta t^*\) (both dependent on \(\Phi\)) such that \[ e_{\Delta t}^\text{weak}\leq K\Delta t^p;~~\forall~ 0 < \Delta t\leq \Delta t^*, \] more specifically, we take \(p\) to be the largest value for which this holds.
For a wide class of problems the EM method converges with weak order one.
Proof Outline: Update from Higham and Kloeden §9.3.
It turns out that it is not necessary to compute exact Brownian paths to obtain weak order equal to one. It is sufficient for the increments \(\Delta W_{n}\) to have the correct low-order moments, hence we may us the weak EM method \[ X_{n+1}=X_{n}+\Delta tf(X_{n})+\overline{\Delta W_{n}} g(X_{n}), \] where \(\overline{\Delta W_{n}}\) are i.i.d. uniform random variables taking value \(\pm \sqrt{ \Delta t }\).
1.3.2 Strong Error
At each time \(t_{n}\) the absolute value of the error \(\lvert X_{n}-X(t_{n}) \rvert\) is a random variable. The strong error is defined as \[ e_{\Delta t}^\text{strong}:=\sup_{0\leq t_{n}\leq T}\mathbb{E}[\lvert X_{n}-X(t_{n}) \rvert ], \] i.e. the path-wise error averaged over all paths. Further, we say that a method converges strongly if \[ e_{\Delta t}^\text{strong}\to 0 ~~\text{as}~~ \Delta t \to 0, \] and we say that a method has strong order \(p\) on a class of SDEs if there exists \(K\) and \(\Delta t^*\) such that \[ e_{\Delta t}^\text{strong}\leq K \Delta t^p;~~\forall~0<\Delta t \leq\Delta t^*. \]
In general the EM method has strong order \(p=\frac{1}{2}\).
Proof: Update from proof outline given in Chapter 9.
The definition of strong convergence has important implications for individual paths, namely if we have that \[ \mathbb{E}\lvert X_{n}-X(t_{n}) \rvert \leq K\Delta t^{1/2}, \] then by the Markov Inequality we have that \[ \mathbb{P}(\lvert X_{n}-X(t_{n}) \rvert \geq \Delta t^{1/4})\leq K\Delta t^{1/4}, \] or in other words, along any path the error with be small with probability close to 1.
1.4 SDE Simulation Examples
1.4.1 Multiple Black-Scholes SDE Simulation
Suppose we wish to use the Euler scheme to simulate the following system of SDEs \[ \begin{cases} dS_{t}^1 = rS_{t}^1dt+\sigma_{1}S_{t}^1dW_{t}^1 \\ dS_{t}^2=rS_{t}^2dt+\sigma_{2}S_{t}^2dW_{t}^2, \end{cases} \] where we have that \(d\left< W^1, W^1 \right>_{t}=\rho~dt\). For discretized time step \(\Delta t\) we have \[ \begin{cases} \check{S}_{n+1}^1=\check{S}_{n}^1+r\check{S}_{n}^1\Delta t+\sigma_{1}\check{S}_{n}^1\Delta W_{n+1}^1 \\ \check{S}_{n+1}^2=\check{S}_{n}^2+r\check{S}_{n}^2\Delta t+\sigma_{1}\check{S}_{n}^2 \underbrace{\Delta W_{n+1}^2}_{\rho\Delta W_{n+1}^1+\sqrt{ 1-\rho^2 }\Delta \tilde{W}_{n+1}}, \\ \end{cases} \] where \(\tilde{W}_{n}\perp W_{n}^1\). The idea is that to avoid simulating the correlated Weiner processes we instead represent \(W_{n}^2\) in terms of \(W^1_{n}\) and some independent Weiner process \(\tilde{W}_{n}\). If we write the correlation matrix of \(W_{n}^1\) and \(W_{n}^2\) the Cholesky decomposition theorem states there exists a lower triangular matrix \(A\) such that \[ \begin{bmatrix} 1 & \rho \\ \rho & 1 \end{bmatrix} = AA^T=\begin{bmatrix} a_{11} & 0 \\ a_{12} & a_{22} \end{bmatrix} \begin{bmatrix} a_{11} & a_{12} \\ 0 & a_{22} \end{bmatrix}=\begin{bmatrix} a_{11}^2 & a_{11}a_{12} \\ a_{12}a_{11} & a_{11}^2 + a_{12}^2 \end{bmatrix}, \] which gives that \(a_{11}=1\), \(a_{12}=\rho\) and \(a_{22}=\sqrt{ 1-\rho^2 }\). Then we have that \[ \begin{bmatrix} W^1 \\ W^2 \end{bmatrix} \stackrel{\mathcal{L}}{=}A\cdot \begin{bmatrix} W^1 \\ \tilde{W} \end{bmatrix}, \] leading to our definition of our independent process and our new system \[ \begin{cases} \check{S}_{n+1}^1 = S_{n}^1 + r\check{S}_{n}^1\Delta t+\sigma_{1}\check{S}_{n}^1 \Delta W_{n+1}^1 \\ \check{S}_{n+1}^2=S_{n}^2+r\check{S}_{n}^2\Delta t +\sigma_{2}\check{S}_{n}^2\left( \rho \Delta W_{n+1}^1 + \sqrt{ 1-\rho^2 }\Delta \tilde{W}_{n+1} \right). \end{cases} \]
1.4.2 Heston Model SDE Simulation
A famous example of Multi-Dimenstional SDE simulation is the Heston Model for stochastic volatility whereby an asset price process \((S_{t})\) and corresponding volatility process \((V_{t})\) have dynamics under the risk-neutral measure given b \[ \begin{cases} dS_{t}=rS_{t}dt+S_{t}\sqrt{V_{t}}dW_{t}^1 \\ dV_{t}=\kappa(\theta-V_{t})dt+\eta \sqrt{ V_{t} }dW_{t}^2, \end{cases} \] where \(W^1, W^2\) are dependent Wiener processes with constant \(d\left< W^{1},W^2 \right>_{t}=\rho\). Suppose we wish to simulate these processes using the Euler scheme. Then for discretized step \(\Delta t\) we would simulate \[ \begin{cases} \check{S}_{n+1} = \check{S}_{n} + r\check{S}_{n}^1\Delta t+\check{S}_{n}\sqrt{ \check{V}_{n}} \Delta W_{n+1}^1 \\ \check{V}_{n+1}=\check{V}_{n}+\kappa (\theta-\check{V}_{n})\Delta t +\eta \sqrt{\check{V}_{n}}\left( \rho \Delta W_{n+1}^1 + \sqrt{ 1-\rho^2 }\Delta \tilde{W}_{n+1} \right). \end{cases} \] We can define the following function in R that simulates stochastic volatility via the Euler scheme:
We note that there are some subtle challenges to consider if we attempt the same applying the Milstein scheme and this is something we covered in assignment 1.
1.4.3 Jump Diffusion Processes SDE Simulation
The [[merton-jump-diffusion-model]] for stock price process \((S_{t})\) allows the process to jump by incorporating a Poisson process \((Nt)\) with rate \(\lambda\) that is independent of the driving standard Weiner process \((W_{t})\) \[ dS_{t}=rS_{t}dt+\sigma S_{t}dW_{t}+\eta S_{t}dN_{t}. \] First assume that the jump diffusion has homogeneous rate. To simulate such a process we would:
- Simulate the jump times; and then
- Simulate the process (e.g. with Euler or Milstein) between the jump points.
With a Poisson process times between jumps are exponentially distributed and so we can first simulate a collection of jump times that span over the time interval using exponential random variables and while loops.
We then need to discretize the interval \([t_{0},T]\) with time step \(\Delta t\) and proceed to apply Euler’s method between jumps (accounting for the fact that the jump points will intersect the discretization steps) and adding our jump impact at the correct points.
Below we include an example function generating a jump diffusion process with time homogeneous rate \(\lambda\):
# Jump Diffusion Process (Euler's Scheme)
jump_diff_Euler <- function(S0, t0, T, delta_t, r, sigma, eta, lambda){
# Simulate jumps on the interval
t <- rexp(1, lambda)
jump_times <- c()
while(t<T){
jump_times <- append(jump_times, t)
t <- t + rexp(1, lambda)
}
# Apply Eulers scheme between jumps
times <- append(seq(t0, T, delta_t), jump_times) |> sort()
S <- S0
for(n in 1:(length(times)-1)){
delta_times <- times[n+1] - times[n]
jump_component <- ifelse(times[n] %in% jump_times, 1, 0)
S[n+1] <- S[n] + r*S[n]*delta_times + sigma*S[n]*rnorm(1,0,delta_times) + eta*S[n]*jump_component
}
results <- data.frame(
Time = times,
Price = S
)
return(results)
}Next we consider the case of a process dependent rate \(\lambda(S_{t})\) which raises the question of whether we are still able pre-simulate the jump times. The answer is yes if there is some upper bound \[ \lambda(S_{t})\leq \overline{\lambda}, \] through the use of thinned Poisson processes. The idea is that we can simulate the jumps for this maximum rate \(\overline\lambda\) and then at each proposed jump point \(j_{i}\) we can accept or reject the joint with probabilities \[ f(\tau_{i})=\begin{cases} 1 & \text{w.p.} \quad \frac{\lambda (S_{t})\tau_{i}}{\overline\lambda} \\ 0 & \text{o.w.} \end{cases} \] ### Brownian Bridges SDE Simulation
The Brownian Bridge is an SDE \((X_{t})\) characterized by the dynamics \[ \begin{cases}dX_{t}=\frac{b-X_{t}}{1-t}dt+dB_{t}\\ X_{0}=a,\end{cases} \] for \(a,b \in \mathbb{R}\) and \(0 \leq t < 1\). The solution to this SDE is given by \[ X_{t}=a(1-t)+bt+(1-t)\int_{0}^t \frac{1}{1-s}dB_{s}. \] Suppose we have some skeletal process (e.g. end of day prices) and we wish to estimate intraday price movements assuming they have Brownian motion dynamics. In this case we can simulate a Brownian bridge (e.g. using Euler’s scheme) \[ \check{X}_{n+1} = \check{X}_{n}+\frac{b-\check{X}_{n}}{1-n}\Delta t + \Delta W_{n+1} \] where we specify \(\check{X}_0=a\), updating \(a\) and \(b\) for each new estimation interval.