Simulation Estimation of Option Prices

Philip H. Dybvig

Washington University

Saint Louis, Missouri

* 1 V = E [- V ] 0 r 1This assumes no dividends. E* indicates expectations under the risk-neutral probabilities that equalize expected returns across assets, and r is one plus the riskless rate. Intuitively, investing in any asset is a fair gamble after discounting (by 1/r) adjusts for time preference and changing from actual probabilities to risk-neutral probabilities adjusts for risk preferences.

* 1 V = E [-- V ] 0 0 r0 1 * 1 * 1 = E [-- E [-- V ]] 0 r0 1 r1 2 * 1 = E [----- V ] 0 r0 r1 2 * 1 = ... = E [-------------- V ] 0 r0 r1 ... rT-1 T

- If the the interest rate is nonrandom, we can take the discounting out of the expectation, and we simply discount the expected cash flows in the risk-neutral probabilities. In a simulation, we simulate the underlying processes using the risk-neutral probability distributions, compute the cash flows, and then take net present values.
- If the interest rate moves randomly, we need to discount
using the realized spot rates, and we take the expectation
*after*discounting.*Do not*discount discount expected terminal cash flows at a multiperiod riskless rate. - For multiple cash flows, we can treat each maturity separately and add the results. Or, we can accumulate cash flows and reinvest them in the riskless asset which will give the same answer.

- x1 = a + (b-a) z1: uniform on the interval [a,b]
- x2 = (z1 < p ? u : d), where 0 < p < 1: Bernoulli random variable that is u with probability p and d with probability 1-p
- x3 = (z1 < p1 ? u : (z1 < p1 + p2 ? m : d)), where 0 < p1 < p1 + p2 < 1: three-valued random variable giving u with probability p1, m with probability p2, and d with probability p3 = 1 - p1 - p2.
- x4 = sqrt(-2.0 * log(z1))*cos(2*pi*z2) and x5 = sqrt(-2.0 * log(z1))*sin(2*pi*z2): two independent normal random variables each having mean zero and variance 1. (The constant pi ~ 3.1415926535.)
- x6 = z1 + z2 + z3 + z4 + z5 + z6 - z7 - z8 - z9 - z10 - z11 - z12: a different variable approximately normal with mean zero and variance 1.
- x7 = -k * log(z1): exponential distribution with mean k
- x8 = m + s * x4: a normal random variable with mean m and standard deviation s and therefore variance s^2. (Obviously, we can use x5 or x6 in place of x4.)
- x9 = m + s * (z1-0.5) + sqrt((double) 12): uniform random variable with mean of m and standard deviation s.

// // simu2.h Stochastic Volatility Option Pricing Model // class svprice { public: svprice(double ttm=0.25,int nper=12,double r=.05,double initstd=.15, double meanstd=.2,double k=6.0,double sigstd=.5); double eurcall(double stock,double strike,long int nsimu=1000); private: int npers; double tinc, r1per, stockP, sigma0, sigma, meansigma, sigmasigma, kappa, c0, c1, c2, c3, c4, c5; double stocktotret();};

// // simu2test.cc Stochastic Volatility Option Pricing Test File // #include <iostream.h> #include <math.h> #include "simu2.h" main() { long int i,nsimus; svprice sim2; for(nsimus=10;nsimus<=100000l;nsimus *= 10) cout << nsimus << " " << floor(sim2.eurcall((double) 100.0,(double) 100.0,nsimus)*100.0+0.5)/100.0 << endl; cout << endl; for(i=0;i<20;i++) cout << "100000 " << floor(sim2.eurcall((double) 100.0,(double) 100.0,(long int) 100000l)*100.0+0.5)/100.0 << endl; return(0);}

// // simu2.cc Stochastic Volatility Option Pricing Model // #include <math.h> #include <stdlib.h> #include <iostream.h> #include "simu2.h" #define unifrand() ((double) rand()/((double) RAND_MAX)) #define MAX(a,b) (((a) > (b)) ? (a) : (b))

svprice::svprice(double ttm,int nper,double r,double initstd, double meanstd, double k, double sigstd) { npers = nper; tinc = ttm/(double) nper; r1per = 1.0 + r*tinc; sigma0 = initstd; meansigma = meanstd; sigmasigma = sigstd; kappa = k; c0 = kappa * tinc * meansigma; c1 = 1.0 - kappa * tinc; c2 = 1.0 - sigmasigma * sqrt(tinc)*0.5*sqrt((double) 12); c3 = sigmasigma * sqrt(tinc) * sqrt((double) 12); c4 = sqrt(tinc)*sqrt((double) 12);}

double svprice::eurcall(double stock,double strike,long int nsimu) { long int i,n; double x; x=0.0; for(n=0;n<nsimu;n++) { stockP = stock; sigma = sigma0; for(i=0;i<npers;i++) { stockP *= stocktotret(); } x += MAX(stockP-strike,0);} return(x/(double) nsimu);} double svprice::stocktotret() { double stockret; // // The following straightforward computations are algebraically the same as // the ones used below but are much slower because many more calculations are // performed in each pass through the loop. // // stockret = r1per + sigma * sqrt(tinc) * (unifrand()-0.5) * sqrt((double) 12); // sigma = (kappa*tinc * meansigma + (1.0 - kappa * tinc) * sigma) * // (1 + sigmasigma * sqrt(tinc) * (unifrand()-0.5) * sqrt((double) 12)); // return(stockret);} // stockret = r1per + sigma * c4 * (unifrand()-0.5); sigma = (c0 + c1 * sigma) * (c2 + c3 * unifrand()); return(stockret);}