Monte Carlo, Fourier and Finite Element Methods

6D.9 Finite Difference Methods solved PDE on grids. This chapter covers the two other great families of numerical methods, which rest on the two other pillars of Course 6. Monte Carlo methods use the probabilistic representation: by Feynman–Kac (6D.3 Stochastic Differential Equations and Feynman–Kac), a solution of a parabolic equation is an expectation over random paths, and an expectation can be estimated by simulating many paths and averaging. The error decreases like $N^{-1/2}$ with the number of paths, slowly, but independently of the dimension, which makes Monte Carlo the method of choice for options on many assets and for path-dependent contracts. Finite element methods use the weak formulation: by Lax–Milgram (4A.5 Hilbert Spaces and Lax–Milgram, 6A.6 Weak Solutions and Elliptic Regularity), a solution of an elliptic equation is the solution of a variational problem in a Hilbert space, and restricting that problem to a finite-dimensional subspace of piecewise polynomials gives a linear system whose solution is, up to a constant, the best approximation available in the subspace. In between sits Fourier pricing, which uses the explicit characteristic functions of many models and the FFT.

By the end of this chapter you will be able to:

• estimate an expectation by Monte Carlo with a confidence interval, and reduce its variance with antithetic and control variates;
• simulate SDEs by the Euler–Maruyama scheme, prove its strong order $\frac12$, and explain weak order $1$ and multilevel Monte Carlo;
• describe how option prices are computed from characteristic functions by the FFT;
• set up the finite element method for a second-order elliptic problem, and prove Céa's lemma, the interpolation estimate and the Aubin–Nitsche $L^2$ estimate in one dimension;
• compute finite element stiffness matrices, including the cotangent formula on triangles.

10.1Solitaire and the bomb

10.2Monte Carlo

To estimate $m = \mathbb{E}X$ for an integrable random variable $X$ that we can simulate, draw independent copies $X_1, \ldots, X_N$ and form the sample mean $\bar X_N = \frac1N\sum_iX_i$.

The error is $O(N^{-1/2})$ whatever the dimension: to gain one decimal digit, multiply the work by $100$. Every improvement therefore aims at the constant $\sigma$.

Antithetic variates. If $X = f(Z)$ with $Z$ standard normal, $-Z$ has the same law, and $Y = \frac12\big(f(Z) + f(-Z)\big)$ also has mean $m$, with

When $f$ is monotone, $f(Z)$ and $f(-Z)$ are negatively correlated (Exercise 10.8), and averaging the pair costs two evaluations but reduces variance by more than a factor $2$. For the call above the reduction per evaluation was about a factor $2.0$: modest, because the payoff is far from linear.

Control variates. If $W$ is a simulated quantity with known mean $\mathbb{E}W = w$, then $X - b(W - w)$ has mean $m$ for every $b$, and variance $\operatorname{Var}X - 2b\operatorname{Cov}(X, W) + b^2\operatorname{Var}W$, minimised at

where $\rho$ is the correlation of $X$ and $W$. For the call, $W = S_T$ with $\mathbb{E}^{\mathbb{Q}}S_T = S_0e^{rT}$ is a natural control: the estimated correlation gave $\rho^2 \approx 0.85$ and a variance reduction by a factor of about $7$, equivalent to seven times as many samples. The best controls are prices of similar contracts with known formulas, for example a geometric-average Asian option as a control for an arithmetic-average one.

10.3Simulating SDEs: Euler–Maruyama

For $dX = b(X)\,dt + \sigma(X)\,dB$, with step $\Delta t = T/n$ and grid times $t_k = k\Delta t$, the Euler–Maruyama scheme is

the stochastic analogue of Euler's method, with the noise increments sampled exactly. Two kinds of error matter. The strong error $\big(\mathbb{E}|X_T - Y_n|^2\big)^{1/2}$ measures how close individual paths are (needed when the paths themselves matter, or for multilevel methods); the weak error $|\mathbb{E}f(X_T) - \mathbb{E}f(Y_n)|$ measures how close expectations are (all that matters for pricing).

The order $\frac12$ comes from the term $\mathbb{E}|X_s - X_{\eta(s)}|^2 \approx \Delta t$: Brownian motion moves a distance $\sqrt{\Delta t}$ in time $\Delta t$ (6D.1 Brownian Motion). Including the next term of the Itô–Taylor expansion (the Milstein scheme, which adds $\frac12\sigma\sigma'\big((\Delta B)^2 - \Delta t\big)$ in one dimension) raises the strong order to $1$.

The weak order of Euler–Maruyama is $1$: for smooth $f$ and coefficients, $\mathbb{E}f(X_T) - \mathbb{E}f(Y_n) = C_f\Delta t + O(\Delta t^2)$ (statement; Denis Talay and Luciano Tubaro, 1990, who also proved the expansion that justifies Richardson extrapolation). For geometric Brownian motion it can be seen directly: $\mathbb{E}Y_n = x(1 + \mu\Delta t)^n$ while $\mathbb{E}X_T = xe^{\mu T}$, and the difference is about $-\frac12x\mu^2T\Delta te^{\mu T}$, halving when $\Delta t$ halves (Exercise 10.10).

Multilevel Monte Carlo. To reach accuracy $\varepsilon$ with Euler–Maruyama and plain Monte Carlo one needs $\Delta t \sim \varepsilon$ (weak error) and $N \sim \varepsilon^{-2}$ paths (statistical error), a cost of order $\varepsilon^{-3}$. Michael Giles (2008) observed that the expectation on the finest grid can be written as a telescoping sum $\mathbb{E}P_L = \mathbb{E}P_0 + \sum_{\ell=1}^L\mathbb{E}(P_\ell - P_{\ell-1})$ over grids with $\Delta t_\ell = 2^{-\ell}T$, and that the corrections $P_\ell - P_{\ell-1}$, computed on the same Brownian path at two resolutions, have variance of order $\Delta t_\ell$ by the strong error estimate. So few samples are needed on the expensive fine levels, and the total cost falls to order $\varepsilon^{-2}(\log\varepsilon)^2$ (statement): strong convergence, useless on its own for pricing, buys a speedup for weak problems.

10.4Fourier pricing

Many models used in practice are not solvable for prices but have explicit characteristic functions $\phi_T(v) = \mathbb{E}^{\mathbb{Q}}\big[e^{iv\log S_T}\big]$: Lévy (jump) models, and Steven Heston's stochastic volatility model (1993). Prices can then be computed by Fourier inversion. A call price, as a function of the log-strike $k = \log K$, is not integrable (it tends to $S_0$ as $k \to -\infty$), so Peter Carr and Dilip Madan (1999) damp it: for $\alpha > 0$ with $\mathbb{E}S_T^{1+\alpha} < \infty$, the function $e^{\alpha k}C(k)$ is integrable, and its Fourier transform is explicit in terms of $\phi_T$,

(statement; derived by Fubini from $C(k) = e^{-rT}\mathbb{E}(S_T - e^k)^+$, Exercise 10.11). Inverting, $C(k) = \frac{e^{-\alpha k}}{\pi}\int_0^\infty\operatorname{Re}\big(e^{-ivk}\psi(v)\big)\,dv$, and discretising the integral on a grid of $N$ frequencies gives prices on a grid of $N$ strikes by a single FFT (6D.9 Finite Difference Methods), in $O(N\log N)$ operations. This is how the parameters of stochastic volatility models are calibrated to hundreds of quoted options in seconds.

10.5Finite elements

Consider the model problem $-(pu')' + qu = f$ on $(0, 1)$ with $u(0) = u(1) = 0$, $p \geq p_0 > 0$, $q \geq 0$ (the method works the same way in any dimension). Its weak form (6A.6 Weak Solutions and Elliptic Regularity) is: find $u \in V = H^1_0(0, 1)$ with

The form $a$ is bounded, $|a(u, v)| \leq M\|u\|_{H^1}\|v\|_{H^1}$, and coercive, $a(v, v) \geq \alpha\|v\|_{H^1}^2$ (by Poincaré's inequality), so Lax–Milgram gives a unique solution (4A.5 Hilbert Spaces and Lax–Milgram).

The Galerkin method chooses a finite-dimensional subspace $V_h \subset V$ and asks for $u_h \in V_h$ with $a(u_h, v_h) = \ell(v_h)$ for all $v_h \in V_h$. In a basis $\varphi_1, \ldots, \varphi_N$ of $V_h$, with $u_h = \sum_jc_j\varphi_j$, this is the linear system

with the stiffness matrix $A_{ij} = a(\varphi_j, \varphi_i)$, which is positive definite by coercivity. The finite element method chooses $V_h$ to be piecewise polynomials on a mesh, with a local basis, so that $A$ is sparse. The simplest choice in one dimension: nodes $x_j = jh$, $h = \frac{1}{N+1}$, and the hat functions $\varphi_j$, piecewise linear, equal to $1$ at $x_j$ and $0$ at the other nodes. For $p = 1$, $q = 0$, the stiffness matrix is $\frac1h\operatorname{tridiag}(-1, 2, -1)$, the finite difference matrix of 6D.9 Finite Difference Methods times $h$ (Exercise 10.12).

So the Galerkin solution is quasi-optimal: up to the constant $\frac M\alpha$, it is as good as the best approximation of $u$ from $V_h$. (For symmetric $a$ it is exactly the best approximation in the energy norm $a(v, v)^{1/2}$: the Galerkin solution is the orthogonal projection of $u$ onto $V_h$ for the inner product $a$.) Error estimates therefore reduce to approximation theory.

With Céa's lemma: $\|u - u_h\|_{H^1} \leq \frac M\alpha\|u - I_hu\|_{H^1} \leq Ch\|u''\|_{L^2}$, first order in the energy norm. But the $L^2$ error of the Galerkin solution is better than the $H^1$ error by a full power of $h$, by a duality argument of Jean-Pierre Aubin and Joachim Nitsche (1967–68).

The same structure works in higher dimensions on meshes of triangles or tetrahedra: Céa's lemma unchanged; interpolation estimates $\|u - I_hu\|_{H^1} \leq Ch|u|_{H^2}$ for piecewise linear elements on meshes whose triangles are not too thin (statement; Brenner and Scott, chapter 4); and Aubin–Nitsche on convex polygonal domains, where $H^2$ regularity of the dual problem holds. Polynomials of degree $p$ give $O(h^p)$ and $O(h^{p+1})$. For parabolic problems, finite elements in space are combined with the time-stepping schemes of 6D.9 Finite Difference Methods (the method of lines). Unlike finite differences, finite elements handle curved boundaries, unstructured meshes and variable coefficients naturally, and their error theory is the Hilbert space theory of 4A.5 Hilbert Spaces and Lax–Milgram, which is why they dominate engineering simulation.

10.6History

Ulam, von Neumann and Metropolis created the Monte Carlo method at Los Alamos in 1946–49, and Boyle applied it to options in 1977. Gisiro Maruyama analysed the Euler scheme for SDEs in 1955 and Grigori Milstein his higher-order scheme in 1974; Talay and Tubaro's weak error expansion dates from 1990 and Giles's multilevel method from 2008. Richard Courant proposed piecewise linear elements in 1943, engineers developed the method in the 1950s, Clough named it in 1960, and the mathematical theory (Céa's lemma, 1964; Aubin and Nitsche, 1967–68; Ciarlet's book, 1978) followed. Carr and Madan's FFT method dates from 1999 and Heston's model from 1993.

10.7Exercises

Solution

(a) $1.96\times14.7/\sqrt N = 0.05$ gives $N \approx (576)^2 \approx 3.3\times10^5$. (b) The standard deviation becomes $14.7\sqrt{0.15} \approx 5.7$, so $N \approx (1.96\times5.7/0.05)^2 \approx 5.0\times10^4$, about seven times fewer.

Solution

$(f(Z) - f(Z'))(g(Z) - g(Z')) \geq 0$ pointwise, since both factors have the sign of $Z - Z'$. Taking expectations, $2\operatorname{Cov}(f(Z), g(Z)) \geq 0$. With $g(z) = -f(-z)$: $\operatorname{Cov}(f(Z), -f(-Z)) \geq 0$, i.e. $\operatorname{Cov}(f(Z), f(-Z)) \leq 0$.

Solution

$V(b) = \operatorname{Var}X - 2b\operatorname{Cov}(X, W) + b^2\operatorname{Var}W$ is a quadratic in $b$ minimised at $b^* = \operatorname{Cov}(X, W)/\operatorname{Var}W$, with $V(b^*) = \operatorname{Var}X - \frac{\operatorname{Cov}(X, W)^2}{\operatorname{Var}W} = (1 - \rho^2)\operatorname{Var}X$. The factor is $\operatorname{Var}X/V(b^*) = (1 - \rho^2)^{-1}$; for $\rho^2 = 0.85$, about $6.7$.

Solution

$\mathbb{E}[Y_{k+1}\mid Y_k] = Y_k(1 + \mu\Delta t)$ since $\Delta B_k$ has mean zero and is independent of $Y_k$; iterate. The errors are $100(1.005)^{10} - 100e^{0.05} \approx -0.0131$, then $-0.0066$ and $-0.0033$: halving with $\Delta t$. They match $-\frac12x\mu^2T\Delta te^{\mu T} = -0.0131$ for $\Delta t = 0.1$.

Solution

$\int_{-\infty}^se^{(\alpha + iv)k}(e^s - e^k)\,dk = \frac{e^{(\alpha + 1 + iv)s}}{\alpha + iv} - \frac{e^{(\alpha + 1 + iv)s}}{\alpha + 1 + iv} = \frac{e^{(\alpha + 1 + iv)s}}{(\alpha + iv)(\alpha + 1 + iv)}$. Integrating against $e^{-rT}q(s)\,ds$ gives $\frac{e^{-rT}\mathbb{E}e^{i(v - (\alpha + 1)i)s}}{(\alpha + iv)(\alpha + 1 + iv)}$, and $(\alpha + iv)(\alpha + 1 + iv) = \alpha^2 + \alpha - v^2 + i(2\alpha + 1)v$.

Solution

$\varphi_i' = \frac1h$ on $(x_{i-1}, x_i)$ and $-\frac1h$ on $(x_i, x_{i+1})$. So $\int(\varphi_i')^2 = \frac2h$, $\int\varphi_i'\varphi_{i+1}' = -\frac1h$ (overlap on one element), and $0$ otherwise. Dividing the equations $\frac1h(-c_{i-1} + 2c_i - c_{i+1}) = hf(x_i)$ by $h$ gives the difference scheme, with $c_i = u_h(x_i)$.

Solution

In the first case $a(v, v) = \|v\|_{H^1}^2$ and $|a(u, v)| \leq \|u\|\|v\|$, so $M = \alpha = 1$. In the second, $\|v\|_{H^1}^2 \leq (1 + \frac{1}{\pi^2})\int v'^2$, so $a(v, v) \geq \frac{\pi^2}{\pi^2 + 1}\|v\|_{H^1}^2$, $\alpha = \frac{\pi^2}{1 + \pi^2}$ and $M = 1$.

Solution

(a) $\varphi_i$ is constant ($0$) along the edge $jk$, so its gradient is perpendicular to it; it increases from $0$ to $1$ over the distance $h_i$ towards $i$. (b) The gradients are inward normals to the edges $jk$ and $ik$, and the angle between them is $\pi - \theta_k$, so $\nabla\varphi_i\cdot\nabla\varphi_j = -\frac{\cos\theta_k}{h_ih_j}$. Multiplying by the area, $|T|\cdot\frac{-\cos\theta_k}{h_ih_j} = \frac{\frac12|ik||jk|\sin\theta_k(-\cos\theta_k)}{|ik||jk|\sin^2\theta_k} = -\frac12\cot\theta_k$. (c) The edge $ij$ belongs to two triangles, with opposite angles $\alpha_{ij}$ and $\beta_{ij}$, and their contributions add.