Recall that we model a stock price $S$ as a stochastic process $\{S_{t}:0\leq t\leq T\}$ governed by the stochastic differential equation (geometric Brownian motion)
$$(1)\;\;\;\;dS_{t}=\mu S_{t}dt+\sigma S_{t}dW_{t}$$
where $\mu$ is the drift or expected return on the stock, $W_{t}$ is a Weiner process $(dW_{t}=N(0,1)\sqrt{dt}$ where $N(0,1)$ is the standard normal distribution), and $\sigma$ is the volatility of the stock (a measure of the variance since $\sigma N(0,1)=N(\sigma^{2},1)$). Using arbitrage arguments (Merton's method) we arrive at the (linear) Black-Scholes-Merton partial differential equation
$$(2)\;\;\;\;V_{t}+\frac{1}{2}\sigma^{2}S^{2}V_{ss}+rSV_{s}-rV=0.$$ The price $V$ of any derivative must satisfy this equation, and conversely, any solution to (2) gives the price of some derivative based on the boundary conditions used (for European call options, the boundary conditions at the expiry time $T$ are $V(\cdot,T)=\max(S_{T}-K,0)$ where $K$ is the strike price of the option and $S_{T}$ is the value of the (random) stock price at time $T$). Pricing methods for derivatives based on (2) are known as PDE methods in mathematical finance, and usually consist of numerically solving (2) using finite difference methods (since the domain of definition is rectangular in the $(S,t)$ coordinate system). This can be difficult, however, since the varied boundary conditions must be taken into account and issues of convergence, stability, etc. enter in.
A different (and more recently popular) approach is probabilistic and involves sampling the randomness in the geometric Brownian motion model (1) of the stock price multiple times and taking an average. Indeed, an application of the Feynman-Kac formula shows that when discounted appropriately, solutions to (1) are martingales and this implies that $V$ is simply the expected value of the discounted payoff of the derivative. For a European call option, the payoff is $f(S)=\max(S-K,0)$ and so we get
$$(3)\;\;\;\;V_{\text{EuroCall}}=e^{-rT}\mathcal{E}[f(S_{T})]$$
where $r$ is the riskless return, $\mathcal{E}$ is taken under the risk-neutral probability measure, and $S_{T}$ is the final value of the stochastic process $\{S_{t}\}$ at expiry of the option $T$. (From now on, $V$ is for the price of a European call option.) Thus, beginning with (1), passing to the logarithm, applying Ito's Lemma and then applying (3) we get
$$\begin{align*}
&dS_{t}=rS_{t}\;dt+\sigma S_{t}dW_{t}\\
\Longrightarrow&d\log S_{t}=(r-\frac{1}{2}\sigma^{2})dt+\sigma dW_{t}\\
\Longrightarrow&\log S_{t}=\log S_{0}+(r-\frac{1}{2}\sigma^{2})t+\sigma W_{t}\\
\Longrightarrow&S_{t}=S_{0}\exp\left\{\left(r-\frac{1}{2}\sigma^{2}\right)t+\sigma\sqrt{T}N(0,1)\right\}\\
\Longrightarrow&V=e^{-rT}\mathcal{E}\left[f\left(S_{0}\exp\left\{\left(r-\frac{1}{2}\sigma^{2}\right)t+\sigma\sqrt{T}N(0,1)\right\}\right)\right]\\
(4) \Longrightarrow&V=e^{-rT}\mathcal{E}\left[\max\left(S_{0}\exp\left\{\left(r-\frac{1}{2}\sigma^{2}\right)t+\sigma\sqrt{T}N(0,1)\right\}-K,0\right)\right]\\
\end{align*}$$
(Recall that $dW_{t}=N(0,1)\sqrt{dt}$ and so $W_{T}=N(0,1)\sqrt{T}$). The above computations are somewhat complicated when carried out in detail, and while the explanation for how (2) is derived from (1) is relatively straight-forward, the derivation of (3) from (1) (and thus the previous computations) is significantly more subtle but far more important in the theory of mathematical finance. A full explanation of the concepts involved (properties of Weiner processes, stochastic integration and Ito's lemma, risk-neutral valuation and the risk-neutral measure, change of numeraire, etc.) is beyond the scope of this post but will be elaborated on in subsequent posts. In any event, we only need to take faith that (4) is equivalent to (2) for pricing European call options.
The idea of Monte Carlo simulation is now evident. We simulate $N$ paths of $S_{T}$ by sampling from the standard normal distribution $N(0,1)$ and then compute $f(S^{n}_{T})$. Since each sampling is independent, the random variables $\{f(S_{T}^{n})\}_{n}$ are independent and identically distributed and so the law of large numbers implies
$$\frac{1}{N}\sum_{n=1}^{N}f(S^{n}_{T})\to\mathcal{E}[f(S_{T})]\;\text{as}\;N\to\infty$$
pointwise (and of course, in probability). We then then discount at $e^{-rT}$ and this is our estimate on the price $V$.
The following implementation is coded procedurally in C++ and prices a European call option using the method explained above. The only part which may require explanation is the method in which $N(0,1)$ is sampled. While there are several ways to do this, a simple yet accurate/efficient method is known as the Box Muller algorithm, implemented here.
MonteCarloMethod_Source.cpp
#include <iostream>
#include <cmath>
using namespace std;
double BoxMullerGaussian();
double MonteCarloSimulator(double Expiry, double Strike, double Spot, double Vol, double Riskless, unsigned long N);
int main(void) {
double Expiry, Strike, Spot, Vol, Riskless, N = 0;
cout << "Enter: Expiry, Strike, Spot, Vol, Riskless, # Trials" << endl;
cin >> Expiry >> Strike >> Spot >> Vol >> Riskless >> N;
cout << "The price of the option is "
<< MonteCarloSimulator(Expiry, Strike, Spot, Vol, Riskless, N)
<< endl;
double pause;
cin >> pause;
return 0;
}
double BoxMullerGaussian() {
double x,y;
double sizeSquared;
do {
x = 2.0*rand()/(double)(RAND_MAX) - 1;
y = 2.0*rand()/(double)(RAND_MAX) - 1;
sizeSquared = x*x + y*y;
} while(sizeSquared >= 1.0);
return x*sqrt(-2*log(sizeSquared)/sizeSquared);
}
double MonteCarloSimulator(double Expiry, double Strike, double Spot, double Vol, double Riskless, unsigned long N) {
double var = Vol*Vol*Expiry;
double std = sqrt(var);
double itoCorrection = -0.5*var;
double movedSpot = Spot*exp(Riskless*Expiry + itoCorrection);
double _Spot = 0;
double runningSum = 0;
for (unsigned long i=0; i<N; i++) {
double _Gaussian = GetOneGaussianByBoxMuller();
_Spot = movedSpot*exp(std*_Gaussian);
double _Payoff = _Spot - Strike;
_Payoff = _Payoff>0 ? _Payoff : 0;
runningSum += _Payoff;
}
return exp(-Riskless*Expiry)*(runningSum/N); // expectation
}
As an example, we consider Microsoft's stock. Today it closed at 41.35, so $\text{Spot}=41.35$ (the spot price in option pricing is $S_{0}$, which is not random, but $S_{t}$ is for $t>0$ upto $t=T=\text{Expiry}$). We write the option with a strike price $\text{Strike}=40.00$ so that the option is in the money initially (for the buyer). If the time to expiry is 6 mo., then $\text{Expiry}=0.5$. The 6mo. risk-free rate is $\text{Riskless}=0.33$ (taken to be the LIBOR rate quoted this week). The implied 6mo. volatility for Microsoft options is around $\text{Vol}=0.22$. With these parameters and $N=100$ trials, we get
$$V_{\text{MSFTEuroCall}}=8.98$$
