Chapter 18

Numerical Applications to Derivative Pricing

18.1 Overview of Deterministic Numerical Methods

18.1.1 Quadrature Formulae

One of the fundamental problems in numerical analysis is the evaluation of integrals. Consider a one-dimensional definite integral of a function f : [a, b] → ℝ,

$I(f)\equiv I(f;[a,b]):={\displaystyle {\int }_a^{b}f(x)dx\text{with}-\infty <a<b<\infty .}$

The function f is said to be integrable (and the integral of f is defined) if the integral I(|f|) exists and is finite. Assuming that a closed-form expression for I(f) is unavailable or intractable, we rely on a numerical evaluation. Any explicit formula that is suitable for providing an approximation of I(f) is called a quadrature formula, or a quadrature rule, or a numerical integration formula. A typical quadrature takes the form

$Q_n(f)\equiv Q_n(f;[a,b]):={\displaystyle \sum _{i=1}^n{w}_{i}f(x_i),}\text{ (18.1)}$

where the nodes x1, x2,..., xn ∈ [a, b] and positive weights w1, w2,..., wn are selected such that Q(f) approximates I(f) for any sufficiently smooth function f. In this section we focus on the approximation of one-dimensional integrals with bounded integration intervals and sufficiently smooth integrands. We consider two main approaches to the construction of a quadrature formula (18.1), namely, the Newton–Cotes quadrature formula and the Gaussian quadrature formula, which are, respectively, covered in the following two subsections.

18.1.1.1 Newton–Cotes Quadrature Formulae

In the Newton–Cotes quadrature formulae, the nodes x1, x2,..., xn are fixed and equally spaced within [a, b]. A closed Newton–Cotes formula includes the endpoints of the integration interval among the nodes, whereas all nodes of an open Newton–Cotes formula are internal points of the integration interval. The weights w1, w2,..., wn, are selected in such a way that the quadrature rule (18.1) is precise for polynomials of the highest possible degree. Suppose that (18.1) is exact for every polynomial of degree no greater than d and is not exact for some polynomial of degree d + 1. Then, the quadrature formula Q is said to have degree of precision d. Since the quadrature rule Q(f) in (18.1) is linear, i.e., Q(f + g) = Q(f)+ Q(g) holds for any functions f and g, we need only test it on monomials xk, k = 0, 1, 2,... to find d. Thus, the degree of precision of a quadrature formula Q is the largest integer d such that I(xk) = Q(xk) for all k = 0, 1,..., d and I(xd+1) ≠ Q(xd+1).

Example 18.1

Construct Newton–Cotes formulae with (a) n = 1 node and (b) n = 2 nodes.

Solution.

(a) There only exists an open Newton–Cotes rule with one node $x_1=\frac{a+b}{2}$ located in the middle of the integration interval [a, b]. We find the weight w1 so that the rule Q1(f) = w1f(x1) is exact for monomials of the highest possible degree. Let us verify whether ${\displaystyle {\int }_a^b{x}^k\text{d}x=w_1{x}_1^k}$ holds for k = 0, 1, 2.

1. If k = 0, then $I(\text{l})={\displaystyle {\int }_a^{b}1}\text{ d}x=b-a=Q_1(1)=w_1$. Thus, w1 = b − a.
2. If k = 1, then $I(x)={\displaystyle {\int }_a^{b}x}\text{ d}x=\frac{1}{2}(b^2-a^2)=Q_1(x)=(b-a)\frac{a+b}{2}$.
3. If k = 2, then $I(x^2)={\displaystyle {\int }_a^b{x}^2}\text{ d}x=\frac{1}{3}(b^3-a^3)\ne Q_1(x^2)=(b-a){\left(\frac{a+b}{2}\right)}^2$.

So the degree of precision of the open Newton–Cotes formula with one node is 1. This formula is called the midpoint rule. It is given by

$I(f)\approx Q_1(f)=(b-a)f\left(\frac{a+b}{2}\right).$

(b) There exist two Newton–Cotes formulae with two nodes, namely, an open formula and a closed formula. Let us consider the case of the closed quadrature formula with the nodes x1 = a and x2 = b. The weights w1 and w2 are as follows:

$\begin{array}{ccccc}\{\begin{array}{c}I(1)=Q_1(1)\\ I(x)=Q_1(x)\end{array}& \Rightarrow & \{\begin{array}{c}{w}_1+w_2=b-a\\ w_{1}a+w_{2}b=\frac{b^2-a^2}{2}\end{array}& \Rightarrow & \{\begin{array}{c}{w}_1=\frac{b-a}{2}\\ w_2=\frac{b-a}{2}\end{array}\end{array}$

The closed Newton–Cotes formula with two nodes denoted Q2(f) is called the trapezoidal rule and is given by

$I(f)\approx Q_2(f)=\frac{b-a}{2}(f(a)+f(b)).$

Let us apply it to f (x) = x2:

$I(x^2)=\frac{b^3-a^3}{3}\ne Q_2(x^2)=\frac{b-a}{2}(a^2+b^2).$

Thus, the degree of precision of the trapezoidal rule is 1.

Alternatively, quadrature formulae can be constructed with the use of the Lagrange interpolation method. Every sufficiently smooth function f can be approximated by the Lagrange polynomial

$p(x){\displaystyle \sum _{i=1}^n{L}_{n,i}(x)f(x_i),\text{ where }{L}_{n,i}(x)={\displaystyle \prod _{\begin{array}{c}j=1\\ j\ne i\end{array}}^n\frac{x-x_j}{x_i-x_j}.}}$

It interpolates f at the nodes x1, x2,..., xn, i.e., p(xi) = f (xi) for all i = 1, 2,..., n. The integral of f can be approximated by the integral of the interpolating polynomial:

$I\left(f\right)\approx I\left(p\right)=I\left({\displaystyle \sum _{i=1}^n{L}_{n,i}\left(x\right)f}\left(x_i\right)\right)={\displaystyle \sum _{i=1}^{n}I\left(L_{n,i}\left(x\right)\right)f\left(x_i\right)={\displaystyle \sum _{i=1}^n{w}_{i}f\left(x_i\right).}}$

As a result, a quadrature formula with weights $w_i:={\displaystyle {\int }_a^b{L}_{n,i}(x)\text{d}x}$ is defined.

For a general integrand f, a quadrature formula is not exact. The error term I(f) − Q(f) can be found by expanding the function f in a Taylor series. The following theorem generalizes our findings related to the degree of precision of a quadrature formula from the above example and provides the formula for the error term.

Theorem 18.1.

Let Qn(f) denote a Newton–Cotes quadrature rule (open or closed) with n nodes.

1. (a) If n is odd and f has n + 1 continuous derivatives, then the degree of precision of Qn is equal to n and the error is I(f) − Qn(f) = C(b − a)n+2 f (n+1)(ξ) for some constant C independent of f and some point ξ ∈ (a, b).
2. (b) If n is even and f has n continuous derivatives, then the degree of precision of Qn is equal to n − 1 and the error is I(f) − Qn(f) = C(b − a)n+1f(n)(ξ ) for some constant C independent of f and some point ξ ∈ (a, b).

The constant of the error term I(f) − Qn(f) can be found using the fact that the quadrature formula is not exact for f (x) = xd+1:

$\begin{array}{c}{I(x^{d+1})-Q_n(x^{d+1})=C{(b-a)}^{d+1}{(x^{d+1})}^{(d+1)}|}_{x=\xi }=C{(b-a)}^{d+1}(d+1)!\\ \Rightarrow C=\frac{1}{{(b-a)}^{d+1}(d+1)!}\left(\frac{b^{d+2}-a^{d+2}}{d+2}-Q_n(x^{d+1})\right).\end{array}$

Let us list some of the Newton–Cotes quadrature rules together with their error terms.

• The midpoint rule (n = 1):

  $I(f)=(b-a)f\left(\frac{a+b}{2}\right)+\frac{{(b-a)}^3}{24}{f}^{″}(\xi ).\text{ (18.2)}$

• The trapezoidal rule (n = 2):

  $I(f)=\frac{{(b-a)}^2}{2}[f(a)+f(b)]-\frac{{(b-a)}^3}{12}{f}^{″}(\xi ).\text{ (18.3)}$

• Simpson’s rule (n = 3):

  $I(f)=\frac{{(b-a)}^2}{6}\left[f(a)+4f\left(\frac{a+b}{2}\right)+f(b)\right]-\frac{{(b-a)}^5}{2880}{f}^{(4)}(\xi ).\text{ (18.4)}$

18.1.1.2 Gaussian Quadrature Formulae

For the Gaussian quadrature formulae, the nodes and weights in (18.1) are selected so as to achieve the highest possible degree of precision. Fixing n, we determine the nodes x1,..., xn and the weights w1,..., wn such that the rule Qn(f) provides the exact value of I(f) for all f (x) ∈ {1, x, x2,..., x2n − 1}. That is, the degree of precision of the quadrature formula is 2n − 1. To determine the 2n unknowns, we form a system of 2n equations (nonlinear equations, in general). For example, if n = 1, then we obtain two simultaneous equations

$w_1=b-a\text{and}(b-a)x_1=\frac{{(b-a)}^2}{2}.$

Solving them, we obtain the midpoint rule.

Example 18.2

Find the Gaussian quadrature rule for n = 2.

Solution. First, we apply the change of variables $x=\frac{b-a}{2}t+\frac{a+b}{2}$ to obtain an integral over [−1, 1]:

${\displaystyle {\int }_a^{b}f(x)dx={\displaystyle {\int }_{-1}^{1}g(t)dt,\text{where}g(t):=\frac{b-a}{2}\cdot }}f\left(\frac{b-a}{2}t+\frac{a+b}{2}\right).$

Second, find the nodes t1, t2 and weights w1, w2 such that ${\displaystyle {\int }_{-1}^{1}g(t)\text{d}t=w_{1}g(t_1)+w_{2}g(t_2)}$ holds for all g(t) ∈ {1, t, t2, t3}:

$\begin{array}{l}g(t)=1\to 2=w_1+w_2,\\ g(t)=t\to 0=w_1{t}_1+w_2{t}_2,\\ g(t)=t^2\to \frac{2}{3}=w_1{t}_1^2+w_2{t}_2^2,\\ g(t)=t^3\to 0=w_1{t}_1^3+w_2{t}_2^3,\end{array}$

Solving the above simultaneous equations gives w1 = w2 = 1, $t_1=-\frac{1}{\sqrt{3}}$, and $t_2=\frac{1}{\sqrt{3}}$. Applying the inverse change of variables, we obtain the following quadrature rule for the integral of f (x) over [a, b]:

${\displaystyle {\int }_a^{b}f(x)dx}\approx \frac{b-a}{2}\left[f\left(\frac{a+b}{2}-\frac{b-a}{\sqrt{3}}\right)+f\left(\frac{a+b}{2}+\frac{b-a}{\sqrt{3}}\right)\right].$

18.1.1.3 Composite Quadrature Formulae

Newton–Cotes integration formulae and Gaussian integration rules can only be applied to integrals with small integration intervals since the error term grows as a power function of the length b − a. Large integration intervals can be first partitioned into smaller subintervals. After this, some simple rule is applied on each subinterval and the final result is given by a sum of approximations for the subintervals:

${\displaystyle {\int }_a^{b}f\left(x\right)\text{d}x}={\displaystyle \sum _{j=1}^m{\displaystyle {\int }_{c_{j-1}}^{c_j}f\left(x\right)\text{d}x\approx {\displaystyle \sum _{j=1}^{m}Q\left(f;\left[c_{j-1},c_j\right]\right),}}}$

where c1, c2,..., cm are points partitioning [a, b] so that a = c0 < c1 <··· < cm = b. The resulting formula is called a composite quadrature rule. For example, let us obtain the composite trapezoidal rule. By splitting the integration interval into m equal subintervals of step size $h=\frac{b-a}{m}$ and applying the trapezoidal rule on each subinterval [cj − 1, cj] with cj = a + jh, j = 0, 1,..., m, we have

$\begin{array}{cc}\hfill {\displaystyle {\int }_a^{b}f(x)dx}=& {\displaystyle \sum _{j=1}^m{\displaystyle {\int }_{c_{j-1}}^{c_j}f(x)dx}}\hfill \\ \hfill =& {\displaystyle \sum _{j=1}^m\frac{c_j-c_{j-1}}{2}[f(c_{j-1})+f(c_j)]-{\displaystyle \sum _{j=1}^m\frac{{(c_j-c_{j-1})}^3}{12}{f}^{″}({\xi }_j)}}\hfill \\ \hfill =& \underset{=T_h(f)(\text{the composite trapezoidal rule})}{\underbrace{\frac{h}{2}\left[f(c_0)+2{\displaystyle \sum _{j=1}^{m-1}f(c_j)+f(c_m)}\right]}}-\underset{\text{the error term}}{\underbrace{\frac{h^3}{12}{\displaystyle \sum _{j=1}^m{f}^{″}({\xi }_j),}}}\hfill \end{array}$

where ξ j ∈ (cj −1, cj) for each j = 1, 2,..., m. By applying the intermediate value theorem, we obtain a more compact formula for the error term: $-\frac{(b-a)h^2}{12}{f}^{″}(\xi )=\mathcal{O}(h^2)$, as h → 0, where ξ ∈ (a, b).

18.1.1.4 Extrapolation and Romberg Integration

Suppose that an integral I(f) is approximated by some quadrature formula Qh(f) which is a function of a small parameter h (e.g., the distance between two neighbouring nodes in a Newton–Cotes integration formula with uniformly spaced nodes). The Richardson extrapolation method is a general procedure that takes two approximations computed with different values of parameter h and combines them in such a way that the resulting approximation is of higher order than the original one. To apply this method, the exponent of h in the leading term of the error of the original approximations must be known. For example, the error term of the composite trapezoidal rule admits the following asymptotic expansion:

$I(f)={\displaystyle {\int }_a^{b}f(x)dx=T_h(f)+C_2{h}^2+C_4{h}^4+C_6{h}^6+\cdots }$

provided that f is a sufficiently smooth function. Here C2,C4,C6,... are some constants independent of h, and Th denotes the composite trapezoidal rule with step size h. Let us take two trapezoidal approximations with steps h and 2h, respectively,

$\begin{array}{l}I(f)=T_h(f)+C_2{h}^2+C_4{h}^4+C_6{h}^6+\cdots \\ I(f)=T_{2h}(f)+4C_2{h}^2+16C_4{h}^4+64C_6{h}^6+\cdots \end{array}$

Multiplying the former by 4 and subtracting the latter from the product, we have

$3I(f)=[4T_h(f)-T_{2h}(f)]+(4C_2-4C_2)h^2+(4C_4-16C_4)h^4+(4C_6-64C_6)h^6+\cdots$

Dividing both sides by 3 gives

$I(f)=\underset{=S_h(f)}{\underbrace{\left[\frac{4T_h(f)-T_{2h}(f)}{3}\right]}}-4C_4{h}^2-20C_6{h}^6+\cdots$

As a result, we have derived a new approximation denoted Sh(f), whose approximation error is of order $\mathcal{O}(h^4)$, instead of $\mathcal{O}(h^2)$, as that of the composite trapezoidal rule. In fact, the combination of two trapezoidal rules with respective step sizes h and 2h gives us Simpson’s rule. Since the exponent of h appearing in the leading term of the approximation error is known, we can apply the extrapolation procedure again to obtain another approximation with an error of order $\mathcal{O}(h^6)$, as h → 0. By applying the Richardson extrapolation method repeatedly on the composite trapezoidal quadrature rule, a sequence of improved approximations of the integral I(f) is generated. This method is known as Romberg integration.

Define the approximation R(k, 0) to be the composite trapezoidal rule with n = 2k + 1 quadrature nodes, i.e., R(k, 0) = Thk(f) with step size $h_k=\frac{b-a}{2^k}$. We begin with the approximation R(L, 0) for some positive integer L. By successively halving the step size, we obtain a series of trapezoidal approximations R(L, 0),R(L + 1, 0),R(L + 2, 0),.... Each successive approximation R(k, 0) is obtained from its predecessor R(k − 1, 0) by adding contributions from the midpoints of every subinterval:

$R(k,0)=\frac{1}{2}R(k-1,0)+h_k{\displaystyle \sum _{i=1}^{2^{k-1}}f(a+(2i-1)h_k)}.$

The approximation error of the trapezoidal approximations is $I(f)-R(k,0)=\mathcal{O}(h_k^2)$. As for the next step, we generate a series of improvements. First, we calculate the extrapolations of the composite trapezoidal approximations:

$R(k,1)=\frac{4^{1}R(k,0)-R(k-1,0)}{4^1-1}=\frac{4R(k,0)-R(k-1,0)}{3},k=L+1,L+2,...$

The rule R(k, 1), is equivalent to the composite Simpson rule with 2k + 1 points. The approximation error of the first improvement is $I(f)-R(k,1)=\mathcal{O}(h_k^4)$. Now generate the second improvement:

$R(k,2)=\frac{4^{2}R(k,1)-R(k-1,1)}{4^2-1}=\frac{16R(k,1)-R(k-1,1)}{15},k=L+2,L+3,...$

The second extrapolation, R(k, 2), is equivalent to Boole’s rule with 2k + 1 points. The approximation error is $I(f)-R(k,2)=\mathcal{O}(h_k^6)$. The third improvement is

$R(k,3)=\frac{4^{3}R(k,2)-R(k-1,2)}{4^3-1}=\frac{48R(k,2)-R(k-1,2)}{47},k=L+3,L+4,...$

The approximation error is $I(f)-R(k,3)=\mathcal{O}(h_k^8)$. In a similar manner, we can generate further improvements provided that the integrand f has continuous high-order derivatives. The jth improvement is given by

$R(k,j)=\frac{4^{j}R(k,j-1)-R(k-1,j-1)}{4^j-1},k=L+j,L+j+1,...$

It is convenient to organize the Romberg method in the form of a lower-triangular table, where the first column is filled with the composite trapezoidal approximations R(L + j, 0), j ≥ 0, and the next column is defined as the extrapolations of the previous one:

$\begin{array}{l}n=2^L+1\hfill \\ n=2^{L+1}+1\hfill \\ n=2^{L+2}+1\hfill \\ n=2^{L+3}+1\hfill \\ ⋮\hfill \end{array}|\begin{array}{l}R(L,0)\hfill \\ R(L+1,0)R(L+1,1)\hfill \\ R(L+2,0)R(L+2,1)R(L+2,2)\hfill \\ R(L+3,0)R(L+3,1)R(L+3,2)R(L+3,3)\hfill \\ ⋮\ddots \hfill \end{array}$

In the above table, n denotes the number of nodes. This table can be filled row by row. It is a common practice to compute new rows until two successive diagonal entries, say R(L + j, j) and R(L + j − 1, j − 1), differ by less than a given error tolerance.

18.1.2 Finite-Difference Methods

18.1.2.1 Finite-Difference Approximations for ODEs

Finite differences are used to numerically solve (initial-)boundary value problems for differential equations by reducing them to systems of linear or nonlinear algebraic equations. As an illustrative example, let us consider the following boundary value problem (BVP) for a second-order ordinary differential equation (ODE):

$\begin{array}{r}\hfill y^{″}(x)=f(x,y(x),y^{\prime }(x)),l<x<r,\text{ (18.5)}\\ \hfill y(l)={\tilde{y}}_l,y(r)={\tilde{y}}_r,\text{ (18.6)}\end{array}$

where f is a continuous function of its arguments, ỹl and ỹr are constant boundary values, and y(x) is an unknown function to be computed on the interval [l, r]. Let us discretize the BVP (18.5)–(18.6) as follows.

1. Introduce a mesh xj = l + j δx, j = 0, 1,..., n + 1, with the step size $\delta x=\frac{r-l}{n+1}$. The exact solution y(·) is approximated by a mesh solution y. at the points xj; that is, y(xj)≈ yj for j = 0, 1,..., n +1.
2. Discretize the ODE (18.5) using finite-difference (FD) approximations of derivatives (see Table 18.1):

   Table 18.1

   Commonly used finite-difference approximations of first and second derivatives.

   Derivative	Finite Difference + Order	Type
   f′(x)	$\frac{f(x+\delta x)-f(x)}{\delta x}+\mathcal{O}(\delta x)$	Forward
   f′(x)	$\frac{f(x)-f(x-\delta x)}{\delta x}+\mathcal{O}(\delta x)$	Backward
   f′(x)	$\frac{f(x+\delta x)-f(x-\delta x)}{2h}+\mathcal{O}(\delta x^2)$	Central
   f″(x)	$\frac{f(x+\delta x)-2f(x)+f(x-\delta x)}{\delta x^2}+\mathcal{O}(\delta x^2)$	Central

   $\begin{array}{c}{y}^{″}(x)=f(x,y(x),y^{\prime }(x))\\ \downarrow \begin{array}{l}\text{(replace all derivatives}\\ \text{by their FD approximations}\end{array}\hfill \\ \frac{y(x_j+\delta x)-2y(x_j)+y(x_j-\delta x)}{\delta x^2}\approx f\left(x_j,y(x_j),\frac{y(x_j+\delta x)-y(x_j-\delta x)}{2\delta x}\right)\\ \downarrow \hfill \\ \frac{y(x_{j+1})-2y(x_j)+y(x_{j-1})}{\delta x^2}\approx f\left(x_j,y(x_j),\frac{y(x_{j+1})-y(x_{j-1})}{2\delta x}\right)\\ \downarrow \begin{array}{l}\text{(obtain a difference equation}\\ \text{for the mesh solution)}\end{array}\hfill \\ \frac{y_{j+1}-2y_j+y_{j-1}}{\delta x^2}=f\left(x_j,y_j,\frac{y_{j+1}-y_{j-1}}{2\delta x}\right).\end{array}$

3. The boundary conditions (18.6) give the values of the mesh solution y0 and yn+1:

   $y_0={\tilde{y}}_l,y_{n+1}={\tilde{y}}_r.$

As a result, we obtain a system of n difference equations with n unknowns y1, y2,..., yn. In general, the resulting difference equations are nonlinear, and the system can be solved by the multivariate Newton method. If the function f in (18.5) is a linear function of y and y′ (i.e., in the case of a linear second order ODE), the resulting system is linear and it can be solved by the Gaussian elimination approach or by an iterative method such the Gauss–Seidel and Jacobi methods. For example, consider the following linear ODE: y″(x) = a(x)+ b(x)y(x)+ c(x)y′(x). The respective system of difference equations for the mesh solution has a triangular matrix of coefficients and is given by

$\{\begin{array}{l}-(2+b_1\delta x^2)y_1+(1-c_1\delta x/2)=\delta x^2{a}_1-(1+c_1\delta x/2){\tilde{y}}_l,\hfill \\ (1+c_i\delta x/2)y_{i-1}-(2+b_i\delta x^2)y_i+(1-c_i\delta x/2)y_{i+1}=\delta x^2{a}_i,i=2,3,...,n-1,\hfill \\ (1+c_n\delta x/2)y_{n-1}-(2+b_n\delta x^2)y_n=\delta x^2{a}_n-(1-c_n\delta x/2){\tilde{y}}_r,\hfill \end{array}$

where ai := a(xi), bi := b(xi), ci := c(xi) for i = 1, 2,..., n. The system of linear equations (for the mesh solution y·) can be solved by the Gaussian elimination method. Since the matrix of coefficients is tridiagonal, the application of the Gaussian elimination method requires $\mathcal{O}(n)$ arithmetic operations.

18.1.2.2 Second-Order Linear PDEs

The finite-difference approach can be applied to solve numerically a partial differential equation (PDE) that involves more than one independent variable. By replacing all partial derivatives with their respective finite-difference approximations, we can reduce a differential equation to a system of difference equations. Consider a PDE of the form

$A\frac{{\partial }^{2}u}{\partial x^2}+2B\frac{{\partial }^{2}u}{\partial x\partial y}+C\frac{{\partial }^{2}u}{\partial y^2}+D\frac{\partial u}{\partial x}+E\frac{\partial u}{\partial y}+Fu=G(x,y),\text{ (18.7)}$

where A,B,C,...,F are constants, G(x, y) is a known function, and u = u(x, y) is the solution to be determined. The PDE (18.7) is called

• an elliptic-type PDE, if AC − B2 > 0 (e.g., the Laplace equation $\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}=0$);
• a parabolic-type PDE, if AC − B2 = 0 (e.g., the heat equation $\frac{\partial u}{\partial t}=c^2\frac{{\partial }^{2}u}{\partial x^2}$);
• a hyperbolic-type PDE, if AC − B2 < 0 (e.g., the wave equation $\frac{{\partial }^{2}u}{\partial t^2}=c^2\frac{{\partial }^{2}u}{\partial x^2}$.

Let B = 0 and C = 0, then Equation (18.7) is of a parabolic type. Relabel y→ t to obtain the following equation:

$\frac{\partial u}{\partial t}=\tilde{A}\frac{{\partial }^{2}u}{\partial x^2}+\tilde{D}\frac{\partial u}{\partial x}+\tilde{F}u+\tilde{G}(x,y).$

The simplest parabolic-type PDE is the heat equation

$\frac{\partial u}{\partial t}=c^2\frac{{\partial }^{2}u}{\partial x^2}.\text{ (18.8)}$

This equation models the temperature in a thin insulated rod. Here x and t are called spacial and temporal variables, respectively. As is well known, the Black–Scholes equation for the GBM model can be reduced to the heat equation by applying a change of variables and a scale transformation. The heat equation can be solved numerically or analytically using finite-differences or spectral expansions, respectively.

A differential equation may have multiple solutions, and hence a differential problem becomes well-formulated only when a PDE is combined with initial and boundary conditions. We consider a typical initial boundary value (IBV) problem for the heat equation:

$\{\begin{array}{ll}\frac{\partial u(t,x)}{\partial t}=c^2\frac{{\partial }^{2}u(t,x)}{\partial x^2}\hfill & \text{for }0\le x\le \ell ,0<t\le T,\hfill \\ u(0,x)=f(x)\hfill & \text{for }0\le x\le \ell ,\hfill \\ u(t,0)=g_0(t)\hfill & \text{for }0\le t\le T,\hfill \\ u(t,\ell )=g_1(t)\hfill & \text{for }0\le t\le T,\hfill \end{array}\text{ (18.9)}$

where f, g0, and g1 are continuous functions, and c is a constant. The unknown solution is a function u(t, x) defined on the rectangle D := [0,T]× [0, ℓ]. In the following section, we describe three finite-difference schemes for solving (18.9) numerically. The computational solution is a function defined on a two-dimensional grid (or mesh) of points in the domain D. Being given a mesh solution, the approximation solution can be calculated at every point of the domain D using interpolation methods.

18.1.2.3 Finite-Difference Approximations for the Heat Equation

Our goal is to approximate the solution u to the IBV problem (18.9) on a spacial-temporal mesh of points. First, we define a spacial mesh at which the solution is to be approximated. Divide [0, ℓ] into n + 1 equal subintervals of length $\delta x=\frac{\ell }{n+1}$. Introduce the spacial mesh points xj = j δx, j = 0, 1,..., n + 1. Second, the derivatives with respect to the spacial variable x are approximated by their central finite-difference approximations (with the order of approximation $\mathcal{O}(\delta x^2))$. Inserting the finite-difference formula for uxx in the heat equation (18.8) gives

$\frac{\partial u(t,x)}{\partial t}=c^2\frac{u(t,x-\delta x)-2u(t,x)+u(t,x+\delta x)}{\delta x^2}+\mathcal{O}(\delta x^2).\text{ (18.10)}$

Let uj(t) denote the semi-discrete approximation to u(t, xj) for all j = 0, 1,..., n + 1. Each uj is a function of time, not a function of the spacial variable x. Discarding the error term and writing (18.10) at every internal spacial node xj, we have a system of ordinary differential equations

$\frac{\text{d}{u}_j(t)}{\text{d}t}=c^2\frac{u_{j-1}(t)-2u_j(t)+u_{j+1}(t)}{\delta x^2},j=1,2,...,n,\text{ (18.11)}$

supplemented by the boundary conditions u0(t) = g0(t) and un+1(t) = g1(t). Employing the initial condition on u(t, x) gives

$u_j(0)=f(x_j),j=0,1,...,n+1.$

As a result, we obtain a Cauchy problem for a system of ODEs which can be written in matrix-vector form as follows:

$\begin{array}{ll}\frac{\text{du}(t)}{\text{d}t}=-\frac{c^2}{\delta x^2}\text{Au}(t)+\frac{c^2}{\delta x^2}\text{g}(t)\text{ subject to u}(0)=\text{f, where}\hfill & \hfill \text{(18.12)}\\ \text{u}(t):=\left[\begin{array}{c}{u}_1(t)\\ u_2(t)\\ ⋮\\ u_n(t)\end{array}\right]\in {ℝ}^n,\text{ f}:=\left[\begin{array}{c}f(x_1)\\ f(x_2)\\ ⋮\\ f(x_n)\end{array}\right]\in {ℝ}^n,\text{ g}(t):=\left[\begin{array}{c}{g}_0(t)\\ 0\\ ⋮\\ 0\\ g_1(t)\end{array}\right]\in {ℝ}^n,\hfill & \hfill \text{(18.13)}\\ \text{A}=\left[\begin{array}{cccccc}2& -1& 0& \cdots & 0& 0\\ -1& 2& -1& \cdots & 0& 0\\ 0& -1& 2& \cdots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& \cdots & 2& -1\\ 0& 0& 0& \cdots & -1& 2\end{array}\right]\in {ℝ}^{n\times n}.\hfill & \hfill \text{(18.14)}\end{array}$

The coefficient matrix A is called the discretization matrix. The key feature of A is that it is a sparse matrix. All its nonzero elements are located on its main diagonal and upper and lower adjacent diagonals. Such a matrix is called a tridiagonal matrix or a band matrix with bandwidth 3.

For the numerical solution of the system of ODEs (18.12), we can use a time-discretization scheme by introducing a temporal mesh ti = i δt, i = 0, 1,..., m with a mesh size $\delta t:=\frac{T}{m}$. Using different time discretizations, we derive three computational schemes, namely, the explicit Euler method, the implicit Euler method, and the Crank–Nicolson method. Applying one of these three schemes, we obtain a fully discretized approximation solution

$u_{i,j}\approx u_j(t_i)\approx u(t_j,x_j),i=0,1,...,m,j=0,1,...,n+1,$

defined on the two-dimensional temporal-spatial mesh (see Figure 18.2). This is called a mesh solution to the IBV problem (18.9).

The Explicit Method

The explicit Euler method is obtained by replacing the time derivative by a forward finitedifference approximation. For a continuously differentiable function f, we have

$\frac{\text{d}f(t)}{\text{d}t}=\frac{f(t+\delta t)-f(t)}{\delta t}+\mathcal{O}(\delta t).$

Applying this formula for $\frac{\partial u}{\partial t}$ in (18.10) gives

$\frac{u(t_{i+1},x_j)-u(t_i,x_j)}{\delta t}=c^2\frac{u(t_i,x_{j-1})-2u(t_i,x_j)+u(t_i,x_{j+1})}{\delta x^2}+\mathcal{O}(\delta t+\delta x^2)\text{ (18.15)}$

for each i = 0, 1,..., m − 1 and j = 1, 2,..., n. Discarding the error term $\mathcal{O}(\delta t+\delta x^2)$ and replacing every u(ti, xj) by its approximation ui, j, we obtain the following finite-difference approximation to the heat equation (18.8):

$u_{i+1,j}=\alpha u_{i,j-1}+(1-2\alpha )u_{i,j}+\alpha u_{i,j+1},i=1,2,...,n,j=0,1,...,m-1,\text{ (18.16)}$

where $\alpha :=c^2\frac{\delta t}{\delta x^2}$. The explicit method starts with the values u0, j (for i = 0), which are given by the initial time condition:

$u_{0,j}=f(x_j),j=0,1,...,n+1.\text{ (18.17)}$

Additionally, the values ui,0 and ui,n+1 for all i = 1, 2,..., m are known thanks to the boundary conditions:

$u_{i,0}=g_0(t_i),u_{i,n+1}=g_1(t_i),i=1,2,...,m.\text{ (18.18)}$

The formula (18.16) is an explicit expression for every approximation ui+1, j at time ti+1 in terms of the solution at the previous time ti. Therefore, this method is called the explicit Euler method. For a fixed i, the values ui+1, j for all j = 0, 1,..., n + 1 are computed from the values {ui, j}j = 0,1,..., n+1 of the ith time level. So the explicit rule (18.16) allows us to advance from the ith time level to the (i + 1)st time level. The accuracy of the explicit method is $\mathcal{O}(\delta t+\delta x^2)$.

The explicit method (18.16) can also be written in matrix-vector form. Let us denote u(i) := [ui,1, ui,2, ..., ui,n]T. For each i = 0, 1, 2,..., m, the vector u(i) is an approximation to u(ti). Inserting the forward finite-difference approximation of the time derivative $\frac{\text{du}(t)}{\text{d}t}$,

${\frac{\text{d}\text{u}(t)}{\text{d}t}|}_{t=t_i}=\frac{\text{u}(t_{i+1})-\text{u}(t_i)}{\delta t}+\mathcal{O}(\delta t)\approx \frac{{\text{u}}^{(i+1)}-{\text{u}}^{(i)}}{\delta t},$

into Equation (18.12) written at t = ti gives

$\frac{{\text{u}}^{(i+1)}-{\text{u}}^{(i)}}{\delta t}=-\frac{c^2}{\delta x^2}{\text{Au}}^{(i)}+\frac{c^2}{\delta x^2}\text{g}(t_i),$

which leads to the matrix-vector equations

${\text{u}}^{(i+1)}=({\text{I}}_n-\alpha \text{A}){\text{u}}^{(i)}+\alpha {\text{g}}^{(i)},i\ge 0,\text{ (18.19)}$

where In denotes the n-by-n identity matrix and g(i) := g(ti). The formula (18.19) is an iterative rule which starts at

${\text{u}}^{\text{(0)}}=\text{f}.\text{ (18.20)}$

The problem (18.19)-(18.20) admits a closed-form solution

${\text{u}}^{\text{(k)}}={({\text{I}}_n-\alpha \text{A})}^k\text{f}+\alpha {\displaystyle \sum _{i=1}^{k-1}{({\text{I}}_n-\alpha \text{A})}^{k-i-1}{g}^{(i)},}0\le k\le m.$

The Implicit Method

To derive the implicit method, we replace the time derivative of the solution by a backward finite-difference approximation:

$\frac{\text{d}f(t)}{\text{d}t}=\frac{f(t)-f(t-\delta t)}{\delta t}\mathcal{O}(\delta t).$

Applying this formula for $\frac{\partial f(t)}{\partial t}$ in (18.10) gives

$\frac{u(t_{i+1},x_j)-u(t_i,x_j)}{\delta t}=c^2\frac{u(t_{i+1},x_{j-1})-2u(t_{i+1},x_j)+u(t_{i+1},x_{j+1})}{\delta x^2}+\mathcal{O}(\delta t+\delta x^2)\text{ (18.21)}$

for each i = 0, 1,..., m − 1 and j = 1, 2,..., n. Discarding the error term $\mathcal{O}(\delta t+\delta x^2)$ and replacing every u(ti, xj) by its approximation ui, j, we have

$-\alpha u_{i+1,j-1}+(2\alpha +1)u_{i+1,j}-\alpha u_{i+1,j+1}=u_{i,j},i=1,2,...,n,j=0,1,...,m-1,\text{ (18.22)}$

subject to the initial condition (18.17) and boundary condition (18.18). As is seen from (18.22), the approximations ui+1,j of the (i + 1)st time level are coupled and cannot be calculated explicitly from approximations ui,j of the ith time level. Hence, this scheme is called the implicit method.

The implicit method (18.22) can be written in matrix-vector form. Inserting the backward finite-difference approximation of the time derivative $\frac{\text{d}\text{u}(t)}{\text{d}t}$,

${\frac{\text{d}\text{u}(t)}{\text{d}t}|}_{t=t_{i+1}}=\frac{\text{u}(t_{i+1})-\text{u}(t_i)}{\delta t}+\mathcal{O}(\delta t)\approx \frac{{\text{u}}^{(i+1)}-{\text{u}}^{(i)}}{\delta t},$

into Equation (18.12) written at t = ti+1 gives

$\frac{{\text{u}}^{(i+1)}-{\text{u}}^{(i)}}{\delta t}=-\frac{c^2}{\delta x^2}{\text{Au}}^{(i+1)}+\frac{c^2}{\delta x^2}{\text{g}}^{(i+1)}.$

As a result, we obtain the following computational formula:

$({\text{I}}_n+\alpha \text{A}){\text{u}}^{(i+1)}={\text{u}}^{(i)}+\alpha {\text{g}}^{(i+1)},i\ge 0,\text{ (18.23)}$

subject to the initial condition (18.20). The order of approximation is $\mathcal{O}(\delta t+\delta x^2)$ which is the same as that of the explicit method.

To advance to the (i + 1)st time level and find the mesh solution u(i+1), we need to solve a system of linear equations with a tridiagonal coefficient matrix C = In + αA. In general, the numerical solution of a system with n equations and a band coefficient matrix requires $\mathcal{O}(n)$ operations. So the linear system Cx = b for x = u(i+1) with the tridiagonal matrix C and vector b = u(i) + αg(i+1) can be solved by the Gaussian elimination method with n − 1 eliminations (during the forward step) and n − 1 substitutions (during the backward step).

The Crank–Nicolson Method

The Crank–Nicolson method is one of the most well-known implicit finite-difference methods thanks to its stability properties and high order of approximation. To derive the Crank– Nicolson method, we first approximate the time derivative of u(t) at point t = ti + δt/2 with the central finite-difference formula that has the second order of approximation:

$\frac{\text{u}(t_{i+1})-\text{u}(t_i)}{\delta t}+\mathcal{O}(\delta t^2)=\frac{\text{d}\text{u}(t)}{\text{d}t}{|}_{t=t_i+\delta t/2}=-\frac{c^2}{\delta x^2}\text{Au}(t_i+\delta t/2)+\frac{c^2}{\delta x^2}\text{g}(t_i+\delta t/2).$

Second, we use the average of the spatial approximations at times ti+1 and ti in the right-hand side of the above equation:

$\frac{\text{u}(t_{i+1})-\text{u}(t_i)}{\delta t}=-\frac{c^2}{2\delta x^2}\text{A}(\text{u}(t_i)+\text{u}(t_{i+1}))+\frac{c^2}{2\delta x^2}(\text{g}(t_i)+\text{g}(t_{i+1}))+\mathcal{O}(\delta t^2+\delta x^2).$

As a result, we obtain the following computational scheme, called the Crank–Nicolson method:

$\left({\text{I}}_n+\frac{\alpha }{2}\text{A}\right){\text{u}}^{(i+1)}=\left({\text{I}}_n-\frac{\alpha }{2}\text{A}\right){\text{u}}^{(i)}+\frac{\alpha }{2}\left({\text{g}}^{(i)}+{\text{g}}^{\text{(i+1)}}\right),i\ge 0.\text{ (18.24)}$

The method starts with the initial condition (18.20). As is seen from (18.24), the Crank– Nicolson method is an implicit computational scheme since we need to solve a system of linear equations to advance the mesh solution. However, the order of approximation of the Crank–Nicolson method is $\mathcal{O}(\delta t^2+\delta x^2)$, which is higher than the order of approximation of both explicit and implicit Euler methods. Additionally, the Crank–Nicolson method is unconditionally stable, meaning that it is not sensitive to roundoff errors regardless of the values of the mesh sizes δt and δx.

18.1.2.4 Stability Analysis

The stability of a computational method is closely associated with numerical errors. A finite-difference scheme is stable if the roundoff error made at one iteration does not amplify as the computations are continued. That is, if the errors decay or remain bounded, the numerical method is stable. If, on the contrary, the errors grow with time, the numerical method is said to be unstable. To investigate the stability of the explicit and implicit Euler methods and the Crank-Nicolson scheme, we represent them in the form of an iteration rule:

${\text{u}}^{(i+1)}={\text{Qu}}^{(i)}+{\text{c}}^{(i)},i\ge 0,\text{ (18.25)}$

where

$\begin{array}{c}\text{Q}=\{\begin{array}{ll}{\text{I}}_n-\alpha \text{A}\hfill & (\text{the explicit method})\hfill \\ {({\text{I}}_n+\alpha \text{A})}^{-1}\hfill & (\text{the implicit method})\hfill \\ {\left({\text{I}}_n+\frac{\alpha }{2}\text{A}\right)}^{-1}\left({\text{I}}_n-\frac{\alpha }{2}\text{A}\right)\hfill & (\text{the Crank-Nicolson method})\hfill \end{array}\\ {\text{c}}^{(i)}=\{\begin{array}{ll}\alpha {\text{g}}^{(i)}\hfill & (\text{the explicit method})\hfill \\ \alpha {({\text{I}}_n+\alpha \text{A})}^{-1}{\text{g}}^{(i+1)}\hfill & (\text{the implicit method})\hfill \\ \frac{\alpha }{2}{({\text{I}}_n+\alpha \text{A})}^{-1}\left({\text{g}}^{(i)}+{\text{g}}^{(i+1)}\right)\hfill & (\text{the Crank-Nicolson method})\hfill \end{array}\end{array}$

Let e(0) be the error introduced to the method initially. That is, in actual computations we start with the vector ũ(0) = u(0) + e(0), which includes the introduced error. Of interest is the propagation of this error during computations as we advance the solution in time. Let ũ(i) denote the mesh solution (of the ith time step) that includes the error. This perturbed mesh solution satisfies the rule (18.25):

${\tilde{u}}^{\left(i+1\right)}=Q{\tilde{u}}^{\left(i\right)}+{\text{c}}^{\left(i\right)}.$

Subtracting (18.25) from this equation gives us the following recurrence formula for the accumulated error e(i) := ũ(i) − u(i):

${\text{e}}^{(i+1)}={\text{Q}\tilde{\text{u}}}^{(i)}-{\text{Qu}}^{(i)}={\text{Qe}}^{(i)}.$

Hence, the error at the kth time step is

${\text{e}}^{(k)}={\text{Q}}^k{\text{e}}^{(0)},k\ge 0.\text{ (18.26)}$

Since the approximation error is a vector, we need matrix and vector norms to measure the magnitude of the error. Examples of matrix norms include

• ||Q||∞ := maxi∑j | Qij | = ||Q┬||1 (the biggest row-sum of magnitudes);
• ||Q||1 := maxj∑j | Qij | = ||Q┬||∞ (the biggest column-sum of magnitudes);
• ${||\text{Q}||}_2:=\sqrt{{\lambda }_{\max }({\text{Q}}^{{\rm T}}\text{Q})}$ (the spectral norm given by the square root of the largest eigenvalue of the positive-semidefinite matrix Q┬Q).

Here Q is an n-by-n matrix of real elements Qij. Recall that a matrix norm || · ||* on ℝn×n is said to be compatible with a vector norm || · || on ℝn if

$||\text{Qx}||\le {||\text{Q}||}_{*}\cdot ||\text{x}||$

holds for all Q ∈ ℝn×n and all x ∈ ℝn.

According to the definition given in the beginning of this subsection, a finite-difference method with the coefficient matrix Q is stable if the errors decay during the computations. For the error e(k) to decay and eventually vanish as k → ∞, that is, limk→∞ ||e(k) || = 0, it is necessary and sufficient that

${||{\text{Q}}^k||}_{*}<1\text{ (18.27)}$

for some matrix norm ||· ||* and some positive integer K. Indeed, using properties of matrix and vector norms and the representation (18.26), we obtain

$0\le ||{\text{e}}^{(k)}||\le {||\text{Q}||}_{*}^b\cdot {||{\text{Q}}^k||}_{*}^a\cdot ||{\text{e}}^{(0)}||,\text{ (18.28)}$

where the nonnegative integers a and b < K are defined so that k = a· K + b. Provided that (18.27) holds, the right-hand side of (18.28) converges to 0 as k→ ∞. Note that if QK ||* = 1, then the errors remain bounded as k→ ∞, i.e., lim supk→∞ ||e(k) || <∞.

The stability condition can be formulated with the use of the spectral radius. The spectral radius ρ(Q) of an n-by-n matrix Q is defined by the maximum magnitude of the eigenvalues λ1,λ2,...,λn of Q. That is,

$\rho (\text{Q}):=\text{max}\{|{\lambda }_1|,|{\lambda }_2|,...,|{\lambda }_n|\}.$

Note that for a symmetric matrix, the 2-norm is equal to the spectral radius of the matrix, i.e.,||Q||2 = ρ(Q) if Q = QT.

Proposition 18.2.

A finite-difference method is stable if the spectral radius of the matrix Q be less than one, i.e.,

$\rho (\text{Q})<1.\text{ (18.29)}$

Proof. Let us prove that the conditions (18.27) and (18.29) are equivalent. First, use the fact that the spectral radius does not exceed ${||{\text{Q}}^k||}_{*}^{1/k}$ for any matrix norm || · ||* and any positive integer k. That is, $\rho (\text{Q})\le {||{\text{Q}}^k||}_{*}^{1/k}$. Therefore, ||Qk||* < 1 implies that ρ(Q) < 1. On the other hand, according to Gelfand’s formula, for any matrix norm || · ||* we have ${\lim }_{k\to \infty }{||{\text{Q}}^k||}_{*}^{1/k}=\rho (\text{Q})$. Let ρ(Q) = α0 < 1. Then, there exists k so that we have

$|{||{\text{Q}}^k||}_{*}^{1/k}-\rho (\text{Q})|<\frac{1-{\alpha }_0}{2}\Rightarrow {||{\text{Q}}^k||}_{*}^{1/k}<\rho (\text{Q})+\frac{1-{\alpha }_0}{2}=\frac{1+{\alpha }_0}{2}<1.$

Stability Analysis of the Explicit Method

Let us calculate the spectral radius of Q = In − αA, where $\alpha =c^2\frac{\delta t}{\delta x^2}>0$ and A is given by (18.14). Let λ be an eigenvalue of A and x be a respective eigenvector, i.e., Ax = λx holds.1 We have

$\text{Qx}=({\text{I}}_n-\alpha \text{A})\text{x}={\text{I}}_n\text{x}-\alpha \text{Ax}=\text{x}-\alpha \lambda \text{x}=(1-\alpha \lambda )\text{x}\text{.}$

That is, all eigenvalues of the matrix Q are of the form μ = 1 − αλ. Therefore, to calculate ρ(Q) we only need to know the smallest and largest eigenvalues of the matrix A.

Proposition 18.3.

The eigenvalues of the n-by-n matrix A defined by (18.14) are

${\lambda }_j=2-2\cos {\theta }_j,$

which are associated with the respective n-by-1 eigenvectors

${\text{x}}_j={[\sin ({\theta }_j),\sin (2{\theta }_j),...,\sin (n{\theta }_j)]}^{{\rm T}},\text{ (18.30)}$

where ${\theta }_j=\frac{j\pi }{n+1}$ for all j = 1, 2,..., n.

Proof. It is left as an exercise for the reader to verify this.

Due to Proposition 18.3, the eigenvalues of Q are

${\mu }_j=1-\alpha {\lambda }_j=1-\alpha (2-2\cos {\theta }_j)=1-4\alpha {\sin }^2\left(\frac{j\pi }{2(n+1)}\right),j=1,2,...,n.$

Hence, the stability requirement, ρ(Q) = maxj|μj| < 1, can be written as follows:

$-1<1-4\alpha {\sin }^2\left(\frac{j\pi }{2(n+1)}\right)<1,j=1,2,...,n.$

Since α > 0, the right inequality is automatically satisfied. Rearranging the left inequality gives

$\frac{1}{2}>\alpha {\sin }^2\left(\frac{j\pi }{2(n+1)}\right),j=1,2,...,n.$

The sequence $\left\{{\sin }^2\left(\frac{j\pi }{2(n+1)}\right);j=1,2,...,n\right\}$ attains its maximum value at j = n. The maximum approaches 1 from below as n→ ∞. Thus, the above restrictions on α are satisfied for any n ≥ 1 iff $\alpha \le \frac{1}{2}$. To insure the stability of the explicit method, we require that

$c^2\frac{\delta t}{\delta x^2}\le \frac{1}{2}\text{ (18.31)}$

holds. The explicit method is said to be conditionally stable. The condition (18.31) can also be interpreted as follows: every time we double the number of mesh points along the x-axis (and hence halve the mesh size δx), we need to quadruple the number of mesh points along the t-axis.

The stability of an explicit finite-difference method can also be investigated using the matrix norm. The explicit Euler method is stable if ||Q||* ≤ 1 for some norm || · ||*. Calculating the biggest row-sum of magnitudes of Qij we have

${||\text{Q}||}_{\infty }=|1-2\alpha |+2|\alpha |.$

If $\alpha >\frac{1}{2}$, then ||Q||∞ = 4α − 1 > 1 and the scheme is unstable. If $\alpha \le \frac{1}{2}$, then ||Q||∞ = 1 and the scheme is stable.

Stability Analysis of the Implicit Method

The eigenvalues νj of the matrix In + αA can be calculated from those of A as follows:

${\nu }_j=1+\alpha {\lambda }_j=1+4\alpha {\sin }^2\left(\frac{j\pi }{2(n+1)}\right),j=1,2,...,n.$

All of them are strictly greater than one. In particular this fact ensures that the inverse Q = (In + αA)−1 exists. Since Q is a nonsingular matrix, its eigenvalues μj are reciprocals of the eigenvalues νj of In + αA. Therefore, the implicit method is stable since

$0<\frac{1}{1+4\alpha {\sin }^2\left(\frac{j\pi }{2(n+1)}\right)}<1$

for all j = 1, 2,..., n and for all α > 0. The implicit method is said to be unconditionally stable since there is no restriction on the choice of the mesh sizes δx and δt.

Stability Analysis of the Crank–Nicolson Method

An important observation is that the matrices $\text{B}\equiv {\text{I}}_n+\frac{\alpha }{2}\text{A}$ and $\text{C}\equiv {\text{I}}_n-\frac{\alpha }{2}\text{A}$ have the same eigenvectors as those of the matrix A. The eigenvalues ${\nu }_j^B$ of B and ${\nu }_j^C$ of C associated with the eigenvector xj defined by (18.30) are, respectively, given by

${\nu }_j^B=1+\frac{\alpha }{2}{\lambda }_j\text{ and }{\nu }_j^C=1-\frac{\alpha }{2}{\lambda }_j.$

Moreover, since α > 0, all eigenvalues of B are strictly positive; in fact, they all lie in the interval (1, 1+ α). Hence the matrix B is nonsingular. The eigenvalues of the inverse B− 1 are of the form $\frac{1}{{\nu }_j^B}$ since

${\text{Bx}}_j={\nu }_j^B\text{x}\iff \frac{1}{{\nu }_j^B}{\text{x}}_j={\text{B}}^{-1}{\text{x}}_{j.}$

Therefore, the eigenvectors of Q = B−1C are the same as those of A. The eigenvalues μj of Q are

${\mu }_j=\frac{1-\frac{\alpha }{2}{\lambda }_j}{1+\frac{\alpha }{2}{\lambda }_j}=\frac{1-2\alpha {\sin }^2\left(\frac{j\pi }{2(n+1)}\right)}{1+2\alpha {\sin }^2\left(\frac{j\pi }{2(n+1)}\right)}$

for j = 1, 2,..., n. Indeed, we have

${\text{Qx}}_j={\text{B}}^{-1}({\text{Cx}}_j)={\nu }_j^C{\text{B}}^{-1}{\text{x}}_j=\frac{{\nu }_j^C}{v_j^B}{\text{x}}_j.$

The magnitude of every eigenvalue μj is strictly less than one. Therefore, the Crank– Nicolson method is unconditionally stable.

18.2 Pricing European Options

Valuation of a European derivative is probably the most elementary problem of computational finance. In the case of a single asset price model, this problem reduces to the numerical evaluation of a one-dimensional integral. Hence, the result can be achieved with the use of quadrature rules. On the other hand, according to the risk-neutral pricing formula, the no-arbitrage derivative value is the expected value of the discounted payoff function. The mathematical expectation can be evaluated by the Monte Carlo method, which is especially efficient for pricing basket (i.e., multi-asset) options and path-dependent derivatives, the numerical evaluation of which is discussed in later sections. The PDE formulation of the derivative price problem allows for the application of finite-difference schemes. Finally, we may use multinomial tree methods (also known as lattice methods) that can be derived as a result of a discrete-time approximation of a continuous-time asset price model. In fact, tree methods are closely related to Monte Carlo methods. Furthermore, they can also be considered as a special case of explicit finite-difference methods.

18.2.1 Pricing European Options by Quadrature Rules

Let us consider a standard European-style option maturing at time T with payoff function Λ(S(T)), where S(T) is the time-T price of the underlying asset. Here, we assume that the underlying price process {S(t)}t≥0 is time homogeneous and Markovian and that the (risk-neutral) probability law of the price S(T) conditional on S(t) = S for 0 ≤ t < T is known. We can therefore express the option value v(τ, S) as function of the time to maturity τ = T − t and spot S = S(t). As we have seen in previous chapters, the no-arbitrage option value is given by the risk-neutral pricing formula and can be evaluated as a one-dimensional integral (using the risk-free bank account B(t) = ert as numéraire):

$\upsilon (\tau ,S)={\text{e}}^{-rT}{\tilde{\text{E}}}_{t,S}[\Lambda (S(T))]={\text{e}}^{-rT}{\displaystyle {\int }_0^{\infty }\Lambda (S^{\prime })\tilde{p}(\tau ;S,S^{\prime })\text{d}{S}^{\prime }}\text{ (18.32)}$

Here, $\tilde{p}$(τ ; S, S′) denotes the risk-neutral transition PDF of the asset price process. Some examples of processes for which the PDF $\tilde{p}$ is known in closed form include the log-normal (GBM) model and the constant elasticity of variance (CEV) diffusion model. There are also other interesting families of more generally state-dependent nonlinear volatility diffusion models that admit analytically closed-form expressions for the transition PDF. Some jump-diffusion processes such as the variance gamma (VG) model and the normal inverse Gaussian (NIG) model, which are constructed from Brownian motion with the use of a stochastic time change, also admit transition PDFs in closed form.

We can also compute the Greeks of the contract such as the delta $\Delta (\tau ,S)=\frac{\partial \upsilon }{\partial S}(\tau ,S)$. Provided that we can change the order of integration and differentiation and that the integrals obtained exist, we have

$\Delta (\tau ,S)={\text{e}}^{-rT}\frac{\partial }{\partial S}{\displaystyle {\int }_0^{\infty }\Lambda (S^{\prime })\tilde{p}(\tau ;S,S^{\prime })\text{d}{S}^{\prime }}={\text{e}}^{-rT}{\displaystyle {\int }_0^{\infty }\Lambda (S^{\prime })\frac{\partial \tilde{p}}{\partial S}(\tau ;S,S^{\prime })\text{d}{S}^{\prime }}.\text{ (18.33)}$

In some cases, the integrals in (18.32)–(18.33) can be computed in closed form. Explicit formulas for the Greeks in the case of the log-normal (GBM) asset price model were derived in previous chapters. In general, we rely on numerical methods. The integral in (18.32) can be computed with the use of quadrature formulae like those of Subsection 18.1.1:

$\upsilon (\tau ,S)\approx {\text{e}}^{-rT}{\displaystyle \sum _{i=1}^n{\mathcal{w}}_i\Lambda (s_i)\tilde{p}(\tau ;S,s_i),}$

where the summation is taken over a set of nodes s1, s2,..., sn ∈ ℝ+. The integral in (18.32) is an improper one. To use a quadrature formula, we can first apply a change of variables to transform the integration interval into a finite one. Alternatively, we can truncate the integration region and integrate the function Λ(S′)$\tilde{p}$(τ ; S, S′) w.r.t. S′ over a bounded interval (0,Smax). One case is where the payoff Λ is nonzero on a finite interval such as the put option payoff Λ(S′) = (K − S′)+ = (K − S)?{S′≤K}. If Λ(S′) = 0 for all S′ ≥ K, then we take Smax = K. In many other cases we may observe that the integrand goes to zero exponentially fast as S′→ ∞. If this is the case, then for any given ε > 0 one can find Smax > 0 such that ${\displaystyle {\int }_{S_{\max }}^{\infty }\Lambda (S^{\prime })\tilde{p}(\tau ;S,S^{\prime })\text{d}{S}^{\prime }<\epsilon }$.

Let us calculate the initial (set t = 0) value of a European call option under two asset price models that are more involved than the GBM model, namely, the constant elasticity of variance (CEV) diffusion model and the geometric variance gamma (GVG) model. First, we need to identify the transition probability density of each asset price process. Second, we apply a quadrature rule to evaluate the integral (18.32), which can be written for the case with the European call option with time to maturity T , strike price K, and spot S0 as

$\begin{array}{ll}\upsilon (T,S_0)\hfill & ={\text{e}}^{-rT}{\displaystyle {\int }_0^{\infty }\Lambda (S)\tilde{p}(T;S_0,S)\text{d}S\text{ (18.34)}}\hfill \\ \hfill & ={\text{e}}^{-rT}{\displaystyle {\int }_K^{\infty }(S-K)\tilde{p}(T;S_0,S)\text{d}S.\text{ (18.35)}}\hfill \end{array}$

In the following calculations, we assume that the risk-free interest rate r is 5%, the dividend yield q is 0, the initial asset price S0 is 100, and the time to maturity T is 0.5 (years).

Example 18.3

(Option Pricing in the CEV model). Consider the CEV diffusion model {S(t) ∈ ℝ+}t≥0, which follows

$\text{d}S(t)=\nu S(t)\text{d}t+\alpha S{(t)}^{\beta +1}\text{d}\tilde{W}(t)\text{ (18.36)}$

with real constants α > 0, β, ν = r − q and where ${\left\{\tilde{W}\left(t\right)\right\}}_{t\ge 0}$ is a standard Brownian motion in a risk-neutral measure ℙ. The local volatility is a nonlinear (power) function: σ(S)/S = αSβ. We consider the case with β < 0. Hence, the origin is attainable in finite time, i.e., the stock can attain a zero price in finite time. In particular, the origin is an exit boundary for β ∈ [−1/2, 0) and it is a regular boundary (which we specify as killing) for β < − 1/2. Note that the discounted process {e−νtS(t)}t≥0 is ℙ-martingale. For numerical illustration, we take the following parameter values: β =−2 and α = 2500. The instantaneous volatility at S0 = 100 is $\alpha S_0^{\beta }=\frac{2500}{{100}^2}=25\%$. The transition density has the known explicit representation

$\tilde{p}(t;S_0,S)={\text{e}}^{-\nu t}\frac{{({\text{e}}^{-\nu t}S)}^{-2\beta -\frac{3}{2}}{S}_0^{\frac{1}{2}}}{{\alpha }^2|\beta |{\tau }_t}{\text{e}}^{-\frac{{({\text{e}}^{-\nu t}S)}^{-2\beta }+S_0^{-2\beta }}{2{\alpha }^2{\beta }^2{\tau }_t}}{I}_{\frac{1}{2|\beta |}}\left(\frac{{({\text{e}}^{-\nu t}S)}^{-\beta }{S}_0^{-\beta }}{{\alpha }^2{\beta }^2{\tau }_t}\right),\text{ (18.37)}$

for S0,S, t > 0, ${\tau }_t=\frac{1}{2\nu \beta }(e^{2\nu \beta t}-1)$ if ν ≠ 0 and τt = t if ν = 0. This PDF involves Iα(z) which is the modified Bessel function of the first kind of order α and argument z.

To simplify our calculations, we recall how the above CEV process is related to the simpler known squared Bessel process (SQB){X(t) ∈ ℝ+}t≥0 with index μ≡ 1/(2β) and having the transition PDF

$p_u(t;x_0,x)={\left(\frac{x}{x_0}\right)}^{\frac{\mu }{2}}\frac{{\text{e}}^{-(x+x_0)/2t}}{2t}{I}_{|\mu |}\left(\sqrt{x{x}_0}/t\right).\text{ (18.38)}$

By defining the invertible smooth mapping X: ℝ+→ ℝ+

${\rm X}(S):=S^{-2\beta }/({\alpha }^2{\beta }^2)=S^{2|\beta |}/({\alpha }^2{\beta }^2),$

with inverse X− 1(x)≡ F(x) := (α2β2x)1/2|β|, we obtain the driftless, i.e., ν = 0, CEV process {S(0)(t) := F(X(t)), t ≥ 0}. The CEV processes with and without drift are simply related by a scale and time change:

$S(t)={\text{e}}^{\nu t}{S}^{(0)}({\tau }_t)={\text{e}}^{\nu t}\text{F}(X({\tau }_t)).\text{ (18.39)}$

It can be verified that the risk-neutral transition PDF of the CEV process, $\tilde{p}$, is related to that of the SQB process as

$\tilde{p}(t;S_0,S)={\text{e}}^{-\nu t}|{{\rm X}}^{\prime }({\text{e}}^{-\nu t}S)|p_{\mu }({\tau }_t;{\rm X}(S_0),{\rm X}({\text{e}}^{-\nu t}S)).\text{ (18.40)}$

By substituting the expression in (18.38) for pμ, we see that (18.40) is exactly Equation (18.37).

Using (18.40), for t = T and ν = r, and changing integration variables with S = erT F(x), x = X(e−rT S), and $\text{d}S=\frac{{\text{e}}^{rT}}{{{\rm X}}^{\prime }(\text{F}(x))}\text{d}x,$, the integral in (18.35) takes the form

$\upsilon (T,S_0)={\displaystyle {\int }_{{\rm X}({\text{e}}^{-rT}K)}^{\infty }(\text{F}(x)-{\text{e}}^{-rT}K)p_u({\tau }_T;{\rm X}(S_0),x)\text{d}x,\text{ (18.41)}}$

where $\mu =\frac{1}{2\beta }$. This is an improper integral over a half-line which can be transformed into an integral over a finite interval by a change of variables. For example,

${\displaystyle {\int }_a^{\infty }f(x)\text{d}x}={\displaystyle {\int }_0^{1}f\left(a+\frac{y}{1-y}\right)\frac{\text{d}y}{{(1-y)}^2}}.\text{ (18.42)}$

To compute the call price, the integral in (18.41) is transformed as in (18.42) where

$a={\rm X}({\text{e}}^{-rt}K)\text{ and }f(x)=\left(\text{F}(x)-{\text{e}}^{-rT}K\right)p_{\mu }({\tau }_T;{\rm X}(S_0),x).$

Alternatively, we can first evaluate the corresponding European put option (with same strike and maturity) and then obtain the call option value with the use of the put-call parity. The evaluation of a put option reduces to computing the integral (18.34) with Λ(S) = (K − S)+ over a finite interval (0, K).

With the above chosen model parameters, we employ Romberg integration with 101 initial nodes and four improvements. First, we evaluate the call option with strike K = 100. The output of the Romberg method is in shown Table 18.3, where the best approximation for the definite integral is the bottom rightmost boxed entry. The approximation error is given by the magnitude of the difference between the diagonal elements in the two bottom rows of the table. Second, we evaluate the call option for strikes K ∈ {90, 100, 110} and compare the results obtained with the call option values calculated under the GBM model with volatility σ = 0.25 (see Table 18.4). For comparison, we also calculated the Black–Scholes implied volatility for the CEV call option values. As is seen, the implied volatility is a decreasing function of strike K.

Example 18.4

(Pricing in the GVG model). The variance gamma (VG) process {X(t)}t≥0 is obtained by evaluating a scaled Brownian motion with drift W (θ,σ)(t) at a random time given by a gamma process G(t):

$X(t):=W^{(\theta ,\sigma )}(G(t))=\theta G(t)+\sigma W(G(t)),t\ge 0.$

Recall that the gamma process is characterized by two positive parameters, μ and υ. The parameter μ controls the rate of jump arrivals and the scaling parameter υ inversely controls the jump size. It is assumed that the process starts from 0 at time 0. The marginal distribution of a gamma process at time t is a gamma distribution with shape parameter $\frac{\mu t}{\upsilon }$ and rate parameter $\frac{1}{\upsilon }$. The mean and variance of G(t) are, respectively, $\frac{\mu t}{\upsilon }$ and $\frac{\mu t}{{\upsilon }^2}$. In the VG model, the parameter μ is set to 1 thus G(t) ~ Gamma $(\frac{t}{\upsilon },\frac{1}{\upsilon })$.

The risk-neutral asset price dynamics is given by the geometric variance gamma (GVG) process as

$S(t)=S_0{\text{e}}^{(r-q+w)t+X(t)},$

where $\omega =\frac{1}{\upsilon }\ln (1-\theta \upsilon -{\sigma }^2\upsilon /2)$ is a correction constant ensuring that the mean rate of return on the asset is r − q. The risk-neutral transition probability distribution is a mixture of the normal and gamma probability distributions. The PDF of Z(t) := ln(S(t)/S(0)) is given by

$p_{z(t)}(z)=\frac{2{\text{e}}^{\theta x/{\sigma }^2}}{{\upsilon }^{t/\upsilon }\sqrt{2\pi }\sigma \text{Γ}(t/\upsilon )}{\left(\frac{x^2}{2{\sigma }^2/\upsilon +{\theta }^2}\right)}^{\frac{t}{2\upsilon }-\frac{1}{4}}{K}_{\frac{t}{\upsilon }-\frac{1}{2}}\left(\frac{|x|}{{\sigma }^2}\sqrt{(2{\sigma }^2/\upsilon +{\theta }^2)}\right),\text{ (18.43)}$

for all t > 0, z ∈ ℝ, where Γ is the gamma function, Kα is the modified Bessel function of the second kind of order α, and x = z − (r − q)t − ωt. Since the asset price is given by an exponential of the log price ratio Z(t), i.e., S(t) = S0eZ(t), the value of a European option maturing at time T with payoff Λ(S) is given by

$\upsilon (T,S_0)={\text{e}}^{-rT}{\displaystyle {\int }_{-\infty }^{\infty }\Lambda (S_0{\text{e}}^z)p_{z(T)}(z)\text{d}z.\text{ (18.44)}}$

For a standard European call with strike K we have

$\upsilon (T,S_0)={\text{e}}^{-rT}{\displaystyle {\int }_{\ln (K/S_0)}^{\infty }(S_0{\text{e}}^z-K)p_{z(T)}(z)\text{d}z.}$

By using the transformation in (18.42), we can change the semi-infinite interval into the finite interval (0, 1) and then apply Romberg integration to evaluate the European call with K = 100. For numerical illustration, we take parameter values from Madan et al. (1998): θ = 0.1436, σ = 0.12136, and υ = 0.3. The other parameters are the same as those in the previous example. The output of Romberg method is given in Table 18.5, where the best approximation for the definite integral is the bottom rightmost boxed value.

18.2.2 Pricing European Options by the Monte Carlo Method

Let{S(t)}t≥0 be a stochastic asset price process starting at S0. Suppose that it can be simulated exactly from its transition distribution. The no-arbitrage price v(T, S0) of a European-style contract with the time to maturing T > 0 and payoff Λ(S(T)) can be estimated by the Monte Carlo method as described in Algorithm 18.1.

To illustrate the use of the Monte Carlo algorithm, we apply Algorithm 18.1 to the CEV diffusion asset price model, the geometric variance gamma pure jump model, and a multidimensional geometric Brownian motion model.

Example 18.5

(Pricing Options under the CEV Diffusion Model). The Monte Carlo simulation of the constant elasticity of variance (CEV) diffusion process that follows the SDE (18.36) is based on the reduction of its transition probability distribution to one of the randomized gamma distributions studied in Chapter 17. Assume the risk-neutral probability measure so that ν = r − q in (18.36). Recall from (18.39) that the CEV diffusion

Example 18.6

(Pricing Options under the Variance Gamma Model). The simulation of the VG process relies on sampling normal and gamma random variables. To simulate X(t) = B(G(t)) at some time t > 0, we first sample the gamma time G(t)~ Gamma(t/ν, 1/ν) and then sample a Brownian increment B(g)~ Norm(θg,σ2 g) given g = G(t). Thus, we can express the GVG asset price S(T) at maturity T (under the risk-neutral measure) as a function of independent gamma and normal random variables:

$S(T)\stackrel{d}{=}{S}_0{\text{e}}^{(r-q+w)T+\theta G/\nu +\sigma \sqrt{G/\nu }Z},\text{where }Z~Norm(0,1)\text{ and }G~Gamma(T/\nu ,1).$

In the Monte Carlo experiment, we estimate the value of a European call option with maturity $T=\frac{1}{2}$ and strike price K = 100. The GVG model parameters used here are the same as those in Example 18.4. In Table 18.8 we report the Monte Carlo estimate, standard error, and 95% confidence interval computed for several values of the sample size n. Figure 18.9 demonstrates the convergence of the sample estimate to the exact value of the call option as n increases.

Example 18.7

(Pricing Basket Options under the Multi-Asset GBM Model). Consider the pricing of derivative contracts under an N-asset GBM model with M Brownian factors. Let us recall the risk-free dynamics of the ith asset price described by

$S_i(t)=S_i(0)\exp \left(\left(r-q_i-{\sigma }_i^2/2\right)t+{\sigma }_i{\displaystyle \sum _{j=1}^M{\ell }_{ij}\tilde{W}(t)}\right),t\ge 0,$

where σi is the volatility of the ith asset, qi is the dividend yield on the ith asset, ${\left\{{\tilde{W}}_j\left(t\right)\right\}}_{t\ge 0}$ are i.i.d. Brownian motions (within the assumed risk-neutral measure $\tilde{\mathbb{P}}$) labelled by j = 1, 2,...,M, and the coefficients ℓij are chosen so that ${\displaystyle {\sum }_{j=1}^M{\ell }_{ij}^2=1}$ for all i = 1, 2,...,N. For each i = 1, 2,...,N, the log-return $X_i(t):=\ln \left(\frac{S_i(t)}{S_i(0)}\right)$ is a normal random variable with mean $(r-q_i-{\sigma }_i^2/2)t$ and variance ${\sigma }_i^{2}t$. The joint distribution of the log-returns X1(t),X2(t) ...,XN (t) is multivariate normal with N-by-N correlation matrix ρ = LLT, where L is an N-by-M matrix with entries ℓij. The covariance among the log-returns is given by the N-by-N matrix with entries

$\text{Cov}\left(X_i(t),X_j(t)\right)=t{\displaystyle {\Sigma }_{ij}}\text{ where }{\displaystyle \Sigma =\text{D}\rho \text{D} \text{and} \text{D}=\text{diag}({\sigma }_1,{\sigma }_2,...,{\sigma }_N).}$

Consider a multi-asset (basket) European derivative contract with time to maturity T and payoff Λ (S1(T),S2(T),...,SN (T) . We have priced some examples of such options, including options on the maximum with payoff maxi{wiSi(T)}, or the minimum with payoff mini{wiSi(T)}, etc., where w1, w2,..., wN are positive coefficients. Given time-0 spot values Si(0) = si > 0, i = 1, 2,...,N, the initial no-arbitrage value of this derivative is

$V(s_1,s_2,...,s_N;T)={\text{e}}^{-rT}{\tilde{\text{E}}}_0\left[\Lambda (S_1(T),S_2(T),...,S_N(T))|S_i(0)=s_i,1\le i\le N\right].$

For large values of N and M, we can only rely on Monte Carlo methods for pricing. To apply Algorithm 18.1, we need to generate n i.i.d. sample values of the payoff function. The completion of part 1 in Algorithm 18.1 involves the following three steps:

1. i. Sample n× M i.i.d. random variables $Z_j^{(k)}$, 1 ≤ k ≤ n, 1 ≤ j ≤ M, from the standard normal distribution.

1. ii. For each i = 1, 2,...,N, obtain n sample values of the terminal stock value Si(T),

   $S_i^{(k)}=s_i\exp \left((r-q_i-{\sigma }_i^2/2)T+{\sigma }_i\sqrt{T}{\displaystyle \sum _{j=1}^M{\ell }_{ij}{Z}_j^{(k)}}\right),k=1,2,...,n.$

1. iii. Calculate the sample values of the discounted payoff function

   ${\Lambda }_k={\text{e}}^{-rT}\Lambda \left(S_1^{(k)},S_2^{(k)},...,S_N^{(k)}\right),k=1,2,...,n.$

18.2.3 Pricing European Options by Tree Methods

Assume a log-normal model for stock prices in an economy with fixed interest rate r and zero stock dividend. Recall that under the risk-neutral measure $\tilde{\mathbb{P}}$ with $\tilde{W}$ as standard Brownian motion, the stock price has the representation

$S(t)=S_0{\text{e}}^{(r-{\sigma }^2/2)t+\sigma \tilde{W}(t)}.$

As was discussed in Chapter 2, the log-normal continuous-time model can be derived as a limiting case of discrete-time binomial tree models. Therefore, we can use a binomial tree to approximate the GBM dynamics. To guarantee a good order of approximation, the length of time periods has to be small enough. Multinomial tree models allow for improving both the accuracy and speed of the binomial method. In an n-nomial tree model, the stock price can jump to one of n possible branches over each time period. As a result, we can attain a higher order of approximation without having to increase the computational work by much.

18.2.3.1 Binomial Model

Consider a recombining binomial tree model with N + 1 trading dates tn = n δt, n = 0, 1,...,N uniformly spread throughout the interval [0, T] with step size $\delta t=\frac{T}{N}$ (i.e., the number of periods per year equals $\frac{1}{\delta t}$). The upward and downward factors, u and d, of the standard binomial tree model can be obtained by matching the first two moments of the stock price returns. The following solution is the most commonly used in binomial models:

$u=e^{\sigma \sqrt{\delta t}}\text{and }d=\frac{1}{u}=e^{-\sigma \sqrt{\delta t}}.\text{ (18.46)}$

The risk-neutral probability is given by the usual formula:

$\tilde{p}=\frac{{\text{e}}^{r\delta t}-d}{u-d}.\text{ (18.47)}$

For this choice of parameters, the variances of a single period return in the continuous- and discrete-time models agree up to $\mathcal{O}(\delta t^2)$. Note that other solutions for u, d, $\tilde{p}$ are possible; however, their analytical expressions are more cumbersome.

Consider the application of a binomial tree to pricing a European option with maturity T and payoff Λ(S). The continuous-time log-normal stock price process is approximated by a discrete-time binomial price process on the N-period lattice with the nodes Sn, k = S0undn − k, k = 0, 1,..., n and n = 0, 1,...,N. We calculate the option values Vn,k := V (n, Sn,k) using the backward-in-time recursion as follows.

1. At maturity, set VN, k = Λ(SN, k) for each k = 0, 1,...,N.
2. For each n = N − 1,N − 2,..., 0, compute the option values Vn, k as follows:

   $V_{n,k}={\text{e}}^{-r\delta t}(\tilde{p}{V}_{n+1,k+1}+(1-\tilde{p})V_{n+1,k}),k=0,1,...,n.$

Example 18.8

In this example, we study whether the binomial option price formula approximates the respective continuous-time solution with accuracy $\mathcal{O}(\delta t)$. Consider a butterfly option centred at 100 with the spread [80, 120]. The payoff is

$\Lambda (S)=\{\begin{array}{ll}S-80\hfill & \text{if }80<S\le 100,\hfill \\ 120-S\hfill & \text{if }100<S<120,\hfill \\ 0\hfill & \text{if }S\le 80\text{ or }S\ge 120.\hfill \end{array}\text{ (18.48)}$

Assume that the annual volatility σ is 0.25 and the annual risk-free rate of interest r is 0.05. Let us compute and compare the no-arbitrage prices and deltas of the butterfly option with $T=\frac{1}{2}$ for spot value S0 = 100 under the following asset price models:

1. (a) the log-normal (GBM) model,
2. (b) the binomial model with N = 50, 100, 200, 400 where u and d are given in (18.46).

The payoff of a butterfly option is equivalent to that of a portfolio in three standard European calls. In particular, the payoff in (18.48) is equivalent to a portfolio in one long call with strike 80, two short calls with strike 100, and one long call with strike 120. Therefore, the initial value and delta of a butterfly option with the payoff in (18.48) calculated under the log-normal model are, respectively,

$\begin{array}{l}V=C_{BS}(K=80)-2C_{BS}(K=100)+C_{BS}(K=120),\hfill \\ \Delta ={\Delta }_{BS}(K=80)-2{\Delta }_{BS}(K=100)+{\Delta }_{BS}(K=120),\hfill \end{array}$

where CBS(K) and ΔBS(K) are, respectively, the Black–Scholes call value and call delta for strike K and the same maturity and spot value as those used to price the butterfly option. For T = 0.5 and S0 = 100, we have the following value and delta of the butterfly option under the GBM model:

$V=7.97318602436266\text{ and }\Delta =-\mathrm{0.0381920926996022.}$

The errors of binomial approximations are reported in Table 18.10. As is seen from the table, the binomial model approximation error is proportional to δt.

18.2.3.2 Multinomial Models

In the binomial model, the asset price can jump to one of two possible states in each time step. To improve the accuracy of option valuation, we allow a jump to multiple states for the asset price. Consider a European-style option with a payoff function Λ(S) and maturity date T > 0. The no-arbitrage option value V (t, S) as a function of (actual) calendar time t and spot S satisfies

$V(t,S)={\text{e}}^{-r\delta t}\tilde{\text{E}}[V(t+\delta ,S(t+\delta t))|S(t)=S],0\le t<t+\delta t\le T.\text{ (18.49)}$

This equation can be used to recursively compute option prices

$V(T-\delta t,\cdot ),V(T-2\delta t,\cdot ),...,V(0,\cdot )$

at every time t = n δt starting at maturity T , where V (T, S) = Λ(S), and then advancing backward in time until time 0 is reached. To make this approach feasible, we need to introduce a discrete grid in the spatial variable S and then approximate the mathematical expectation in (18.49) by a finite sum evaluated on the grid.

Let us fix the time step δt and introduce the time grid tn = n δt for n = 0, 1,...,N, where $N=\frac{T}{\delta t}$. For every period [tn, tn+1], the ratio of asset prices S(tn+1) and S(tn) is a function of the Brownian increment $\delta {\tilde{W}}_{t_{n+1}}:=\tilde{W}\left(t_{n+1}\right)-\tilde{W}\left(t_n\right)~Norm\left(0,\delta t\right)$:

$\frac{S(t_{n+1})}{S(t_n)}={\text{e}}^{(r-{\sigma }^2/2)\delta t+\sigma \delta t{\tilde{W}}_{t_{n+1}}}.$

Each Brownian increment $\delta {\tilde{W}}_{t_n}$ is approximated by a discrete random variable $\delta {\tilde{W}}_n$ that takes values − b δx, − (b − 1) δ x,..., b δx with respective probabilities p− b, p− b+1,..., pb, summing to unity, where b ∈ ℕ is a branching parameter and δx > 0 is a mesh size. As a result, the Brownian motion ${\left\{\tilde{W}\left(t\right)\right\}}_{t\ge 0}$ is approximated by a Markov chain ${\left\{{\tilde{W}}_{t_n}\right\}}_{n\ge 0}$, which is defined at discrete times tn = n δt, n ≥ 0, and takes discrete values k δx, k ∈ ℤ. For each time tn, the variable ${\tilde{W}}_{t_n}$ is given by a sum of n i.i.d. increments,

${\tilde{W}}_{t_n}=\delta {\tilde{W}}_1+\delta {\tilde{W}}_2+\cdots +\delta {\tilde{W}}_n;$

it takes values in the set {k δx| k = − nb, − nb + 1 ,..., nb}. The value of ${\tilde{W}}_{t_{n+1}}$ conditional on ${\tilde{W}}_{t_n}=k\delta x$ is in the set {(k − b) δx, (k − b +1) δx,..., (k + b) δx}. The process ${\left\{{\tilde{W}}_{t_n}\right\}}_{n\ge 0}$ is a random walk whose realizations are paths in a recombining multinomial tree. The tree is organized such that there are 2bn + 1 nodes at time tn and each node has 2b + 1 branches. By setting b equal to 1 or 2, we can obtain a trinomial or pentanomial tree, respectively. The latter is shown in Figure 18.11.

The approximation of the process ${\left\{\tilde{W}\left(t\right)\right\}}_{t\ge 0}$ by ${\left\{{\tilde{W}}_{t_n}\right\}}_{n\ge 0}$ leads to an approximation of the continuous-time asset price process {S(t)}t≥0 by a discrete-time Markov chain ${\left\{S_{t_n}\right\}}_{n\ge 0}$ defined by

$S_{t_n}=\exp \left((r-{\sigma }^2/2)t_n+\sigma {\tilde{W}}_{t_n}\right).$

Denote the value of $S_{t_n}$ at the node (n, k) of the multinomial tree by

$S_{n,k}:=\exp ((r-{\sigma }^2/2)n\delta t+\sigma k\delta x).$

The conditional expectation in (18.49) can now be approximated as follows:

$\begin{array}{ll}{\text{e}}^{r\delta t}V(t_n,S_{n,k})\hfill & =\tilde{\text{E}}[V(t_{n+1},S(t_{n+1}))|S(t_n)=S_{n,k}]\hfill \\ \hfill & =\tilde{\text{E}}\left[V\left(t_{n+1},S_{n,k}{\text{e}}^{(r-{\sigma }^2/2)\delta t+\sigma \delta {\tilde{W}}_{t_{n+1}}}\right)\right]\hfill \\ \hfill & \approx \tilde{\text{E}}\left[V\left(t_{n+1},S_{n,k}{\text{e}}^{(r-{\sigma }^2/2)\delta t+\sigma \delta {\tilde{W}}_{t_{n+1}}}\right)\right]\hfill \\ \hfill & \approx {\displaystyle \sum _{j=-b}^b{p}_{j}V(t_{n+1,}{S}_{n+1,k+j}).}\text{ (18.50)}\hfill \end{array}$

This gives us the following computational scheme:

$V_{n,k}={\text{e}}^{-r\delta t}{\displaystyle \sum _{j=-b}^b{p}_j{V}_{n+1,k+j},n=N-1,N-2,...,0,k=-nb,-nb+1,...,nb}$

where Vn, k denotes the approximation of the option price V (tn,Sn,k) at the node (n, k). As usual, the option values at maturity T = N δt are given by the payoff function:

$V_{N,k}=\Lambda (S_{N,k}),k=-Nb,-Nb+1,...,Nb.$

To attain a high order of convergence in multinomial methods, the probabilities {pj; j = − b, − b + 1 ,..., b} should be chosen such that ${\tilde{W}}_{t_n}$ is a good approximation of $\tilde{W}\left(t_n\right)$. This can be achieved by using the moment matching method. According to Heston and Zhou (2000), if the first k moments of the Brownian increment δWn match those of $\delta {\tilde{W}}_{t_n}={\tilde{W}}_{t_n}-{\tilde{W}}_{t_{n-1}}$, that is,

${\displaystyle \sum _{j=-b}^b{(j\delta x)}^l{p}_j=0,\text{for all odd }l\le k\text{ (18.51)}}$

and

${\displaystyle \sum _{j=-b}^b{(j\delta x)}^l{p}_j=\frac{\delta t^{l/2}l!}{2^{l/2}(l/2)!},\text{for all even }l\le k,\text{ (18.52)}}$

and if the terminal payoff function is 2k times continuously differentiable, then the multi-nomial approximation (18.50) has a local error of $\mathcal{O}\left(\delta t^{\left(k+1\right)/2}\right)$, and the associated discrete-time solution converges to the continuous-time solution at a rate of $\mathcal{O}\left(\delta t^{\left(k-1\right)/2}\right)$, that is,

$V(t_n,S_{n,k})=V_{n,k}+\mathcal{O}\left(\delta t^{(k-1)/2}\right).$

We can easily match all the odd moments by putting p−j = pj for all j. To match the first b even moments, we must solve the simultaneous equations

${\displaystyle \sum _{j=0}^{b}2{(j\delta x)}^{2k}{p}_j=\frac{\delta t^k(2k)!}{2^{k}k!}},k=1,2,...,b.\text{ (18.53)}$

Additionally, we require the probabilities pj to sum up to one,

$p_0+2{\displaystyle \sum _{j=1}^b{p}_j=1.}\text{ (18.54)}$

The left- and right-hand sides of (18.53) are, respectively, proportional to δx2k and δtk. Thus, if assume that $\delta x=\sqrt{a\delta t}$ holds for some constant a > 0, then we can cancel out δx and δt in (18.53) to obtain

${\displaystyle \sum _{j=0}^{b}2{(ja)}^{2k}{p}_j=\frac{(2k)!}{2^{k}k!},k=1,2,...,b.\text{ (18.55)}}$

Solving these equations, we find the probabilities pk that do not depend on δt or δx. The parameter a can be chosen to match an additional moment constraint

${\displaystyle \sum _{j=1}^{b}2{(ja)}^{2b+2}{p}_j=\frac{(2(b+1))!}{2^{b+1}(b+1)!}\cdot \delta t^{b+1}.}\text{ (18.56)}$

Combining (18.54), (18.55), and (18.56), we obtain the following matrix equation on the probabilities p0, p1,..., pb and parameter $\ell =\frac{1}{a}$:

$\left[\begin{array}{ccccc}1& 2& 2& \cdots & 2\\ 0& 2& 2\cdot 2^2& \cdots & 2\cdot b^2\\ 0& 2& 2\cdot 2^4& \cdots & 2\cdot b^4\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 2& 2\cdot 2^{2b}& \cdots & 2\cdot b^{2b}\\ 0& 2& 2\cdot 2^{2(b+1)}& \cdots & 2\cdot b^{2(b+1)}\end{array}\right]\left[\begin{array}{c}{p}_0\\ p_1\\ p_2\\ ⋮\\ p_b\end{array}\right]=\left[\begin{array}{c}1\\ \delta t\\ 3\cdot {\ell }^2\\ ⋮\\ \frac{(2b)!}{2^{b}b!}\cdot {\ell }^b\\ \frac{(2b+2)!}{2^{b+1}(b+1)!}\cdot {\ell }^{b+1}\end{array}\right].\text{ (18.57)}$

For small values of b, the system (18.57) can be solved explicitly. For example, for b = 1, we obtain a trinomial tree with Brownian increments

$\delta {\tilde{W}}_n\stackrel{d}{=}\{\begin{array}{ll}\sqrt{a\delta t}\hfill & \text{with probability }{p}_1=\frac{1}{2a},\hfill \\ 0\hfill & \text{with probability }{p}_0=\frac{a-1}{a},\hfill \\ -\sqrt{a\delta t}\hfill & \text{with probability }{p}_{-1}=\frac{1}{2a},\hfill \end{array}$

where a = 3.

In Alford and Webber (2001), this equation was solved numerically for odd b. No feasible solutions appear to exist for even b. It was confirmed that multinomial methods for b = 3, 7, 11, 15, 19 have higher rates of convergence for European options than those of the binomial method. Note that as b increases the moment matching probabilities take values very close to those determined by the transition density function for Brownian increments:

$p_j\approx \frac{\delta x}{\sqrt{\delta t}}n\left(\frac{j\delta x}{\sqrt{\delta t}}\right)=\sqrt{a}n(j\sqrt{a}),$

where a can be found by solving (18.56). This approximate solution satisfies (18.55) with a great degree of accuracy when b is large and δx is small.

Example 18.9

In this example, we compute discrete-time approximations for the Black–Scholes price of a standard European call option with T = 0.5 (years), K = 100, S0 = 100, r = 0.05, q = 0, and σ = 0.25. The exact option value is 8.260014598601677. We apply three models, namely, the binomial, trinomial, and heptanomial tree models. In theory, the highest rate of convergence of discrete-time approximations to the continuous-time solution as the number of periods N goes to ∞ is provided by the heptanomial tree model. However, a higher rate of convergence is achieved for options with smooth payoffs. The derivative of a standard European payoff has a discontinuity at S = K. Thus, we should not expect that the trinomial and heptanomial tree models will outperform the binomial one. Indeed, as is demonstrated in Table 18.13a and Figure 18.12a, the rate of convergence of each of these three models is just linear. To study the rate of convergence of a tree model, we compute approximations of the option price for an increasing sequence of values of the number of periods N. Assuming that the approximation error is proportional to δtα = N−α, we perform a linear regression to estimate the exponent α. The results for the three models are reported in Table 18.13a and Figure 18.12a.

Example 18.10

(Soft-Strike Call Option under GBM). Power options differ from vanilla European options in that the payoff function is not linear but raised to some power function. In a previous chapter we learned how to analytically price such options under the GBM model. Typically, the terminal payoff of a power option is a quadratic function of the underlying. The widest possible application of power options is for addressing the nonlinear risk of option sellers. There was proposed a class of soft-strike options which do not have a single fixed strike price but a range of strike spread over an interval. Such options allow for addressing limitations of a standard delta hedging when the underlying asset is close to the strike at expiration of the option and any re-balancing of the hedging portfolio causes the underlying asset to become “pinned” to the strike price. Recall that this is related to the fact that the delta position for a standard call approaches a unit step (Heaviside) function as time to maturity goes to zero.

Consider the soft strike European call option with payoff

${\Lambda }_a(S)=\{\begin{array}{ll}0\hfill & \text{if }S<K-a,\hfill \\ \frac{1}{4a}{(S-K+a)}^2\hfill & \text{if }K-a\le S<K+a,\hfill \\ S-K\hfill & \text{if }S>K+a.\hfill \end{array}\text{ (18.58)}$

We have obtained an exact analytical (Black–Scholes) pricing formula for this option in a previous chapter, i.e., see equation (12.45). Note that in the limit as a → 0, the soft-strike payoff converges to the call payoff (S − K)+.

Since the payoff of a soft strike option is a continuously differentiable function of spot, the trinomial and heptanomial model will demonstrate a higher rate of convergence than that of the binomial tree model. Let us compute discrete-time approximations of the exact analytical price of the soft strike call option with K = 100 and a = 5. The other parameters are the same as those of the previous example. By evaluating the standard normal CDF expressions of the analytical pricing formula in (12.45), we obtain the exact option value (within 15 decimals) as 8.351245148287005. Notice that it is slightly larger than the value of a standard European call with the same maturity T and strike K from the previous example. This is expected since the payoff Λa(S) is strictly larger than the call payoff (S − K)+ for K − a < S < K + a. The errors for the multinomial approximations are reported in Table 18.13b and Figure 18.12b. Note that the rate of convergence of the binomial model is linear. The trinomial and heptanomial models demonstrate a superlinear rate of convergence.

18.2.4 Pricing European Options by PDEs

Let us assume the stock price process is a time-homogeneous diffusion with linear drift and local volatility function σ(S), i.e., S(t) satisfies the SDE

$\frac{\text{d}S(t)}{S(t)}=(r-q)\text{d}t+\sigma (S(t))\text{d}\tilde{W}(t).$

A standard European-style derivative has price V = V (t, S) as a function of calendar time t and spot S, where V satisfies the corresponding Black–Scholes PDE

$\frac{\partial V}{\partial t}+\frac{1}{2}{\sigma }^2(S)S^2\frac{{\partial }^{2}V}{\partial S^2}+(r-q)S\frac{\partial V}{\partial S}-rV=0,0\le t\le T,0<S<\infty ,\text{ (18.59)}$

subject to the terminal (payoff) condition at maturity as t ↗ T :

$V(T-,S)\equiv V(T,S)=\Lambda (S),0<S<\infty ,\text{ (18.60)}$

Here, Λ is the payoff function, r and q are, respectively, the risk-free interest rate and dividend rate, the volatility term is a state-dependent function expressed as a product of the spot and local volatility: σ(S)S. Note that this generalizes the case σ(S) = σ, where σ > 0 is the usual constant volatility parameter in the GBM model. Below, we focus on the case of pricing standard (vanilla) European options. The values of European call and put options, denoted below by VC and VP , respectively, satisfy the terminal conditions:

$\begin{array}{l}{V}_C(T,S)={(S-K)}^{+}\equiv \max \{S-K,0\},\text{ (18.61)}\hfill \\ V_P(T,S)={(K-S)}^{+}\equiv \max \{K-S,0\}.\text{ (18.62)}\hfill \end{array}$

The Black–Scholes PDE is defined on the domain

$D=\{(t,S)|0\le t\le T,0<S<\infty ,\}\text{ (18.63)}$

which is unbounded as S → ∞. To complete the formulation of the differential problem (18.59)–(18.60), we need to impose a boundary condition on the solution V for S = 0 and an asymptotic condition on V as S → ∞. These conditions depend on the option. For standard European options, we can use the put-call parity and the fact that an option becomes worthless when it is out of the money. The unboundedness of D makes the application of finite-difference numerical schemes unfeasible. One possible solution is to map the semi-infinite interval S ∈ (0, ∞) onto a finite one by applying a suitable change of variables. Another approach is to truncate D to obtain a bounded rectangular domain. Let us introduce a lower bound Smin that represents small asset prices and an upper bound Smax that represents large asset prices. The asset price S is allowed to change in the interval [Smin, Smax]. Consequently, the Black–Scholes PDE is considered in the truncated domain

$\widehat{D}=\{(t,S)|0\le t\le T,S_{\min }\le S\le S_{\max }\}.\text{ (18.64)}$

The selection of Smin and Smax is also dictated by the type of option. For knock-out barrier options we may set Smin = BL and/or Smax = BU , where BL and BU are, respectively, the lower and upper barriers.

We now need to specify the boundary conditions on V at S = Smin and S = Smax. Consider European options for which we set Smin close to zero and Smax very large. A European call option is out of the money (and hence worthless) when S = Smin (as long as Smin≈ 0). Therefore, for the call price VC we adopt the following boundary condition at the left boundary:

$V_C(t,S_{\min })=0\text{ for all }0\le t\le T.\text{ (18.65)}$

Similarly, a European put option is out of the money (and hence worthless) when S = Smax (as long as Smax ≫ K). Hence, we set

$V_P(t,S_{\max })=0\text{ for all }0\le t\le T.\text{ (18.66)}$

The boundary condition on the other side of $\widehat{D}$ can be obtained by considering the put-call parity

$V_C(t,S)-V_P(t,S)=S{\text{e}}^{-q(T-t)}-K{\text{e}}^{-r(T-t)}\text{ for all }0\le t\le T\text{ and }S>0.\text{ (18.67)}$

Combining (18.65) and (18.67) at S = Smin gives us the left boundary condition for the put value:

$V_P(t,S_{\min })=K{\text{e}}^{-r(T-t)}-S_{\min }{\text{e}}^{-q(T-t)}\text{ for all }0\le t\le T.\text{ (18.68)}$

Combining (18.66) and (18.67) at S = Smax gives us the right boundary condition for the call value:

$V_C(t,S_{\max })=S_{\max }{\text{e}}^{-q(T-t)}-K{\text{e}}^{-r(T-t)}\text{ for all }0\le t\le T.\text{ (18.69)}$

Since the payoff functions of a call and put option are linear in the underlying, we can also use the boundary conditions

$\frac{{\partial }^2}{\partial S^2}{V}_C(t,S_{\max })=0\text{ and }\frac{{\partial }^2}{\partial S^2}{V}_P(t,S_{\min })=0.\text{ (18.70)}$

18.2.4.1 Pricing by the Heat Equation

As we already saw, the Black–Scholes PDE (18.59) with constant volatility σ reduces to the classical heat equation (18.8) in the spatial variable x and time to maturity τ,

$\frac{\partial u(T,x)}{\partial T}=\frac{{\partial }^{2}u(T,x)}{\partial x^2},$

upon applying the following transformation to Equation (18.59):

$\begin{array}{l}S=K{\text{e}}^x,t=T-\frac{\tau }{{\sigma }^2/2},\hfill \\ V(t,S)=K{\text{ e}}^{-\gamma x-({\beta }^2+\ell )T}u(\tau ,x),\hfill \end{array}\text{ (18.71)}$

where the constants γ, β, and l are defined as follows:

$\begin{array}{ll}\gamma =\frac{1}{2}(k-1),\hfill & \beta =\frac{1}{2}(k+1)\hfill \\ \ell =\frac{q}{{\sigma }^2/2},\hfill & k=\frac{r-q}{{\sigma }^2/2}.\hfill \end{array}\text{ (18.72)}$

The change of variable in (18.71) maps the original domain D, defined in (18.63), and the truncated domain $\widehat{D}$, defined in (18.64), to

$\begin{array}{l}{D}_u=\{(\tau ,x)|0\le \tau \le \frac{{\sigma }^2}{2}T,-\infty <x<\infty \},\text{ (18.73)}\hfill \\ {\widehat{D}}_u=\left\{(\tau ,x)|0\le \tau \le \frac{{\sigma }^2}{2}T,x_{\min }\le x\le x_{\max }\right\},\text{ (18.74)}\hfill \end{array}$

respectively, where $x_{\min }:=\ln \left(\frac{S_{\min }}{K}\right)$ and $x_{\max }:=\ln \left(\frac{S_{\max }}{K}\right)$. The terminal conditions (18.61) and (18.62) are transformed into the initial conditions

$\begin{array}{ll}{u}_C(0,x)=\max \{{\text{e}}^{\beta x}-{\text{e}}^{\gamma x},0\}\hfill & (\text{for the call option}),\text{ (18.75)}\hfill \\ u_P(0,x)=\max \{{\text{e}}^{\gamma x}-{\text{e}}^{\beta x},0\}\hfill & (\text{for the put option}),\text{ (18.76)}\hfill \end{array}$

respectively, for all x ∈ [xmin, xmax]. Finally, the boundary conditions for the heat equation on the truncated domain ${\widehat{D}}_u$ are as follows:

$\begin{array}{ll}\{\begin{array}{l}{u}_C(\tau ,x_{\min })=0,\hfill \\ u_C(\tau ,x_{\max })={\text{e}}^{\beta x_{\max }+{\beta }^2\tau }-{\text{e}}^{\gamma x_{\max }+{\gamma }^2\tau }\hfill \end{array}\hfill & (\text{for the call option}),\text{ (18.77)}\hfill \\ \{\begin{array}{l}{u}_C(\tau ,x_{\min })={\text{e}}^{\gamma x_{\min }+{\gamma }^2\tau }-{\text{e}}^{\beta x_{\min }+{\beta }^2\tau },\hfill \\ u_P(\tau ,x_{\max })=0\hfill \end{array}\hfill & (\text{for the put option}),\text{ (18.78)}\hfill \end{array}$

for all $\tau \in [0,\frac{{\sigma }^2}{2}T]$. To find the numerical solution, we apply one of the finite-difference schemes to the initial boundary value problem obtained.

18.2.4.2 Pricing by the Black–Scholes PDE

Before we introduce a mesh in the domain $\widehat{D}$ and discretize the PDE (18.59), let us rewrite the terminal value problem (18.59)–(18.60) as an initial value problem. We replace the actual time t by the time to maturity τ = T − t in Equations (18.59), (18.60), (18.68), and (18.69) to obtain the following differential equation:

$\frac{\partial \upsilon }{\partial T}=\frac{1}{2}{\sigma }^2(S)S^2\frac{{\partial }^2\upsilon }{\partial S^2}+(r-q)S\frac{\partial \upsilon }{\partial S}-r\upsilon ,\mathcal{T}>0,S_{\min }<S<S_{\max }\text{ (18.79)}$

subject to the initial condition

$\upsilon (0,S)=\Lambda (S),\text{ (18.80)}$

and the boundary conditions

${\upsilon }_P(\tau ,S_{\min })=K{\text{e}}^{-r\tau }-S_{\min }{\text{e}}^{-q\tau }\text{ and }{\upsilon }_P(\tau ,S_{\max })=0\text{ (18.81)}$

and

${\upsilon }_C(\tau ,S_{\max })=S_{\max }{\text{e}}^{-q\tau }-K{\text{e}}^{-r\tau }\text{ and }{\upsilon }_C(\tau ,S_{\min })=0\text{ (18.82)}$

for the European put and call options, respectively. Here, the option value v = v(τ, S) is a function of the spot price S and the time to maturity τ.

Consider a uniform spacial mesh on the interval [Smin, Smax]:

$s_j=S_{\min }+j\delta s,j=0,1,...,n+1,\text{ where }\delta s=\frac{S_{\max }-S_{\min }}{n+1}.$

Replacing all derivatives with respect to S by their central finite-difference approximations, we obtain the following approximation to the Black–Scholes PDE (18.79):

$\begin{array}{ll}\frac{\partial \upsilon (\tau ,S)}{\partial \tau }\hfill & =\frac{1}{2}{\sigma }^2(S)S^2\frac{\upsilon (\tau ,S+\delta s)-2\upsilon (\tau ,S)+\upsilon (\tau ,S-\delta s)}{\delta s^2}\hfill \\ \hfill & +(r-q)S\frac{\upsilon (\tau ,S+\delta s)-\upsilon (\tau ,S-\delta s)}{2\delta s}-r\upsilon (\tau ,S)+\mathcal{O}(\delta s^2).\text{ (18.83)}\hfill \end{array}$

Let Vj(τ) denote the semi-discrete approximation to v(τ, sj). Applying (18.83) at each internal node sj, we obtain the following system of first-order ordinary differential equations:

$\begin{array}{l}\frac{d{V}_j\left(\tau \right)}{d\tau }=\underset{=:L_{j,j-1}}{\underbrace{\frac{1}{2}\left({\left(\frac{\sigma \left(s_j\right)s_j}{\delta s}\right)}^2-\frac{\left(r-q\right)s_j}{\delta s}\right)}}{V}_{j-1}\left(\tau \right)+\underset{=:L_{j,j}}{\underbrace{\left(-{\left(\frac{\sigma \left(s_j\right)s_j}{\delta s}\right)}^2-r\right)}}{V}_j\left(\tau \right)\hfill \\ +\underset{=:L_{j,j+1}}{\underbrace{\frac{1}{2}\left({\left(\frac{\sigma \left(s_j\right)s_j}{\delta s}\right)}^2+\frac{\left(r-q\right)s_j}{\delta s}\right)}}{V}_{j+1}\left(\tau \right),j=1,2,...,n\hfill \end{array}\text{ (18.84)}$

System (18.84) has n equations in n+2 unknown functions V0(τ),V1(τ),...,Vn(τ),Vn+1(τ). Using the boundary conditions we have the functions V0(τ)and Vn+1(τ), which, respectively, approximate the solution at the boundary nodes s0 = Smin and sn+1 = Smax. As a result, the system of differential equations (18.84) can be written as the following matrix-vector differential equation with an n-by-n tridiagonal coefficient matrix L whose entries are defined in (18.84):

$\frac{\text{dV}(\tau )}{\text{d}\tau }=\text{LV}(\tau )+\text{G}(\tau ),\text{ (18.85)}$

subject to the initial condition

$\text{V}(0)=\Lambda :={[\Lambda (s_1),\Lambda (s_2),...,\Lambda (s_n)]}^{{\rm T}}.\text{ (18.86)}$

Here we use the notation:

$\text{L}=\left[\begin{array}{cccccc}{L}_{11}& L_{12}& 0& \cdots & 0& 0\\ L_{21}& L_{22}& L_{23}& \cdots & 0& 0\\ 0& L_{32}& L_{33}& \cdots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& \cdots & L_{n-1,n-1}& L_{n-1,n}\\ 0& 0& 0& \cdots & L_{n,n-1}& L_{n,n}\end{array}\right],\text{V}(T)=\left[\begin{array}{c}{V}_1(T)\\ V_2(T)\\ ⋮\\ V_{n-1}\\ V_n(T)\end{array}\right].$

The vector G(T) ∈ ℝn is given by

${\left[\left(\frac{{\sigma }^2(s_0)s_0^2}{2\delta s^2}-\frac{(r-q)s_0}{2\delta s}\right)V_0(T),0,...,0,\left(\frac{{\sigma }^2(s_{n+1})s_{n+1}^2}{2\delta s^2}+\frac{(r-q)s_{n+1}}{2\delta s}\right)V_{n+1}(T)\right]}^{{\rm T}}.$

G(τ) contains boundary values of the mesh solution. For vanilla European options we can use the boundary conditions (18.81) and (18.82) to obtain the following:

$V_0(T)=K{\text{e}}^{-rT}-S_{\min }{\text{e}}^{-qT}\text{ and }{V}_{n+1}(T)=0\text{ for European put options}\text{ (18.87)}$

and

$V_{\mathrm{n+1}}(T)=S_{\max }{\text{e}}^{-qT}-K{\text{e}}^{-rT}\text{ and }{V}_0(T)=0\text{ for European call options}\text{ (18.88)}$

Applying finite differences in time, we derive the explicit method, the implicit method, and the Crank–Nicolson method as follows, by introducing the fully discretized approximate solution (the mesh solution) to the Black–Scholes equation:

$V_{i,j}\approx V_j(T_i)\approx V(x_j,T_i),i=0,1,...,m,j=0,1,...,n+1.$

The approximation solution at the ith time level will be denoted by

${\text{V}}^{(i)}:={[V_{i,1},V_{i,2},...,V_{i,n}]}^{{\rm T}}\approx \text{V}(t_i)$

for each i = 0, 1, 2,..., m.

The Explicit Method

Using a forward finite difference in (18.85) gives

$\frac{{\text{V}}^{(i+1)}-{\text{V}}^{(i)}}{\delta t}={\text{LV}}^{(i)}+\text{G}(t_i),$

which leads to the iterative rule

${\text{V}}^{(i+1)}=({\text{I}}_n+\delta t\text{L}){\text{V}}^{(i)}+\delta t\text{G}(t_i),i\ge 0.\text{ (18.89)}$

The Implicit Method

Approximating the time derivative in (18.85) with a backward finite difference gives

$\frac{{\text{V}}^{(i+1)}-{\text{V}}^{(i)}}{\delta t}={\text{LV}}^{(i+1)}+\text{G}(t_{i+1}),$

which leads to the iterative rule

$({\text{I}}_n-\delta t\text{L}){\text{V}}^{(i+1)}={\text{V}}^{(i)}+\delta t\text{G}(t_{i+1}),i\ge 0.\text{ (18.90)}$

The Crank–Nicolson Method

The Crank–Nicolson method is derived by taking the arithmetic average of the explicit (18.90) and the implicit (18.90) methods to give the following:

$\left({\text{I}}_n-\frac{\delta t}{2}\text{L}\right)V^{(i+1)}=\left({\text{I}}_n+\frac{\delta t}{2}\text{L}\right)V^{(i)}+\frac{\delta t}{2}(\text{G}(t_i)+\text{G}(t_{i+1})),i\ge 0.\text{ (18.91)}$

All the above schemes are subject to the initial condition

${\text{V}}^{(0)}=\Lambda$

with Λ defined in (18.86).

Example 18.11

In this example we apply the implicit and Crank–Nicolson schemes to price a standard call option under the Black–Scholes (log-normal) model. The parameters used are the same as those of Example 18.9. The truncated domain $\widehat{D}$ has the lower bound Smin = 0 and the upper bound

$S_{\max }=S_0\exp \left((r-q-{\sigma }^2/2)T+5\sigma \sqrt{T}\right)\cong \mathrm{244.3077.}$

To discretize the PDE problem, we split the space and time intervals in 5001 and 1001 subintervals, respectively. Thus, the time and space discretization steps are

$\delta t=\frac{0.5}{1001}\cong 0.0005\text{ and }\delta s\cong \frac{244.3077}{5001}\cong \mathrm{0.05.}$

The results for the two schemes are reported in Figure 18.14a. As is seen from the graphs in Figure18.14b, the rate of convergence of the Crank–Nicolson method is higher than that of the implicit method.

18.2.5 Calibration of Asset Price Models to Empirical Data

The term calibration refers to the computational process of fitting an asset price model to historical financial data. First, we select an asset price price model that seems to explain well empirical observations. There are two typical calibration procedures. The least squares method (LSM) allows for fitting model derivative values to the respective market prices. To apply this method, we assume that a derivative pricing formula is available in closed form or that the prices can be computed numerically in a fast and efficient manner. The maximum likelihood estimation (MLE) method allows for fitting the probability distribution to historical asset prices (or their returns). To apply it we assume that a transition probability density is known in closed form. Suppose the model can be characterized by a vector of parameters, say ξ = (ξ 1,ξ 2,...,ξ k). The objective is find the best-fitted vector $\widehat{\xi }$ that solves the respective optimization problem. Typically, there are two levels of calibration: first, an initial full calibration of all parameters of the model and, second, a faster re-calibration conditional on the vector $\widehat{\xi }$ computed previously that can be used as soon as new data have arrived. The second calibration scheme may be used throughout the day or even for longer periods, while the full calibration only needs to be executed at the outset or if markets move considerably.

18.2.5.1 Least Squares Method

Suppose we wish to calibrate an asset price model to historical European call and put option prices. Let the option with strike Ki and maturity Ti have an observed price Oi. Thus, we have two arrays of observations, i.e., a pair of variables (Ki,Ti) specifying the option and market prices Oi, where i = 1, 2,..., m. Let the model produce a price Ci = C(Ki,Ti; ξ ) for the same ith option. The goal of the calibration process is to minimize the least squares error between the model and market prices for the m option values:

$F(\xi ):={\displaystyle \sum _{i=1}^m{w}_i{(C(K_i,T_i;\xi )-O_i)}^2\to \underset{\xi }{\min },\text{ (18.92)}}$

where wi is a weight that reflects the relative significance of reproducing the ith option price precisely. The suitable choice of the weights wi, i = 1, 2,..., m, is important for good calibration results. The simplest choice is to set ${\omega }_i=\frac{1}{m}$. The confidence in individual data points is determined by the liquidity of the option. Thus, the weights can also be evaluated from the bid-ask spreads:

$w_i=\frac{1}{|O_i^{\text{ask}}-O_i^{\text{bid}}|}.$

Alternatively, as suggested by Cont and Tankov (2003), one may use the Black–Scholes (BS) “Vegas” evaluated at the implied volatilities of the market option prices to compute the weights:

$w_i={(\partial C^{BS}({\sigma }_i^{BS})/\partial \sigma )}^{-2},$

where ∂CBS/∂σ denotes the derivative of the BS option pricing formula with respect to the volatility σ, evaluated at $\sigma ={\sigma }_i^{\text{BS}}$, and ${\sigma }_i^{\text{BS}}={\sigma }_i^{\text{BS}}(O_i,K_i,T_i)$ is the BS implied volatility for the observed market price Oi. For the observed market prices Oi we may also use the midpoint of the bid and ask prices, i.e., $O_i=(O_i^{\text{ask}},O_i^{\text{bid}})/2$.

18.2.5.2 Maximum Likelihood Estimation

Suppose the asset price process {S(t)}t≽0 is a Markov process with a transition PDF pξ (u, t; S, S′). The PDF depends on the vector of parameters ξ . Typically, empirical data are available to us in the form of a truncated series of historical asset prices. Consider the array of observations (ti, ŝi), i = 0, 1, 2,..., m, where ŝi is a market asset (stock) price at calendar time ti. In general, the MLE method selects values of the model parameters that produce a probability distribution that gives the observed data the greatest probability. To use the method of maximum likelihood, we first need to specify the joint density function for all observations. Let the times ti form an increasing sequence. The joint probability density of the model asset prices at S(t0) = s0,S(t1) = s1,...,S(tm) = sm is then simply a product of transition PDFs

${\displaystyle \prod _{i=1}^m{p}_{\xi }(t_{i-1,}{t}_i;s_{i-1},s_i).}$

Now, by substituting the time series of observed asset prices s0 = ŝ0, s1 = ŝ1,..., sm = ŝm, we obtain the likelihood function

$L(\xi )={\displaystyle \prod _{i=1}^m{p}_{\xi }(t_{i-1,}{t}_i;{\widehat{s}}_{i-1},{\widehat{s}}_i).}$

The method of maximum likelihood estimates the best-fitted $\hat{\xi }$ by finding a value of ξ that maximizes L(ξ):

$\widehat{\xi }=\arg \max L(\xi ).\text{ (18.93)}$

In practice it is often more convenient to work with the (natural) logarithm of the likelihood function, called the log-likelihood. The result is the same regardless of whether we maximize the likelihood or the log-likelihood function, since the natural logarithm is a monotonically increasing function.

In solving the problems (18.92) and (18.93), we apply a numerical optimization routine such as the Nelder–Mead method (a nongradient based algorithm) or gradient search methods. In practice, the parameters ξ 1,ξ 2,...,ξ k have lower and/or upper constraints. However, by changing variables, a constrained optimization problem can be transformed into one without constrains.

18.3 Pricing Early-Exercise and Path-Dependent Options

Pricing exotic derivatives such as American and Asian options is a much less trivial problem in comparison to pricing standard European options. There are no closed-form formulae for the American and the arithmetic Asian options within the Black–Scholes model. In this section we focus on pricing options with an underlying asset price following the log-normal model or another Markov process whose transition distribution is known in closed form.

18.3.1 Pricing American and Bermudan Options

18.3.1.1 Pricing American Options by Tree Methods

Consider the application of the binomial tree method to pricing an American option with time to maturity T > 0 and having the intrinsic payoff value Λ(S). The continuous-time asset price process, which follows the log-normal model, is approximated by a discrete-time binomial model with N + 1 trading dates {n δt : n = 0, 1,...,N}, where δt is the length of one period and N δt = T . Let the factors u and d be specified as in (18.46). The one-period interest rate is erδt − 1. The risk-neutral probability $\tilde{p}$ of the upward move is given in (18.47).

Let Vn, k denote the American option value at time n δt for the asset price Sn, k = S0undn − k with k = 0, 1,..., n and n = 0, 1,...,N. From Chapter 7, we know that a backward-in-time recursion relation can be used to compute the option value at each node of the binomial tree as follows.

1. At maturity, set VN, k = Λ(SN, k) for each k = 0, 1,...,N.

2. For each n = N − 1,N − 2,..., 0 and k = 0, 1,..., n, the option value Vn, k is given by

$V_{n,k}=\max \{\Lambda (S_{n,k}),{\text{e}}^{-r\delta t}(\tilde{p}{V}_{n+1,k+1}+(1-\tilde{p})V_{n+1,k})\}.$

Let us see how the above algorithm can be used to compute prices of a Bermudan derivative that can only be exercised at discrete times 0 < t1 < t2 <··· < tM = T with M ≪ N. Assume that the exercise dates are all of the form n δt for some integer n. At each exercise date, the option value is a maximum of the continuation and the intrinsic values. In between the exercise dates, we treat the derivative as a European option. The steps are as follows.

1. At maturity, set VN, k = Λ(SN, k) for each k = 0, 1,...,N.

2. For each n = N − 1,N − 2,..., 0, compute the option values Vn, k, k = 0, 1,..., n, as follows:

(i) Compute the continuation values

$V_{n,k}^{\text{cont}}={\text{e}}^{-r\delta t}(\tilde{p}{V}_{n+1,k+1}+(1-\tilde{p})V_{n+1,k}).$

(ii) If n δt is an exercise date, i.e., n δt = tm for some m = 1, 2,...,M, then

$V_{n,k}=\max \{\Lambda (S_{n,k}),V_{n,k}^{\text{cont}}\};$

otherwise, we set

$V_{n,k}=V_{n,k}^{\text{cont}}.$

Example 18.12

In this example, we compute discrete-time approximations for the Black–Scholes price of American and Bermudan put options with T = 1, K = 95, S0 = 100, r = 0.05, q = 0, and σ = 0.25. As for the Bermudan options, we assume that the option can only be exercised at one of the M times $t_m=m\frac{T}{M}$, m = 1, 2,...,M. The results of our calculations are reported in Tables 18.15a and 18.15b.

18.3.1.2 Pricing Bermudan Options by the Monte Carlo Method

Instead of a continuously exercisable American option, we consider a Bermudan option that can only be exercised at a fixed set of times 0 < t1 < t2 <··· < tM = T . The option expires at time T and the payoff to the holder is Λ(S). The American option can be viewed as a limiting case of Bermudan options as M→ ∞ and max1≼k≼M|tk − tk − 1 | → 0. Suppose that the number M of exercise times is not too large. Even in this case the early-exercise option price problem poses a significant challenge to the Monte Carlo method.

The initial price of a Bermudan option can be calculated using a backward-in-time recursion. Let V (t, S) denote the no-arbitrage price at time t given that S(t) = S and the option has not been exercised before time t. It is more convenient to write the equation on the option values using discounted quantities. Consider the discounted payoff function Λk(S) = e−rtk Λ(S) at the exercise time tk. Let Vk(S) denote the present value of the time-tk option price. That is, Vk(S) = e−rtk V (tk,S) is the no-arbitrage option value calculated at time tk and discounted to time 0. We also denote Sk = S(tk) for all k = 0, 1,...,M. At the maturity time T = tM we have VM (S) = ΛM (S). At any exercise time tk < tM the holder of the option can exercise the option and immediately receive the payoff (the intrinsic value), or continue to hold the option. In the latter case, the average proceeds to the holder are given by the expected time-tk+1 option value (the continuation value). Therefore, the discounted option values Vk satisfy the following dynamic programming problem:

$\begin{array}{ll}{V}_M(S_M)\hfill & ={\Lambda }_M(S_M)\text{ (18.94)}\hfill \\ V_k(S_k)\hfill & =\max \left\{{\Lambda }_k(S_k),\tilde{\text{E}}[V_{k+1}(S_{k+1})|S_k]\right\},0\le k\le M-1.\text{ (18.95)}\hfill \end{array}$

The mathematical expectations considered here are computed under the risk-neutral probability measure $\tilde{\mathbb{P}}$:

$\tilde{\text{E}}[V_{k+1}(S_{k+1})|S_k]={\displaystyle {\int }_0^{\infty }{V}_{k+1}(S_{k+1})\tilde{p}(t_{k+1}-t_k;S_k,S_{k+1})\text{d}{S}_{k+1}}.$

The main difficulty of the problem is the computation of the continuation value. Let us consider three Monte Carlo approaches to dealing with this issue. Note that the asset price process {S(t)}t≽0 does not necessarily follow geometric Brownian motion. The only assumption regarding the underlying asset process we will require is that the Markov chain S0→ S1→ S2→ ···→ SM can be simulated exactly from the transition probability distribution. Although we only discuss the single-asset case, all the methods presented below can also be used for the evaluation of multi-asset Bermudan options.

Stochastic Tree Method. For each tk, the continuation value $\tilde{\text{E}}\left[V_{k+1}\left(S_{k+1}\right)|S_k\right]$ can be approximated by the Monte Carlo method. First, sample b successors $S_{k+1}^1,S_{k+1}^2,...,S_{k+1}^b$ from the exact probability distribution of Sk+1 conditional on Sk. Second, approximate the mathematical expectation by the sample mean:

$\tilde{\text{E}}[V_{k+1}(S_{k+1})|S_k]\approx \frac{1}{b}{\displaystyle \sum _{i=1}^b{V}_{k+1}(S_{k+1}^i)}.$

The main idea of the stochastic tree method is to construct a nonrecombining stochastic tree with M levels (so that each level corresponds to one of M exercise times) and b branches for each node. As a result we obtain a tree with bM branches at the last level (see Figure 18.16). Since the number of branches grows exponentially fast as M increases, this method is feasible only if both b and M are not very large. Algorithm 18.3 describes the construction of a stochastic tree in detail.

Stochastic Mesh Method. The main shortcoming of the stochastic tree method is that the number of branches of a stochastic tree grows exponentially-fast as the number M of exercise times increases. The stochastic mesh method proposed by Broadie and Glasserman (2004) also operates with a random tree, but there are two key distinctions. First, there is no branching at intermediate times and the total number of sample paths does not explicitly depend on the number of exercise times. We start with sampling n paths of the Markov chain S0→ S1→ S2→ ···→ SM , where Sk≡ S(tk) for all k. These paths form a stochastic tree with n branches (and n· M nodes) all originated at time 0. Second, in computing the continuation value at some time tk we use time-tk+1 option values from all branches regardless of whether or not the values are computed for the successor of the current node. Algorithm 18.4 describes the stochastic mesh method step by step. In comparison with the stochastic tree method, we additionally assume that the risk-neutral transition PDF $\tilde{p}\left(t;S,S^{\prime }\right)$ (under the risk-neutral measure) is known in closed form. Note that, for simplicity, we are assuming a time-homogeneous process.

To explain the formula for the weights, let us consider the computation of the continuation value in more detail. In computing the Monte Carlo approximation of the continuation value $\tilde{\text{E}}[V_{k+1}(S_{k+1})|S_k=S_k^i]$, the asset values $S=S_{k+1}^j,j\ne i$ at time tk+1 are not sampled from the single time-step transition PDF $\tilde{p}(t_{k+1}-t_k;S_k^i,S)$. Since all paths begin with S0, the values $S=S_{k+1}^j$ at time tk+1 are sampled from the marginal PDF $\tilde{p}$(tk+1; S0,S). At the node (tk,$S_k^i$) we have

$\tilde{\text{E}}[V_{k+1}(S_{k+1})|S_k=S_k^i]={\displaystyle {\int }_0^{\infty }{V}_{k+1}(S)\tilde{p}(T_{k+1}-t_k;S_k^i,S)\text{d}S}$

(multiply and divide the integrand by p(tk+1; S0,S))

$\begin{array}{l}={\displaystyle {\int }_0^{\infty }{V}_{k+1}(S)\frac{\tilde{p}(t_{k+1}-t_k;S_k^i,S)}{\tilde{p}(t_{k+1};S_0,S)}\tilde{p}(t_{k+1};S_0,S)\text{d}S}\hfill \\ =\tilde{\text{E}}\left[V_{k+1}(S_{k+1})\frac{\tilde{p}(t_{k+1}-t_k;S_k^i,S_{k+1})}{\tilde{p}(t_{k+1};S_0,S_{k+1})}\right],\hfill \end{array}$

where the time-tk+1 asset price Sk+1 has probability density $\tilde{p}\left(t_{k+1};S_0,S_{k+1}\right)$ in the last mathematical expectation. As we can see, in computing the option value at the node (tk,$S_k^i$) we use the option values from all nodes at the next time layer tk+1, as is illustrated in Figure 18.17, which gives a complete picture of all relations between the nodes of the stochastic tree. The tree looks more like a net or mesh. This structure explains the name of the stochastic mesh method.

Theorem 18.4

(Broadie and Glasserman (2004)). The mesh estimator $\hat{V}$0 is biased high, i.e., E[$\hat{V}$0] ≽ V (0,S0) for all n ≽ 1.

Proof. Let us prove the assertion by using Jensen’s inequality and induction. At the terminal time tM = T we have ${\widehat{V}}_M^i={\Lambda }_M(S_M^i)=V(T,S_M^i)$ for all i = 1, 2,..., n. Assume $E[{\widehat{V}}_{k+1}^i]\succcurlyeq V(t_{k+1},S_{k+1}^i)$ for some k = 0, 1,...,M − 1 and all i = 1, 2,..., n. At time tk we have

$\begin{array}{ll}\text{E}[{\widehat{V}}_k^i]\hfill & \text{E}\left[\max \left\{\Lambda (S_k^i),\frac{1}{n}{\displaystyle \sum _{j=1}^n\frac{\tilde{p}(t_{k+1}-t_k;S_k^i,S_{k+1}^j)}{\tilde{p}(t_{k+1};S_0,S_{k+1}^j)}{\widehat{V}}_{k+1}^j}\right\}\right]\hfill \\ \hfill & \succcurlyeq \max \left\{{\Lambda }_k(S_k^i),\text{E}\left[\frac{1}{n}{\displaystyle \sum _{j=1}^n\frac{\tilde{p}(t_{k+1}-t_k;S_k^i,S_{k+1}^j)}{\tilde{p}(t_{k+1};S_0,S_{k+1}^j)}{\widehat{V}}_{k+1}^j}\right]\right\}\hfill \\ \hfill & \succcurlyeq \max \left\{{\Lambda }_k(S_k^i),\text{E}\left[\frac{\tilde{p}(t_{k+1}-t_k;S_k^i,S{}_{k+1})}{\tilde{p}(t_{k+1};S_0,S_{k+1}^{})}V(t_{k+1},S_{k+1})\right]\right\}\hfill \\ \hfill & =\max \{{\Lambda }_k(S_k^i),V(t_k,S_k^i)\}\hfill \\ \hfill & =V(t_k,S_k^i).\hfill \end{array}$

In conclusion, we mention another popular Monte Carlo technique for pricing Bermudan options based on the least squares (regression) method by Longstaff and Schwartz (2001). Regression methods rely on the representation of the continuation value as a linear combination of some basis functions ψℓ : ℝ→ ℝ, ℓ = 1, 2,...,L. That is, the single time-step conditional expectation is written as a linear combination of a chosen set of so-called basis functions:

$\tilde{\text{E}}[V_{k+1}(S_{k+1})|S_k=S]={\displaystyle \sum _{\ell =1}^L{\beta }_{k\ell }{\psi }_{\ell }(S)}.$

18.3.2 Pricing Asian Options

Asian options are averaging options, where the payoff functions depend on some form of average of the underlying asset prices over the life of the option. By considering different types of averaging, such as arithmetic or geometric, continuous or discrete, one can obtain different types of options. The difficulty with the pricing of Asian options depends both on the underlying asset price model and on the type of the option. If we assume the asset price follows the GBM process, then one can derive closed-form pricing formulae for Asian options with geometric averaging. Pricing arithmetic Asian options is a computational challenge even under the GBM model.

18.3.2.1 Pricing Discrete-Time Asian Options by the Monte Carlo Method

Here we analyze only one type of averaging, namely, discrete arithmetic averaging at M time observation points equally spaced throughout [0,T], with random variables

$A_M:=\frac{1}{M}{\displaystyle \sum _{k=1}^M{S}_k,}{S}_k:=S\left(k\cdot \frac{T}{M}\right).\text{ (18.96)}$

By using methods of extrapolation, one can also approximate values of Asian options for continuous-time arithmetic averaging based on option values computed using discrete averaging. Most Asian options are of European style (i.e., not early exercise). There are two main examples of Asian options: the arithmetic average price (AAP) or floating price option and the arithmetic average strike (AAS) or floating strike option. The payoff functions of Asian calls (C) and puts (P) with time to maturity T are as follows:

$\begin{array}{ll}{C}_M^{AAP}:=\max \{A_M-K,0\},\hfill & P_M^{AAP}:=\max \{K-A_M,0\},\hfill \\ C_M^{AAS}:=\max \{S_M-A_M,0\},\hfill & P_M^{AAS}:=\max \{A_M-S_M,0\}.\hfill \end{array}\text{ (18.97)}$

The plain Monte Carlo evaluation of the no-arbitrage price of an Asian option with a payoff function of the form Λ(AM ,SM) is straightforward.

1. Generate n independent sample paths $(S_1^j,S_2^j,...,S_M^j)$, j = 1, 2,..., n, all starting at the initial spot S0.
2. For each path, calculate the arithmetic average

   $A_M^j=\frac{1}{M}{\displaystyle \sum _{i=1}^M{S}_i^j,j=1,2,...,n.}$

3. The no-arbitrage price $V_0=e^{-rT}\tilde{\text{E}}\left[\Lambda \left(A_M,S_M\right)\right]$ of the option is approximated by the discounted sample mean of the payoff values

   $V_0\approx {\overline{V}}_{0,n}={\text{e}}^{-rT}\frac{1}{n}{\displaystyle \sum _{j=1}^n\Lambda (A_M^j,S_M^j)}.$

4. Additionally, we can estimate the variance

   $\text{Var}\left({\text{e}}^{-rT}\Lambda (A_M,S_M)\right)\approx s_n^2={\text{e}}^{-2rT}\frac{1}{n}{\displaystyle \sum _{j=1}^n{\Lambda }^2(A_M^j,S_M^j)-{\overline{V}}_{0,n}^2}$

   to construct an asymptotically valid 100(1 − α)% confidence interval for V0:

   $({\overline{V}}_{0,n}-z_{\alpha /2}{s}_n/\sqrt{n},{\overline{V}}_{0,n}+z_{\alpha /2}{s}_n/\sqrt{n}).$

To implement the above algorithm, we only assume that sample paths can be generated exactly or approximately with high accuracy. Though the plain Monte Carlo approach is easy to implement, it is not very efficient. Several variance reduction methods have been proposed for valuing derivatives. In Boyle et al. (1997), the performances of the antithetic variate method, control variate (CV) method, and moment matching method were compared to that of the direct Monte Carlo method. They showed by simulation that the control variate method is the most efficient for valuing an Asian option. The CV method, considered in Boyle et al. (1997), was used to improve upon a naive Monte Carlo estimator for an Asian call option by choosing as analytical control variate the corresponding Asian call option with geometric averaging. We focus on this method of variance reduction.

We can achieve a larger variance reduction by considering multiple control variates. Suppose that there are two random variables Y (1) and Y (2) with respective means ν1 and ν2 known analytically. We construct the controlled Monte Carlo estimator of the mean of a random variable X, E[X] = μ, by using the random variable defined by

$Z=X+c_1({\nu }_1-Y^{(1)})+c_2({\nu }_2-Y^{(2)}).\text{ (18.98)}$

For all values of the parameters c1 and c2, the mean of Z coincides with the mean of X. The variance of the controlled estimator Z is

$\begin{array}{ll}\text{Var}(Z)\hfill & \begin{array}{l}=\text{Var}(X)-2c_1\text{Cov}(X,Y^{(1)})-2c_2\text{Cov}(X,Y^{(2)})\\ +c_1^2\text{Var}(Y^{(1)})+2c_1{\text{c}}_{\text{2}}\text{Cov}(Y^{(1)},Y^{(2)})+c_2^2\text{Var}(Y^{(2)})\end{array}\hfill \\ \hfill & \text{Var}(X)-2{\text{c}}^{{\rm T}}{\displaystyle {\sum }_{\text{X,Y}}}+{\text{c}}^{{\rm T}}{\displaystyle {\sum }_{{\text{X,Y}}^{\text{c}}}},\text{ (18.99)}\hfill \end{array}$

where

$\text{c}:=\left[\begin{array}{c}{c}_1\\ c_2\end{array}\right],{\displaystyle {\sum }_{X,\text{Y}}}:=\left[\begin{array}{l}\text{Cov}(X,Y^{(1)})\hfill \\ \text{Cov}(X,Y^{(2)})\hfill \end{array}\right],{\displaystyle {\sum }_{\text{Y},\text{Y}}}:=\left[\begin{array}{cc}\text{Var}(Y^{(1)})& \text{Cov}(Y^{(1)},Y^{(2)})\\ \text{Cov}(Y^{(1)},Y^{(2)})& \text{Var}(Y^{(2)})\end{array}\right].$

The optimal values of c1 and c2 that minimize Var(Z) can be found by solving simultaneous linear equations $\frac{\partial \text{Var}(Z)}{\partial c_i}=0$ for i = 1, 2:

$\{\begin{array}{l}{c}_1\text{Var}(Y^{(1)})+c_2\text{Cov}(Y^{(1)},Y^{(2)})=\text{Cov}(X,Y^{(1)})\hfill \\ c_1\text{Cov}(Y^{(1)},Y^{(2)})+c_2\text{Var}(Y^{(2)})=\text{Cov}(X,Y^{(2)})\hfill \end{array}\iff {\displaystyle {\sum }_{\text{Y,Y}}\text{C}}={\displaystyle {\sum }_{X,\text{Y}}}.\text{ (18.100)}$

The solution to (18.100) is the optimal vector

${\text{c}}_{\text{opt}}={\displaystyle {\sum }_{\text{Y,Y}}^{-1}}{\displaystyle {\sum }_{X,\text{Y}}}.\text{ (18.101)}$

Substituting it back into (18.99) gives

$\begin{array}{ll}\text{Var}(Z)\hfill & =\text{Var}(X)-2{\displaystyle {\sum }_{X,\text{Y}}^{{\rm T}}}{\displaystyle {\sum }_{\text{Y,Y}}^{-1}}{\displaystyle {\sum }_{X,\text{Y}}}+{\displaystyle {\sum }_{X,\text{Y}}^{{\rm T}}}{\displaystyle {\sum }_{\text{Y,Y}}^{-1}}{\displaystyle {\sum }_{\text{Y,Y}}}{\displaystyle {\sum }_{\text{Y,Y}}^{-1}}{\displaystyle {\sum }_{X,\text{Y}}}\hfill \\ \hfill & =\text{Var}(X)-{\displaystyle {\sum }_{X,\text{Y}}^{{\rm T}}}{\displaystyle {\sum }_{\text{Y,Y}}^{-1}}{\displaystyle {\sum }_{X,\text{Y}}}\hfill \\ \hfill & =\text{Var}(X)\left(-1\frac{{\displaystyle {\sum }_{X,\text{Y}}^{{\rm T}}}{\displaystyle {\sum }_{\text{Y,Y}}^{-1}}{\displaystyle {\sum }_{X,\text{Y}}}}{\text{Var}(X)}\right).\hfill \end{array}$

Here the ratio $\frac{{\displaystyle {\sum }_{X,\text{Y}}^{{\rm T}}}{\displaystyle {\sum }_{\text{Y,Y}}^{-1}}{\displaystyle {\sum }_{X,\text{Y}}}}{\text{Var}(X)}$ measures the strength of the linear relation between X and Y. It generalizes Corr2(X, Y) for X and scalar Y . As in the scalar case, the optimal vector of coefficients copt can be estimated by using sample estimates:

${\text{c}}_{\text{opt}}\approx {\text{S}}_{\text{Y,Y}}^{-1}{\text{S}}_{X\text{,Y}},$

where the (i, k)th entry of the square matrix SY,Y is

$\frac{1}{n}{\displaystyle \sum _{j=1}^n{Y}_j^{(i)}{Y}_j^{(k)}-{\overline{Y}}^{(i)}{\overline{Y}}^{(k)},}$

and the ith entry of the column vector SX,Y is

$\frac{1}{n}{\displaystyle \sum _{j=1}^n{X}_j{Y}_j^{(k)}-\overline{X}{\overline{Y}}^{(i)},}$

where (Xj,Y (1)j,Y (2)j), j = 1, 2,..., n are i.i.d. samples, and $\stackrel{\mathrm{—}}{X}$ and $\stackrel{\mathrm{—}}{Y}$ (i) are sample means given by

$\overline{X}=\frac{1}{n}{\displaystyle \sum _{j=1}^n{X}_j}\text{ and }{\overline{Y}}^{(i)}=\frac{1}{n}{\displaystyle \sum _{j=1}^n{Y}_j^{(i)}}i=1,2,$

respectively.

Let us consider the arithmetic average price call option with the payoff $C_M^{AAP}$ given by (18.97) within the GBM model. The variance of the naive Monte Carlo estimator can be significantly reduced with the use of multiple control variates. One control variate is the geometric Asian option, which can be evaluated analytically. The other control variate is an option whose payoff is an average of European payoffs with different maturity times. Let us find out how the payoffs of these control variates relate to the payoff of the average price call. According to the inequality of arithmetic and geometric means, the arithmetic mean of a set of nonnegative real numbers is greater than or equal to the geometric mean of the same set. That is,

$G_M\preccurlyeq A_M\text{ where }{G}_M:={\left({\displaystyle \prod _{k=1}^M{S}_k}\right)}^{\frac{1}{M}}.$

Using this and the fact that max{x, 0} is increasing in x, we have a lower bound for the payoff function:

$\max \{G_M-K,0\}\preccurlyeq \max \{A_M-K,0\}.\text{ (18.102)}$

If the asset price {S(t)}t≽0 follows GBM, then GM is log-normally distributed and hence the geometric option with payoff max{GM − K, 0} can be priced analytically. Let $C_0^{GAP}$ denote the no-arbitrage time-0 value of the geometric average price call with the payoff $C_M^{GAP}:=\max \{G_M-K,0\}$.

Since (x)+≡ max{x, 0} is a convex function of x, we have the following upper bound for $C_M^{AAP}$:

$\max \{A_M-K,0\}\preccurlyeq \frac{1}{M}{\displaystyle \sum _{k=1}^M\max \{S_k-K,0\},}\text{ (18.103)}$

which leads us to consider a new option with the payoff

$C_M^{AAC}:=\frac{1}{M}{\displaystyle \sum _{k=1}^M\max \{S_k-K,0\}}$

as a control variate for the average price call. As is readily observed, the payoff $C_M^{AAE}$ is just an average of European call payoff functions at exercise times tk, k = 1, 2,...,M. The price $C_0^{AAC}={\text{e}}^{-rT}\tilde{\text{E}}[C_M^{AAC}]$ can be computed numerically using only univariate integrals:

$C_0^{AAC}={\text{e}}^{-rT}\frac{1}{M}{\displaystyle \sum _{k=1}^M{\displaystyle {\int }_{-\infty }^{\infty }{\left({\text{e}}^{(r-q-{\sigma }^2/2)t_k+\sigma \sqrt{t_k}z}-K\right)}^{+}\frac{{\text{e}}^{-z^2/2}}{\sqrt{2\pi }}\text{d}z.}}\text{ (18.104)}$

As a result, we can improve the plain Monte Carlo algorithm, as follows.

1. Generate n independent sample paths $(S_1^j,S_2^j,...,S_M^j),j=1,2,...,n$, all starting at spot S0.
2. For each path, labelled j = 1, 2,..., n, calculate
   1. (a) the arithmetic average $A_M^j=\frac{1}{M}{\displaystyle \sum _{k=1}^M{S}_k^j}$,
   2. (b) the geometric average $G_M^j={\left({\displaystyle \prod _{k=1}^M{S}_k^j}\right)}^{\frac{1}{M}}$,
   3. (c) the arithmetic average price payoff $X_j=\max \left\{A_M^j-K,0\right\}$,
   4. (d) the geometric average price payoff $Y_j^{(1)}=\max \left\{G_M^j-K,0\right\}$,
   5. (e) the arithmetic average of European call payoffs $Y_j^{(2)}=\frac{1}{M}{\displaystyle \sum _{k=1}^M\max \{S_k^j-K,0\}}$.
3. Estimate the optimal coefficients c1 and c2: $\text{c}={\text{S}}_{\text{Y,Y}}^{\text{-1}}{\text{S}}_{X\text{,Y}}$, with the respective 2-by-2 matrix and 2-by-1 vector:

   $\begin{array}{l}{\text{S}}_{\text{Y,Y}}=\left[\begin{array}{ll}\frac{1}{n}{\displaystyle \sum _{j=1}^n{(Y_j^{(1)})}^2-{({\overline{Y}}^{(1)})}^2}\hfill & \frac{1}{n}{\displaystyle \sum _{j=1}^n{Y}_j^{(1)}{Y}_j^{(2)}-{\overline{Y}}^{(1)}{\overline{Y}}^{(2)}}\hfill \\ \frac{1}{n}{\displaystyle \sum _{j=1}^n{Y}_j^{(1)}{Y}_j^{(2)}-{\overline{Y}}^{(1)}{\overline{Y}}^{(2)}}\hfill & \frac{1}{n}{\displaystyle \sum _{j=1}^n{(Y_j^{(2)})}^2-{({\overline{Y}}^{(2)})}^2}\hfill \end{array}\right]\hfill \\ {\text{S}}_{X\text{,Y}}=\left[\begin{array}{l}\frac{1}{n}{\displaystyle \sum _{j=1}^n{X}_j{Y}_j^{(1)}-\overline{X}{\overline{Y}}^{(1)}}\hfill \\ \frac{1}{n}{\displaystyle \sum _{j=1}^n{X}_j{Y}_j^{(2)}-\overline{X}{\overline{Y}}^{(2)}}\hfill \end{array}\right].\hfill \end{array}$

4. Calculate the sample values of the controlled estimator

   $Z_j=X_j+c_1({\nu }_1-Y_j^{(1)})+c_2({\nu }_2-Y_j^{(2)}),j=1,2,...,n,$

   where ${\nu }_1=C_0^{GAP}$ and ${\nu }_2=C_0^{AAE}$.

5. The approximate no-arbitrage price $C_0^{AAP}$ of the AAP call option is then:

   $C_0^{AAP}\approx {\overline{C}}_{0,n}^{AAP}={\text{e}}^{-rT}\frac{1}{n}{\displaystyle \sum _{j=1}^n{Z}_j}.$

6. Estimate the variance

   $\text{Var}({\text{e}}^{-rT}{C}_M^{AAP})\approx s_n^2={\text{e}}^{-2rT}\frac{1}{n}{\displaystyle \sum _{j=1}^n{Z}_j^2-{\left({\overline{C}}_{0,n}^{AAP}\right)}^2}$

   and construct the confidence interval for $C_0^{AAP}$:

   $\left({\overline{C}}_{0,n}^{AAP}-z_{\alpha /2}{s}_n/\sqrt{n},{\overline{C}}_{0,n}^{AAP}+z_{\alpha /2}{s}_n/\sqrt{n}\right)$

   with an asymptotically valid 100(1 − α)% confidence level.

Example 18.13

Let us study the efficiency of the single and multiple control variate methods for Monte Carlo pricing of an AAP Asian call option. The underlying asset price follows the geometric Brownian motion model with σ = 0.25, r = 0.05, q = 0, T = 1. The spot value and strike price are, respectively, S0 = 100 and K = 100. The number of observation time moments is M = 128. Using the Black–Scholes pricing formula, we first compute the expected values of the control variates:

$C_0^{AAC}=6.755961403884071\text{ and }{C}_0^{GAP}=\mathrm{5.820772970849206.}$

After that, the Asian call option price is estimated by the Monte Carlo method using n = 10,000 sample paths. In Table 18.18 we report the simulation results for the crude estimate and three control variate estimates using, respectively, the geometric average price call (GAP), arithmetic average of European calls (AAC), and both averages as controls. The control variate method is very efficient. This observation is also supported by high correlation between the crude estimate and the GAP/AAC control estimates. The respective correlations coefficients are

${\rho }_{GAP}\approx 99.96\%,{\rho }_{AAC}\approx 99.46\%.$

References

J. Alford and N. Webber. Very high order lattice methods for one factor models. Available at SSRN 259478, 2001.

P. Boyle, M. Broadie, and P. Glasserman. Monte Carlo methods for security pricing. Journal of Economic Dynamics and Control, 21(8):1267–1321, 1997.

M. Broadie and P. Glasserman. A stochastic mesh method for pricing high-dimensional American options. Journal of Computational Finance, 7:35–72, 2004.

R. Cont and P. Tankov. Financial Modelling with Jump Processes. Chapman & Hall/CRC, 2003.

S. Heston and G. Zhou. On the rate of convergence of discrete-time contingent claims. Mathematical Finance, 10(1):53–75, 2000.

F. A. Longstaff and E. S. Schwartz. Valuing American options by simulation: A simple least-squares approach. Review of Financial Studies, 14(1):113–147, 2001.

D. B. Madan, P. Carr, and E. C. Chang. The Variance Gamma process and option pricing. European Finance Review, 2(1):79–105, 1998.

1Recall that the eigenvalues are zeros of the polynomial det(A − λIn).