If we wanted to simulate the SDE

\[dX_t = b(X_t) dt + \sigma(X_t) dW_t\]

on a computer then we truly want to approximate the associated integral equation over a sort time time interval of length \(h\). Namely,

\[X_{t+h} – X_t= \int_t^{t+h}b(X_s) ds + \int_t^{t+h}\sigma(X_s) dW_s\]

It is reasonable for \(s \in [t,t+h]\) to use the approximation

\begin{align}b(X_s) \approx b(X_t) \qquad\text{and}\qquad\sigma(X_s)\approx \sigma(X_t)\end{align}

which implies

\[X_{t+h} – X_t\approx b(X_t) h + \sigma(X_t) \big( W_{t+h}-W_{t}\big)\]

Since \(W_{t+h}-W_{t}\) is a Gaussian random variable with mean zero and variance \(h\), this discussion suggests the following numerical scheme:

\[ X_{n+1} = X_n + b(X_n) h + \sigma(X_n) \sqrt{h} \eta_n\]

where \(h\) is the time step and the \(\{ \eta_n : n=0,\dots\}\) are a collection of mutually independent standard Gaussian random variable (i.e. mean zero and variance 1). This is called the Euler–Maruyama method.

Use this method to numerically approximate (and plot) several trajectories of the following SDEs numerically.

- \[dX_t = -X_t dt + dW_t\]
- For \(r=-1, 0, 1/4, 1/2, 1\)

\[ dX_t = r X_t dt + X_t dW_t\]Does it head to infinity or approach zero ?Look at different small \(h\). Does the solution go negative ? Should it ? - \[dX_t = X dt -X_t^3 dt + \alpha dW_t\]

Try different values of \(\alpha\). For example, \(\alpha = 1/10, 1, 2 ,10\). How does it act ? What is the long time behavior ? How does it compare with what is learned in the “One dimensional stationary measure” problem about the stationary measure of this equation.