Monte Carlo
Monte Carlo for Pricing Option
The first large-scale applications of the Monte Carlo method were in physics and arose from work on the Manhattan Project in the 1940s. Boyle (1976), applied the method innovatively to Finance and triggered decades of further research. The method uses the fact that the distribution of terminal stock prices is determined by the process generating future stock price movements. This process can be simulated on a computer thus generating a series of stock price trajectories. This series determines a set of terminal stock values which can be used to obtain an estimate of the option value. Furthermore, the standard deviation of the estimate can be obtained at the same time so that the accuracy of the results can be established. Monte Carlo is based on simulating Stock Price behaviour. To get some graphical intuition into this process please consider this playlist as a pre-amble to understanding the code below.
VBA code for Black Scholes Monte Carlo
Function MCCall(S, X, T, r, q, v, n)
' Espen Haug with small tweaks
' D is the drift
' r is the risk free rate
' q is the dividend yield
' T is the Maturity
' v is the volatility of the underlying
' asset return or Black Scholes sigma
' S the underlying asset price
' ST is the Terminal asset price at expiry
' X is the exercise price
D = (r - q - v ^ 2 / 2) * T
vT = v * Sqr(T)
For i = 1 To n
ST = S * Exp(D + vT * Application.NormSInv(Rnd()))
bin = bin + Application.Max((ST - X), 0)
Next
MCCall = Exp(-r * T) * bin / n
End Function
Monte Carlo Simulation of Stock Prices see playlist below
Below we build up incrementally the basics governing stock price simulation which applies Geometric Brownian motion. This approach means we can estimate the value of an option using simulation. The Black Scholes (1973) model can be closely approximated using Monte Carlo. Accuracy improves as we in crease the number of pathways or simulations. Below, we explain why
C++ code for Black Scholes Monte Carlo
//Code originally sourced from Fabrice Rouah
//I made some small changes see this video for changes I made
//https://www.youtube.com/watch?v=p3jaWwY8ayU&t=2s
//see this video for explanation on Monte carlo and GBM
//https://goo.gl/d9W9Gs
using namespace std;
#include <math.h>
#include <iostream>
#include <iomanip>
#include <vector>
#include <stdlib.h>
#include <time.h>
#include <algorithm>
// Calculates the mean of a vector of values
double
VecMean (vector < double >x)
{
int n = x.size ();
double sum = 0.0;
for (int i = 0; i <= n - 1; i++)
sum += x[i];
double xbar = sum / n;
return xbar;
}
// N(0,1) density
double
f (double x)
{
double pi = 4.0 * atan (1.0);
return exp (-x * x * 0.5) / sqrt (2 * pi);
}
// Boole's Rule
double
Boole (double StartPoint, double EndPoint, int n)
{
vector < double >X (n + 1, 0.0);
vector < double >Y (n + 1, 0.0);
double delta_x = (EndPoint - StartPoint) / double (n);
for (int i = 0; i <= n; i++)
{
X[i] = StartPoint + i * delta_x;
Y[i] = f (X[i]);
}
double sum = 0;
for (int t = 0; t <= (n - 1) / 4; t++)
{
int ind = 4 * t;
sum +=
(1 / 45.0) * (14 * Y[ind] + 64 * Y[ind + 1] + 24 * Y[ind + 2] +
64 * Y[ind + 3] + 14 * Y[ind + 4]) * delta_x;
}
return sum;
}
// N(0,1) cdf by Boole's Rule
double
N (double x)
{
return Boole (-10.0, x, 240);
}
// Black Scholes closed form price
double
BS (double S, double K, double v, double T, double r, double q, char PutCall)
{
double d1 = (log (S / K) + (r - q + 0.5 * v * v) * T) / v / sqrt (T);
double d2 = d1 - v * sqrt (T);
double Call = S * exp (-q * T) * N (d1) - K * exp (-r * T) * N (d2);
if (PutCall == 'C')
return Call;
else
return Call - S * exp (-q * T) + K * exp (-r * T);
}
int
main ()
{
srand (time (0)); // Set the seed for random number generation
double S = 100.0; // Spot Price
double K = 100.0; // Strike Price
double T = 1; // Maturity in Years
double r = 0.05; // Interest Rate
double q = 0; // Dividend yeild
double v = 0.2; // Volatility
int Nsims = 1e7; // Number of simulations
double Z; // Random standard normal Z(0,1)
vector < double >ST (Nsims, 0.0); // Initialize terminal prices S(T)
vector < double >ST_K (Nsims, 0.0); // Initialize call payoff
vector < double >K_ST (Nsims, 0.0); // Initialize put payoff
double u1, u2;
double pi = 3.141592653589793;
for (int i = 0; i <= Nsims - 1; i++)
{
// Independent uniform random variables
u1 = ((double) rand () / ((double) (RAND_MAX) + (double) (1)));
u2 = ((double) rand () / ((double) (RAND_MAX) + (double) (1)));
// Floor u1 to avoid errors with log function
u1 = max (u1, 1.0e-10);
// Z ~ N(0,1) by Box-Muller transformation
Z = sqrt (-2.0 * log (u1)) * sin (2 * pi * u2);
ST[i] = S * exp ((r - q - 0.5 * v * v) * T + v * sqrt (T) * Z); // Simulated terminal price S(T)
ST_K[i] = max (ST[i] - K, 0.0); // Call payoff
K_ST[i] = max (K - ST[i], 0.0); // Put payoff
}
// Simulated prices as discounted average of terminal prices
double BSCallSim = exp (-r * T) * VecMean (ST_K);
double BSPutSim = exp (-r * T) * VecMean (K_ST);
// Closed form prices
double BSCall = BS (S, K, v, T, r, q, 'C');
double BSPut = BS (S, K, v, T, r, q, 'P');
// Errors
double CallError = BSCall - BSCallSim;
double PutError = BSPut - BSPutSim;
// Output the results
cout << setprecision (4) << fixed;
cout << "Using " << Nsims << " simulations..." << endl;
cout << " " << endl;
cout << "Method CallPrice PutPrice " << endl;
cout << "----------------------------------" << endl;
cout << "Simulation " << BSCallSim << " " << BSPutSim << endl;
cout << "Closed Form " << BSCall << " " << BSPut << endl;
cout << "Error " << CallError << " " << PutError << endl;
cout << "----------------------------------" << endl;
cout << " " << endl;
system ("PAUSE");
}
Monte Carlo using Box Muller
The Box and Muller takes two samples from the uniform distribution on the interval [0, 1] and maps them to two standard, normally distributed samples - in effect mapping them into two normally distributed samples with the use of sine or cosine functions. A basic and hopefully intuitive explanation is provided in excel here:
The relevant segment of C++ code is highlighted below:
// Independent uniform random variables
u1 = ((double) rand () / ((double) (RAND_MAX) + (double) (1)));
u2 = ((double) rand () / ((double) (RAND_MAX) + (double) (1)));
// Floor u1 to avoid errors with log function
u1 = max (u1, 1.0e-10);
// Z ~ N(0,1) by Box-Muller transformation
Z = sqrt (-2.0 * log (u1)) * sin (2 * pi * u2);
Python Code for Black Scholes Monte Carlo
# http://www.codeandfinance.com/pricing-options-monte-carlo.html
import math
from random import gauss
from math import exp, sqrt
def generate_asset_price(S,v,r,q,T):
return S * exp((r - q - 0.5 * v**2) * T + v * sqrt(T) * gauss(0,1.0))
def call_payoff(S_T,K):
return max(0.0,S_T - K)
S = 100 # underlying price
v = 0.20 # volatility
r = 0.05 # risk free rate
q = 0.0 # dividend yield
T = 1 # maturity
K = 100. # exercise
simulations = 10000
payoffs = []
discount_factor = math.exp(-r * T)
for i in range(simulations):
S_T = generate_asset_price(S,v,r,q,T)
payoffs.append(
call_payoff(S_T, K)
)
call_price = discount_factor * (sum(payoffs) / float(simulations))
# From Put Call Parity
put_price = call_price + K * math.exp(-r * T) - S * math.exp(-q * T)
print (call_price, put_price)
Python Code for Black Scholes Monte Carlo and closed form solution
# http://www.codeandfinance.com/pricing-options-monte-carlo.html
import math
from random import gauss
from math import exp, sqrt, log, sqrt, pi
from scipy.stats import norm
def d1(S,K,T,r,q,sigma):
return(log(S/K)+(r - q +v**2/2.)*T)/v*sqrt(T)
def d2(S,K,T,r,q,v):
return d1(S,K,T,r,q,v)-v*sqrt(T)
## define the call options price function
def bs_call(S,K,T,r,q,v):
return S*exp(-q*T)*norm.cdf(d1(S,K,T,r,q,v))-K*exp(-r*T)*norm.cdf(d2(S,K,T,r,q,v))
## define the put options price function
def bs_put(S,K,T,r,q,v):
return K*exp(-r*T)-S*exp(-q*T)+bs_call(S,K,T,r,q,v)
def generate_asset_price(S,v,r,q,T):
return S * exp((r - q - 0.5 * v**2) * T + v * sqrt(T) * gauss(0,1.0))
def call_payoff(S_T,K):
return max(0.0,S_T - K)
S = 100 # underlying price
v = 0.20 # volatility
r = 0.05 # risk free rate
q = 0.0 # dividend yield
T = 1 # maturity
K = 100. # exercise
simulations = 1000000
payoffs = []
discount_factor = math.exp(-r * T)
for i in range(simulations):
S_T = generate_asset_price(S,v,r,q,T)
payoffs.append(
call_payoff(S_T, K)
)
call_price = discount_factor * (sum(payoffs) / float(simulations))
# From Put Call Parity
put_price = call_price + K * math.exp(-r * T) - S * math.exp(-q * T)
cbs = bs_call(S,K,T,r,q,v)
pbs = bs_put(S,K,T,r,q,v)
print(cbs, pbs)
print (call_price, put_price)