Use the left and right arrow keys to navigate the presentation forward and backward respectively. You can also use the arrows at the bottom right of the screen to navigate with a mouse.
FAIR USE ACT DISCLAIMER: This site is for educational purposes only. This website may contain copyrighted material, the use of which has not been specifically authorized by the copyright holders. The material is made available on this website as a way to advance teaching, and copyright-protected materials are used to the extent necessary to make this class function in a distance learning environment. The Fair Use Copyright Disclaimer is under section 107 of the Copyright Act of 1976, allowance is made for “fair use” for purposes such as criticism, comment, news reporting, teaching, scholarship, education and research.
In the last section, we considered the construction of a Gauss-Markov model in continuous time as generated from SDEs.
Although we focused on the scalar case, these concepts generalize naturally to systems of stochastic differential equations on random vectors.
Rather than belaboring the mathematical machinery, in the next two lectures we will consider how one simulates such systems of equations in practice.
We will begin again with the analogy of deterministic systems of ordinary differential equations (ODEs).
After we introduce some core techniques for ODEs, we will introduce how these methods are extended to systems that are simulated with governing laws with random shocks.
Recall the notion of the initial value problem, where
\[ \begin{align} \frac{\mathrm{d}}{\mathrm{d}t}\pmb{x}&: = \pmb{f}(t, \pmb{x}), \\ \pmb{x}(0)&:= \pmb{x}_0. \end{align} \]
We note that the Lipshitz condition discussed last time guarantees the existence and uniqueness of a solution.
Suppose that the solution is defined over an interval \( [0, T] \) where we make a partition as
\[ \begin{align} t_0 = 0 < t_1 < \cdots < t_n = T \end{align} \] where \( t_i - t_{i-1} = h \) for all \( i \).
Let's consider once again a tangent approximation to this solution, where for the perturbation \( t_i + h \) in time,
\[ \begin{align} \pmb{x}(t_{i+1})&= \pmb{x}(t_{i}) + \pmb{x}'(t_i) h +\mathcal{O}(h^2) \\ &=\pmb{x}(t_{i}) + f(t_i, \pmb{x}(t_i)) h +\mathcal{O}(h^2). \end{align} \]
Euler scheme
Let \( \pmb{x}(t) \) satisfy the initial value problem above; the approxmation of this solution by Euler’s scheme is defined recursively as \[ \begin{align} \hat{\pmb{x}}(t_{i+1}) &=\hat{\pmb{x}}(t_{i}) + f(t_i, \hat{\pmb{x}}(t_i)) h, \end{align} \] where \( h \) is known as the step size in time.
We note that Taylor's theorem guarantees that the difference of the approximation
\[ \begin{align}
\pmb{x}(t_{i+1})&=\pmb{x}(t_{i}) + f(t_i, \pmb{x}(t_i)) h +\mathcal{O}(h^2)
\end{align} \]
from the true solution shrinks at order \( \mathcal{O}(h^2) \) as \( h\rightarrow 0 \).
If \( \pmb{f} \) has continuous second derivatives, we can actually write
\[ \begin{align} \pmb{x}(t_{i+1})= &\pmb{x}(t_{i}) + \pmb{x}'(t_i) h + \frac{1}{2}\pmb{x}''(\tau_i) h^2 \end{align} \] for some \( \tau_i \in (t_i, t_{i+1}) \).
However, the repeated approximation using this rule accumulates error, due to the repeated miss-match between the numerical approximation and the true solution.
Consider, for the numerical solution defined again as \( \hat{\pmb{x}}(t) \),
\[ \begin{align} \pmb{x}(t_{i+1}) - \hat{\pmb{x}}(t_{i+1}) = \pmb{x}(t_{i}) - \hat{\pmb{x}}(t_{i}) + h\left[\pmb{f}(t_{i}, \pmb{x}(t_i)) - \pmb{f}(t_i, \hat{\pmb{x}}(t_i) \right]+ \frac{1}{2}\pmb{x}''(\tau_i) h^2 \end{align} \]
The terms \( \pmb{x}(t_{i}) - \hat{\pmb{x}}(t_{i}) + h\left[\pmb{f}(t_{i}, \pmb{x}(t_i)) - \pmb{f}(t_i, \hat{\pmb{x}}(t_i) \right] \) in the above difference are instead known as the propagation error.
Although the exact Taylor expansion means that the derivative approximation differs by terms at second order,
This motivates the following definition for the numerical solution of ODEs.
Numerical convergence of ODEs
Let \( \pmb{x}(t) \) be an exact solution to an initial value problem, and let \( \hat{\pmb{x}}(t) \) be a numerical approximation, using a maximum step size of \( h \). The approximation \( \hat{\pmb{x}} \) is said to converge to \( \pmb{x} \) at order \( \gamma \) if there exists a constant \( C \) independent of \( h \) such that \[ \begin{align} \max \parallel \pmb{x}(t) - \hat{\pmb{x}}(t) \parallel \leq C h^\gamma \end{align} \] for all \( t \in [0, T] \).
With the above definition in mind, it is a classic result that the Euler approximation converges at order \( 1.0 \) to the true solution to an initial value problem.
The Euler scheme is not used widely in practice due to this low order of convergence, but it is extremely useful to understand the principles of numerical simulation.
Particularly, higher-order Taylor expansions give better (higher-order) rates of convergence of the numerical solution to the true solution.
The Euler scheme also has a direct extension to SDEs that will likewise help explain the approximations made in numerical simulation of these random sample paths.
We recall the form of our generic SDE, but this time presented in vector notations,
\[ \begin{align} \mathrm{d}\pmb{x}:= \pmb{f}(t, \pmb{x}(t)) \mathrm{d}t + \mathbf{S}(t, \pmb{x}(t))\mathrm{d}\pmb{W}_t \end{align} \] where in the above:
The extension of Euler's scheme follows directly to the above by formally writing the following.
Euler-Maruyama
For the system of SDEs defined above, the numerical approximation of the sample path by Euler-Maruyama is defined as, \[ \begin{align} \hat{\pmb{x}}(t_{i+1}) := \hat{\pmb{x}}(t_{i}) + \pmb{f}(t_i, \hat{\pmb{x}}(t_i)) h + \mathbf{S}(t_i, \hat{\pmb{x}}) \boldsymbol{\xi} \sqrt{h} \end{align} \] where \( \boldsymbol{\xi} \sim N\left(\pmb{0}, \mathbf{I}_{N_x}\right) \).
Note that with \( \boldsymbol{\xi} \) distributed as the standard normal means that \( \sqrt{h}\boldsymbol{\xi} \sim N\left(\pmb{0}, h \mathbf{I}_{N_x}\right) \), giving a representation of the discretized Wiener process.
Intuitively, this gives the direct extension of the ODE system by analogy where we make a Gaussian perturbation to the mechanistic evolution of the state by the process laws.
However, as usual, convergence is a more subtle concept…
Strong convergence
Let \( \pmb{x}(t) \) be an exact sample path of a generic system of SDEs, and \( \hat{\pmb{x}}(t) \) be some numerical discretization with a maximum step size of \( h \) between time points. The approximation \( \hat{\pmb{x}} \) is said to strongly converge to \( \pmb{x} \) at order \( \gamma \) if there exists a constant \( C \) independent of \( h \) such that \[ \begin{align} \mathbb{E}\left[\parallel \pmb{x}(t) - \hat{\pmb{x}}(t) \parallel \right] \leq C h^\gamma \end{align} \] for all \( t \in [0, T] \).
Notice that this gives a direct analogy to the deterministic convergence;
For a specific sample path, given a particular realization of \( \pmb{W}_t \), the discrepancy may be less than or greater than this bound;
This is contrasted with weak convergence which describes a convergence in probability.
Weak convergence
Let \( \pmb{x}(t) \) be an exact sample path of a generic system of SDEs, and \( \hat{\pmb{x}}(t) \) be some numerical discretization with a maximum step size of \( h \) between time points. The approximation \( \hat{\pmb{x}} \) is said to weakly converge to \( \pmb{x} \) at order \( \gamma \) if there exists a constant \( C \) independent of \( h \) such that, for any \( 2(\gamma+1) \) continuously differentiable function \( g \) of at most polynomial growth \[ \begin{align} \parallel\mathbb{E}\left[ \pmb{g}(\pmb{x}(t)) - \pmb{g}(\hat{\pmb{x}}(t)) \right]\parallel \leq C h^\gamma \end{align} \] for all \( t \in [0, T] \).
In the simple case where \( \pmb{g} \) is the identity, we see a clear distinction between the two modes of convergence:
Therefore, we say that weak convergence measures the accuracy of the empirical statistics generated from an ensemble of approximate path solutions, computing statistics from these realizations.
Strong convergence guarantees weak convergence with at least the same order;
A remarkable fact arises then which differentiates numerical ODEs and SDEs;
The loss of a half order of convergence arises due to the difference between the deterministic Taylor expansion and the Itô-Taylor expansion:
\[ \begin{align} f(W_T) - f(W_0) = \int_{0}^T f'(W_t)\mathrm{d}W_t +\frac{1}{2} \int_{0}^T f''(W_t) \mathrm{d}t. \end{align} \]
When one corrects for the miss-match in the Itô-Taylor expansion including additional terms, one arives at the Milstein scheme.
It is important to note, however, that when the SDE is in terms of additive noise, i.e., \( \mathbf{S}(t, \pmb{x}(t)) := \mathbf{S}(t) \), the correction vanishes and the Euler-Maruyama scheme gains strong convergence order \( 1.0 \).
Nonetheless, the very poor approximation by Euler and Euler-Maruyama leads to the need to derive schemes with better convergence properties.
In the next lecture, we will consider the development of the widely used 4-stage Runge-Kutta scheme and its extension to simulating SDEs.
We will also consider some general practicalities about simulating discrete models, and tangent-linear models, from continuous-time models.