7.3 Numerical Methods

The prices of many financial instruments do not have closed-form analytical solutions. The valuation of these financial instruments relies on a variety of numerical methods. In this section, we discuss the application of Monte Carlo simulation and finite difference methods.

Monte Carlo simulation

Monte Carlo simulation is a method for iteratively evaluating a deterministic model using random numbers with appropriate probabilities as inputs. For derivative pricing, it simulates a large number of price paths of the underlying assets with probability corresponding to the underlying stochastic process (usually under risk-neutral measure), calculates the discounted payoff of the derivative for each path, and averages the discounted payoffs to yield the derivative price. The validity of Monte Carlo simulation relies on the law of large numbers.

Monte-Carlo simulation can be used to estimate derivative prices if the payoffs only depend on the final values of the underlying assets, and it can be adapted to estimate prices if the payoffs are path-dependent as well. Nevertheless, it cannot be directly applied to American options or any other derivatives with early exercise options.

A. Explain how you can use Monte Carlo simulation to price a European call option?

Solution: If we assume that stock price follows a geometric Brownian motion, we can simulate possible stock price paths. We can split the time between $t$ and $T$ into $N$ equally-spaced time steps.¹⁰ So $\displaystyle\Delta t=\frac{T-t}{N}$ and $t_i=t+\Delta t\times i$, for $i=0,1,2,\cdots,N$. We then simulate the stock price paths under risk-neutral probability using equation $S_i=S_{i-1}e^{(r-\sigma^2/2)(\Delta t)+\sigma\sqrt{\Delta t}\varepsilon_i}$, where $\varepsilon_i$’s are IID random variables from standard normal distribution. Let’s say that we simulate $M$ paths and each one yields a stock price $S_{T,k}$, where $k=1,2,\cdots,M$, at maturity date $T$.

The estimated price of the European call is the present value of the expected payoff, which can be calculated as $\displaystyle C=e^{-r(T-t)}\frac{\sum_{k=1}^M\max(S_{T,k}-K,0)}{M}$.

B. How do you generate random variables that follow $N(\mu,\sigma^2)$ (normal distribution with mean $\mu$ and variance $\sigma^2$) if your computer can only generate random variables that follow continuous uniform distribution between 0 and 1?

Solution: This is a great question to test the basic knowledge of random number generation, the foundation of Monte Carlo simulation. The solution to this question can be dissected to two steps:

1. Generate random variable of $x\sim N(0,1)$ from uniform random number generator using inverse transform method and rejection method.
2. Scale $x$ to $\mu+\sigma x$ to generate the final random variables that follow $N(\mu,\sigma^2)$.

The second step is straightforward; the first step deserves some explanations. A popular approach to generating random variables is the inverse transform method: For any continuous random variable $X$ with cumulative density function $F$ ($U=F(X)$), the random variable $X$ can be defined as the inverse function of U: $X=F^{-1}(U),\ 0\le U\le1$. It is obvious that $X=F^{-1}(U)$ is a one-to-one function with $0\le U\le1$. So any continuous random variable can be generated using the following process:

• Generate a random number $u$ from the standard uniform distribution.
• Compute the value $x$ such that $u=F(x)$ as the random number from the distribution described by $F$.

For this model to work, $F^{-1}(U)$ must be computable. For standard normal distribution, $\displaystyle U=F(X)=\int_{-\infty}^X\frac{1}{\sqrt{2\pi}}e^{-x^2/2}dx$. The inverse function has no analytical solution. Theoretically, we can come up with the one-to-one mapping of $X$ to $U$ as the numeric solution of ordinary differential equation $\displaystyle F'(x)=f(x)=\frac{1}{\sqrt{2\pi}}e^{-x^2/2}$ using numerical integration method such as the Euler method.¹¹ Yet this approach is less efficient than the rejection method:

Some random variables have pdf $f(x)$, but no analytical solution for $F^{-1}(U)$. In these cases, we can use a random variable with pdf $g(y)$ and $Y=G^{-1}(U)$ to help generate random variables with pdf $f(x)$. Assume that $M$ is a constant such that $\displaystyle\frac{f(y)}{g(y)}\le M,\ \forall y$. We can implement the following acceptance-rejection method:

• Sampling step: Generate random variable $y$ from $g(y)$ and a random variable $v$ from standard uniform distribution $[0,1]$.
• Acceptance/rejection step: If $\displaystyle v\le\frac{f(y)}{Mg(y)}$, accept $x=y$; otherwise, repeat the sampling step.¹²

An exponential random variable ($g(x)=\lambda e^{-\lambda x}$) with $\lambda=1$ has cdf $u=G(x)=1-e^{-x}$. So the inverse function has analytical solution $x=-\log(1-u)$ and a random variable with exponential distribution can be conveniently simulated. For standard normal distribution, $\displaystyle f(x)=\frac{1}{\sqrt{2\pi}}e^{-x^2/2}$,

So we can choose $M=1.32$ and use the acceptance-rejection method to generate $x\sim N(0,1)$ random variables and scale them to $N(\mu,\sigma^2)$ random variables.

C. Can you explain a few variance reduction techniques to improve the efficiency of Monte Carlo simulation?

Solution: Monte Carlo simulation, in its basic form, is the mean of IID random variables $Y_1,Y_2,\cdots,Y_M:\ \displaystyle\overline Y=\frac1M\sum_{i=1}^MY_i$. Since the expected value of each $Y_i$ is unbiased, the estimator $\overline Y$ is unbiased as well. If $\operatorname{Var}(Y)=\sigma$ and we generate IID $Y_i$, then $\operatorname{Var}(\overline Y)=\sigma/\sqrt M$, where $M$ is the number of simulations. Not surprisingly, Monte Carlo 187 simulation is computationally intensive if $\sigma$ is large. Thousands or even millions of simulations are often required to get the desired accuracy. Depending on the specific problems, a variety of methods have been applied to reduce variance.

Antithetic variable: For each series of $\varepsilon_i$’s, calculate its corresponding payoff $Y(\varepsilon_1,\cdots,\varepsilon_N)$. Then reverse the sign of all $\varepsilon_i$’s and calculate the corresponding payoff $Y(-\varepsilon_1,\cdots,-\varepsilon_N)$. When $Y(\varepsilon_1,\cdots,\varepsilon_N)$ and $Y(-\varepsilon_1,\cdots,-\varepsilon_N)$ are negatively correlated, the variance is reduced.

Moment matching: Specific samples of the random variable may not match the population distribution well. We can draw a large set of samples first and then rescale the samples to make the samples’ moments (mean and variance are the most commonly used) match the desired population moments.

Control variate: If we want to price a derivative $X$ and there is a closely related derivative $Y$ that has an analytical solution, we can generate a series of random numbers and use the same random sequences to price both $X$ and $Y$ to yield $\hat X$ and $\hat Y$. Then $X$ can be estimated as $\hat X+(Y-\hat Y)$. Essentially we use $(Y-\hat Y)$ to correct the estimation error of $\hat X$.

Importance sampling: To estimate the expected value of $h(x)$ from distribution $f(x)$, instead of drawing $x$ from distribution $f(x)$, we can draw $x$ from distribution $g(x)$ and use Monte Carlo simulation to estimate expected value of $\displaystyle\frac{h(x)f(x)}{g(x)}$:

If $\displaystyle\frac{h(x)f(x)}{g(x)}$ has a smaller variance than $h(x)$, then importance sampling can result in a more efficient estimator. This method is better explained using a deep out-of-the-money option as an example. If we directly use risk-neutral $f(S_T)$ as the distribution, most of the simulated paths will yield $h(S_T)=0$ and as a result the estimation variance will be large. If we introduce a distribution $g(S_T)$ that has much wider span (fatter tail for $S_T$), more simulated paths will have positive $h(S_T)$. The scaling factor $\frac{f(x)}{g(x)}$ will keep the estimator unbiased, but the approach will have lower variance.

Low-discrepancy sequence: Instead of using random samples, we can generate a deterministic sequence of “random variable” that represents the distribution. Such low-discrepancy sequences may make the convergence rate $1/M$.

D. If there is no closed-form pricing formula for an option, how would you estimate its delta and gamma?

Solution: As we have discussed in problem A, the prices of options with or without closed-form pricing formulas can be derived using Monte Carlo simulation. The same methods can also be used to estimate delta and gamma by slightly changing the current underlying price from $S$ to $S\pm\delta S$, where $\delta S$ is a small positive value. Run Monte Carlo simulation for all three starting prices $S-\delta S$, $S$ and $S+\delta S$, we will get their corresponding option prices $f(S-\delta S)$, $f(S)$ and $f(S+\delta S)$.

Estimated delta: $\displaystyle\Delta=\frac{\delta f}{\delta S}=\frac{f(S+\delta S)-f(S-\delta S)}{2\delta S}$

Estimated gamma: $\displaystyle\Gamma=\frac{(f(S+\delta S)-f(S))-(f(S)-f(S-\delta S))}{\delta S^2}$

To reduce variance, it’s often better to use the same random number sequences to estimate $f(S-\delta S)$, $f(S)$ and $f(S+\delta S)$.¹⁴

E. How do you use Monte Carlo simulation to estimate $\pi$?

Solution: Estimation of $\pi$ is a classic example of Monte Carlo simulation. One standard method to estimate $\pi$ is to randomly select points in the unit square ($x$ and $y$ are independent uniform random variables between 0 and 1) and determine the ratio of points that are within the circle $x^2+y^2\le1$. For simplicity, we focus on the first quadrant. As shown in Figure 7.1, any points within the circle satisfy the equation $x_i^2+y_i^2\le1$. The percentage of the points within the circle is proportional to its area:

So we generate a large number of independent $(x,y)$ points, estimate the ratio of the points within the circle to the points in the square, and multiply the ratio by 4 to yield an estimation of $\pi$. Figure 7.1 uses only 1000 points for illustration. With today’s 189 computing power, we can easily generate millions of $(x,y)$ pairs to estimate $\pi$ with good precision. 1,000 simulations with 1,000,000 $(x,y)$ points each using Matlab took less than 1 minute on a laptop and gave an average estimation of $\pi$ as 3.1416 with standard deviation 0.0015.

Finite difference method

The finite difference method is another popular numerical technique for derivative pricing. It numerically solves a differential equation to estimate the price of a derivative by discretizing the time and the price of the underlying security. We can convert the Black-Scholes-Merton equation, a second order nonlinear partial differential equation, to a heat diffusion equation (as we did in Chapter 6). This new equation, expressed as a function of $\tau$ (time to maturity) and $x$ (a function of the price of the underlying security), is a general differential equation for derivatives. The difference between various derivatives lies in the boundary conditions. By building a grid of $x$ and $\tau$ and using the boundary conditions, we can recursively calculate $u$ at every $x$ and $\tau$ using finite difference methods.

A. Can you briefly explain finite difference methods?

Solution: There are several version of finite difference methods used in practice. Let’s briefly go over the explicit difference method, the implicit difference method and the 190 Crank-Nicolson method. As shown in Figure 7.2, if we divide the range of $\tau$, $[0,T]$, into $N$ discrete intervals with increment $\Delta\tau=T/N$ and divide the range of $x$, $[x_0,x_J]$, into $J$ discrete intervals with increment $\Delta x=(x_J-x_0)/J$, the time $\tau$ and the space of $x$ can be expressed as a grid of $\tau_n,\ n=1,\cdots,N$ and $x_j,\ j=1,\cdots,J$.

The explicit difference method uses the forward difference at time $\tau_n$ and the second-order central difference at $x_j$:

Rearranging terms, we can express $u_j^{n+1}$ as a linear combination of $u_{j+1}^n$, $u_j^n$ and $u_{j-1}^n$: $u_j^{n+1}=\alpha u_{j-1}^n+(1-2\alpha)u_j^n+\alpha u_{j+1}^n$, where $\alpha=\Delta t/(\Delta x)^2$. Besides, we often have boundary conditions $u_j^0$, $u_0^n$, and $u_J^n$ for all $n=1,\cdots,N;\ j=0,\cdots,J$. Combining the boundary 191 conditions and equation $u_j^{n+1}=\alpha u_{j-1}^n+(1-2\alpha)u_j^n+\alpha u_{j+1}^n$, we can estimate all $u_j^n$’s on the grid.

The implicit difference method uses the backward difference at time $t_{n+1}$ and the second-order central difference at $x_j$:

The Crank-Nicolson method uses the central difference at time $(t_n+t_{n+1})/2$ and the second-order central difference at $x_j$:

B. If you are solving a parabolic partial differential equation using the explicit finite difference method, is it worse to have too many steps in the time dimension or too many steps in the space dimension?

Solution: The equation for $u_j^{n+1}$ in the explicit finite difference method is $u_j^{n+1}=\alpha u_{j-1}^n+(1-2\alpha)u_j^n+\alpha u_{j+1}^n$, where $\alpha=\Delta t/(\Delta x)^2$. For the explicit finite difference method to be stable, we need to have $1-2\alpha>0\Rightarrow\Delta t/(\Delta x)^2\lt1/2$. So a small $\Delta t$ (i.e., many time steps) is desirable, but a small $\Delta x$ (too many space steps) may make $\Delta t/(\Delta x)^2>1/2$ and the results unstable. In that sense, it is worse to have too many steps in the space dimension. In contrast, the implicit difference method is always stable and convergent.