How do I implement the Euler scheme for this SDE?

Question

I have the following problem:

Consider the system of (Ito) SDE's for the price $S_{t}$ of a stock, the stochastic and past-dependent volatility $\sigma_{z}$ and the long-term averaged volatility $\xi_{t}:$

\begin{array}{l} d S_{t}=\alpha S_{t} d t+\sigma_{i} S_{t} d w_{t}^{1} \\ d \sigma_{t}=-\left(\sigma_{t}-\xi_{t}\right) d t+p \sigma_{t} d w_{t}^{2} \\ d \xi_{t}=\frac{1}{\alpha}\left(\sigma_{z}-\xi_{i}\right) d t \end{array}with $S_{0}=€ 50, \quad \sigma_{0}=0.20, \quad \xi_{0}=0.20 \quad$ and $\alpha=0.10$

For $p=0($ and $\alpha \neq 0)$ the well-known Black-Scholes model is obtained: $$ d S_{t}=\alpha S_{i} d t+\sigma_{0} S_{i} d w_{i}^{1} $$

Question: How do I implement the Euler scheme for this SDE?

Nadine :)

Additional Question after answer of :

Unfortunately my knowledge is not sufficient to follow your answer. I found the following derivation on page 2 of: http://people.math.gatech.edu/~meyer/MA6635/chap4.pdf

We can now state Itô's lemma: Let $X$ satisfy $$ \begin{aligned} d X(t) &=a(X, t) d t+b(X, t) d W \\ X(0) &=X \end{aligned} $$

or in the case of my specific SDE:

$d S_{t}=\alpha S_{t} d t+\sigma_{i} S_{t} d w_{t}^{1} \\$

with $a(X,t) = \alpha S_{t}$ and $b(X,t) = \sigma_i S_t$

Then the following is stated:

"Assume that u(x, t) is a smooth function of the independent variables x and t. Then"

$d u(t)=\left[\frac{\partial u}{\partial t}+a(X, t) \frac{\partial u}{\partial x}+\frac{1}{2} \frac{\partial^{2} u}{\partial x^{2}} b(X, t)^{2}\right] d t+b(X, t) \frac{\partial u}{\partial x} d W$

I don't understand why only the $a(X,t)$ term gets replaced by this îto solution.

The second part states that:

If we apply Itô's lemma to $u=\ln s$ where $$ d S(t)=\mu S d t+\sigma S d W $$ then (with $X=S, a(S, t)=\mu S, b(S, t)=\sigma S)$ we find $$ \frac{\partial u}{\partial s}=\frac{1}{s}, \quad \frac{\partial^{2} u}{\partial s^{2}}=-\frac{1}{s^{2}}, \quad \frac{\partial u}{\partial t}=0 $$

I see where the derivatives come from, but I don't see why u = ln(s).

The document states that:

$d u(t)=\left[\mu-\sigma^{2} / 2\right] d t+\sigma d W$

But find that really strange since there are no $S$ terms or $a(X,t)$ or $b(X,t)$ terms are observed.

Thank you very much for your first reply, but it went way too quick. Since I can't go to univeristy due covid, this forum is the only forum of instruction I have.

Would you mind to elaborate the answers with more steps based on the additional questions?

Thankyou, Nadine :)

Answer 1

The Euler method is a rather crude method. The same goes for Euler-Maruyama. However, there is no trivial extension of ODE methods to SDE methods of the same order. The Milshtein method, even in its derivative-free variation, raises the strong approximation order to one. The next expansion terms of the exact solution contain integrals over the Brownian motion that can not be reduced to values at the sample points. This is an obstacle for higher strong order, it is slightly easier to construct methods with higher weak order, see Kloeden, Platen "Numerical Solution of Stochastic Differential Equations", Springer, 1992.

So it would be best to make the equation as simple as possible before starting with the numerics. Here one can do what is also done for the simple geometric Brownian motion and set $S_t=\exp(X_t)$ or $X_t=\log(S_t)$. Then the Ito formula gives $$ dX_t = \frac{dS_t}{S_t}-\frac12\frac{d\langle S_t\rangle}{S_t^2} =α dt+σ_tdW^1_t-\frac12σ_t^2dt $$ For the numerical step you now translate all differentials into differences and use their known properties $$\newcommand{\D}{\mathit{\Delta}} \begin{align} X_{t+\D t}&=X_t+(α-\tfrac12σ_t^2)\D t+σ_tz^1_t\sqrt{\D t}\\ σ_{t+\D t}&=σ_t−(σ_t−ξ_t)\D t+pσ_tz^2_t\sqrt{\D t}\\ ξ_{t+\D t}&=ξ_t+\frac1α(σ_t−ξ_t)\D t \end{align} $$ where $z^1_t$, $z^2_t$ are random numbers from a standard normal distribution. They get computed anew in every step. Alternatively, if for instance you want to compare the "strong" numerical convergence for different step sizes, you can generate the Brownian motions $W^k_t$ for the smallest step size and use the appropriate differences for larger step sizes. (To do that in the other direction "constructively", you would need to fill the gaps with Brownian bridges.)

t = linspace(0,T,N+1); dt = t[1]-t[0]
W = cumsum(randn([2,N+1])*dt**0.5)
# for an integration test with step size M*dt use
do_integration_test(t[::M],W[:,::M],...)

As the Brownian motion can be considered as approximately a random walk with very small step size, simplified numerical methods that might be somewhat faster are possible, for details and more see the always recommendable P. Forsyth: "An introduction to Computational Finance without Agonizing Pain".

---

As to the more general situation, when translating to the given equation you get $X_t=X(t)=S_t$ and the coefficients as you describe. If necessary, which it here is not, you could translate back, $u=\ln s\implies s=exp(u)$. With constant coefficients, there is no more to do, you can integrate directly $U_t=U_0+[μ−σ^2/2]t+σW_t$, so that $S_t=S_0\exp([μ−σ^2/2]t+σW_t)$, the equation for the geometric Brownian motion.