Stochastic Differential Equations: Lab 1



In [2]:

    
from IPython.core.display import HTML
css_file = '../ipython_notebook_styles/ngcmstyle.css'
HTML(open(css_file, "r").read())









    Out[2]:

This background for these exercises is article of D Higham, An Algorithmic Introduction to Numerical Simulation of Stochastic Differential Equations, SIAM Review 43:525-546 (2001). Higham provides Matlab codes illustrating the basic ideas at http://personal.strath.ac.uk/d.j.higham/algfiles.html, which are also given in the paper.

For random processes in python you should look at the numpy.random module. To set the initial seed (which you should not do in a real simulation, but allows for reproducible testing), see numpy.random.seed.

Brownian processes

Simulate a Brownian process over $[0, 1]$ using a step length $\delta t = 1/N$ for $N = 500, 1000, 2000$. Use a fixed seed of 100. Compare the results.



In [ ]:

Evaluate the function $u(B(t)) = \sin^2(t + B(t))$, where $B(t)$ is a Brownian process, on $M$ Brownian paths for $M = 500, 1000, 2000$. Compare the average path for each $M$.



In [ ]:

Stochastic integrals

Write functions to compute the Itô and Stratonovich integrals of a function $h(t, B(t))$ of a given Brownian process $B(t)$ over the interval $[0, 1]$.



In [ ]:

Test the functions on $h = B(t)$ for various $N$. Compare the limiting values of the integrals.



In [ ]:

Euler-Maruyama's method

Apply the Euler-Maruyama method to the stochastic differential equation

$$ \begin{equation} dX(t) = \lambda X(t) + \mu X(t) \ dB(t), \qquad X(0) = X_0. \end{equation} $$

Choose any reasonable values of the free parameters $\lambda, \mu, X_0$.

The exact solution to this equation is $X(t) = X(0) \exp \left[ \left( \lambda - \tfrac{1}{2} \mu^2 \right) t + \mu B(t) \right]$. Fix the timetstep and compare your solution to the exact solution.



In [ ]:

Vary the timestep of the Brownian path and check how the numerical solution compares to the exact solution.

Convergence

Investigate the weak and strong convergence of your method, applied to the problem above.



In [ ]:

Milstein's method

Implement Milstein's method, applied to the problem above.



In [ ]:

Check the convergence again.



In [ ]:

Compare the performance of the Euler-Maruyama and Milstein method using eg timeit. At what point is one method better than the other?



In [ ]:

Population problem

Apply the algorithms, convergence and performance tests to the SDE

$$ \begin{equation} dX(t) = r X(t) (K - X(t)) \ dt + \beta X(t) \ dB(t), \qquad X(0) = X_0. \end{equation} $$

Use the parameters $r = 2, K = 1, \beta = 0.25, X_0 = 0.5$.



In [ ]:

Investigate how the behaviour varies as you change the parameters $r, K, \beta$.



In [ ]: