Monte Carlo Pricing of European Options

From a Risk-Neutral Price to a Simulation

A European option pays a function of the asset price at one fixed maturity. For a call with strike \(K\), the payoff is \((S_T-K)^+\); for a put, it is \((K-S_T)^+\), where \(x^+=\max(x,0)\). Under the risk-neutral measure, the option price is the expected payoff discounted to the present. If the asset follows geometric Brownian motion, its terminal value can be simulated directly from a normal random variable. No time-stepping approximation is needed to generate that value.

Work on a filtered probability space \((\Omega,\mathcal{F},(\mathcal{F}_t)_{t\ge0},\mathbb{P})\) carrying a Brownian motion. For pricing, assume a risk-neutral measure \(\mathbb{Q}\) under which \(B_t\) is Brownian motion, the interest rate \(r\) is constant, the asset pays no dividends, and its dynamics are \(dS_t=rS_t\,dt+\sigma S_t\,dB_t\), with \(S_0>0\) and constant volatility \(\sigma>0\). Here \(rS_t\) is the drift term and \(\sigma S_t\) is the diffusion coefficient. The probability measure \(\mathbb{P}\) used to describe real-world outcomes and the pricing measure \(\mathbb{Q}\) need not be the same.

Whenever an Ito integral is used below, its integrand is progressively measurable and belongs to \(\mathcal{L}^2\) on each finite time interval: \(\mathbb{E}^{\mathbb{Q}}\int_0^T H_s^2\,ds<\infty\). For geometric Brownian motion this condition holds for \(H_s=\sigma S_s\), since \(\mathbb{E}^{\mathbb{Q}}[S_s^2]=S_0^2e^{(2r+\sigma^2)s}\), whose integral over a finite interval is finite.

The Exact Terminal Distribution

For the simulation, it matters that the terminal formula is the solution of the SDE, not merely a process that happens to satisfy the same differential equation. The following argument establishes the formula and identifies the uniqueness condition that connects it to \(S_t\).

Proof. Define \(f(t,x)=S_0\exp((r-\tfrac12\sigma^2)t+\sigma x)\) and \(X_t=f(t,B_t)\). Its derivatives are \(f_t=(r-\tfrac12\sigma^2)f\), \(f_x=\sigma f\), and \(f_{xx}=\sigma^2f\). Ito's formula, including its second-derivative correction, gives

Also \(X_0=S_0\). The coefficients \(b(x)=rx\) and \(a(x)=\sigma x\) are globally Lipschitz, since \(|b(x)-b(y)|=|r||x-y|\) and \(|a(x)-a(y)|=\sigma|x-y|\), and have linear growth. The standard existence-and-uniqueness theorem for SDEs therefore gives pathwise uniqueness: any two solutions driven by the same Brownian motion with the same initial value agree almost surely at every time. Since \(X\) is a solution and \(S\) is the solution under consideration, this uniqueness implies \(S_t=X_t\). Finally, \(B_T\) has the same distribution as \(\sqrt{T}Z\), giving the stated terminal law. \(\square\)

The normal random variable in this representation is the only random input needed for one terminal simulation. For independent standard normal draws \(Z_1,\ldots,Z_N\), compute \(S_T^{(i)}=S_0\exp((r-\tfrac12\sigma^2)T+\sigma\sqrt{T}Z_i)\), then apply the payoff to each simulated terminal value and discount it. Because the terminal distribution is exact, this procedure has sampling error but no time-discretization bias from an Euler or Milstein time grid.

Black–Scholes Prices from the Payoff Expectations

Set \(d_1=\dfrac{\log(S_0/K)+(r+\tfrac12\sigma^2)T}{\sigma\sqrt{T}}\) and \(d_2=d_1-\sigma\sqrt{T}\). Write \(\Phi(x)\) for the standard normal cumulative distribution function. The formulas for the call and put follow by splitting each expectation at the exercise threshold \(S_T=K\).

Proof. Represent \(S_T\) using \(Z\sim N(0,1)\). The event \(S_T>K\) is equivalent to \(Z>-d_2\), because taking logarithms gives \(\sigma\sqrt{T}Z>\log(K/S_0)-(r-\tfrac12\sigma^2)T=-\sigma\sqrt{T}d_2\). Thus \(\mathbb{Q}(S_T>K)=\Phi(d_2)\). For the call, its price is \(e^{-rT}\mathbb{E}^{\mathbb{Q}}[S_T\mathbf{1}_{\{S_T>K\}}]-Ke^{-rT}\mathbb{Q}(S_T>K)\). In the first term, the normal density \(\phi(z)\) satisfies \(e^{\sigma\sqrt{T}z}\phi(z)=e^{\sigma^2T/2}\phi(z-\sigma\sqrt{T})\), by completing the square in the exponent. Changing variable to \(y=z-\sigma\sqrt{T}\), the condition \(z>-d_2\) becomes \(y>-d_1\). Consequently, \(\mathbb{E}^{\mathbb{Q}}[S_T\mathbf{1}_{\{S_T>K\}}]=S_0e^{rT}\Phi(d_1)\). Substitution gives the call formula. For the put, \(S_T<K\) is equivalent to \(Z<-d_2\), and the same density identity gives \(\mathbb{E}^{\mathbb{Q}}[S_T\mathbf{1}_{\{S_T<K\}}]=S_0e^{rT}\Phi(-d_1)\). Writing the put expectation as \(Ke^{-rT}\mathbb{Q}(S_T<K)-e^{-rT}\mathbb{E}^{\mathbb{Q}}[S_T\mathbf{1}_{\{S_T<K\}}]\) proves the put formula. The threshold has probability zero under the continuous normal distribution, so the strict or non-strict inequality at \(S_T=K\) does not change either expectation. \(\square\)

Monte Carlo Estimators and Their Sampling Error

Let \(Y_i=e^{-rT}(S_T^{(i)}-K)^+\) be the discounted call payoff from the \(i\)-th independent simulation. The estimator is \(\widehat C_N=N^{-1}\sum_{i=1}^N Y_i\). For a put, replace the payoff by \(e^{-rT}(K-S_T^{(i)})^+\). These are independent, identically distributed observations. They have finite variance: a call payoff is bounded above by \(S_T\), a put payoff by \(K\), and geometric Brownian motion has a finite second moment at finite \(T\).

Proof. Linearity of expectation gives \(\mathbb{E}^{\mathbb{Q}}[\widehat C_N]=N^{-1}\sum_i\mathbb{E}^{\mathbb{Q}}[Y_i]=C_0\). Independence gives \(\operatorname{Var}(\widehat C_N)=N^{-2}\sum_i\operatorname{Var}(Y_i)=v/N\). To verify the sample-variance statement, use \(\sum_i(Y_i-\widehat C_N)^2=\sum_i(Y_i-C_0)^2-N(\widehat C_N-C_0)^2\), obtained by expanding both sides and using \(\sum_i(Y_i-\widehat C_N)=0\). Taking expectations, the first sum has expectation \(Nv\), while the second term has expectation \(N(v/N)=v\). Hence the expected numerator is \((N-1)v\), and division by \(N-1\) gives \(\mathbb{E}^{\mathbb{Q}}[s_N^2]=v\). The variance formula yields standard deviation \(\sqrt{v/N}\); substituting the estimated variance \(s_N^2\) gives the estimated standard error. \(\square\)

The estimator's standard error measures random fluctuation across repeated sets of \(N\) simulations. For large \(N\), an approximate confidence interval can be formed as \(\widehat C_N\pm1.96s_N/\sqrt{N}\), using the standard central limit theorem for independent identically distributed variables with finite variance. This is an approximation, not an exact finite-sample guarantee.

A Common Mistake and Why Exact Simulation Matters

Apply the payoff to each path before averaging. In the four-path example, the average terminal value is \((e^{-2}+1+e^2+e^{-4})/4=(0.135335+1+7.389056+0.018316)/4=2.135677\). An incorrect calculation would apply the call payoff to this average and report \((2.135677-1)^+=1.135677\). The correct estimate is the average of the four individual payoffs, \(1.597264\). These differ because \(x\mapsto(x-K)^+\) is nonlinear. A payoff is evaluated path by path at maturity.

For geometric Brownian motion, direct terminal simulation also avoids a separate numerical error source. An Euler simulation over a time grid approximates the path and can introduce time-discretization bias; the exact terminal formula does not. Thus, when using the exact formula, the comparison between the Monte Carlo estimate and the Black–Scholes price isolates sampling error (along with ordinary numerical rounding). Increasing \(N\) reduces the standard error at the rate \(N^{-1/2}\), while it does not change any time-discretization error because there is no time discretization in this method.

Calls and puts can be estimated by the same procedure, with only the payoff changed. The distinction between the pricing formula and its Monte Carlo estimate is important: the Black–Scholes expression is the exact expectation under the model, whereas a finite simulation produces a random sample average. A sound report gives the estimate together with its sampling uncertainty and identifies whether terminal values were simulated exactly or approximated on a time grid.