[R-SIG-Finance] Monte Carlo Option Pricing formula

[rex]

Roupell, Darko <Darko.Roupell at cba.com.au> [2012-02-01 22:45]:
>
>I am trying to cross check option implied employee option price that
>was derived using Monte Carlo simulation. Below is code and parameter
>used and after accounting for dividend yield ( 1.46%) the derived
>option price is 206.8843 using the code snippet provided. Approx 1
>cent below 207.95 that is listed in company prospect using their
>Monte Carlo simulation and below parameters.
>
>As we all know number of iteration can also slightly impact the
>average price, but am I rightly concerned that my methodology may not
>be correct?

Hi Darko,

Yes -- or at least different from the company's assumptions. Enrico
has informally shown that there is a statistically significant
difference in the two means. Conceptually, it may be important, but
practically, the difference is insignificant.

I have 3 other observations: first, your code uses a "for" loop which
should be avoided in code that takes a significant amount of time to
run in R, if possible. rnorm() can obviously be moved out of the loop
and called once to generate the eps vector, but delS[i] = max(mv,0) can
also be moved out of the loop, and then there's no need for a loop.

To see the difference eliminating the loop makes, let's compare your
code to a no-loop version:

start = Sys.time()
exercisePrice   = 0;
timeToExpiry    = 3;        #% in years
underlyingPrice = 490;      #% underlying in cents
expectedVol     = 0.5;      #% expected volatility
expectedDiv     = 0.0146;   #% expected dividend in cents
riskFreeRate    = 0.0425;   #% interest rate
itr = 500000                #% number of iterations
delS = 0*array(0,itr)

for( i in (1:itr))
{
      eps = rnorm(1)          #% random number generator
      dS = expectedDiv*underlyingPrice+underlyingPrice*(riskFreeRate*timeToExpiry) + (underlyingPrice*expectedVol*eps*sqrt(timeToExpiry))
      mv = dS - exercisePrice;

      delS[i] = max(mv,0);
}
mean(delS)
#206.3257
Sys.time - start
Time difference of 36.64036 secs

BTW, you multiply the dividend by the stock price, which is
inconsistent with your"expected dividend in cents" code comment. Plus,
it's a percentage per year, right? If so, I think the first term in
your dS expression should be: expectedDiv*underlyingPrice*timeToExpiry.

I'm going to rewrite your code using short variable names; long 
names make the math itself more difficult for me to immediately see.

start <- Sys.time()
E   <- 0             #exercise price
t   <- 3             #option life (years)
S   <- 490           #stock price (cents)
v   <- 0.5           #est. future annual volatility
d   <- 0.0146        #dividend (code suggests this is a percentage of S)
r   <- 0.0425        #annual interest rate
N   <- 5000000
z   <- S*(d+r*t + rnorm(N, 0, v*sqrt(t)))
sum(z[z > E])/N
#206.203
Sys.time() - start
Time difference of 0.5397904 secs

The loop code takes 68 times as long to run.

Second, why are you using a normal distribution? The usual rough
approximation is to assume a lognormal distribution for stock prices,
but that's flawed, too.

Third, what is the annual interest rate? From your code, it appears to
be the drift rate of the process driving the stock price, which has
nothing to do with interest rates. A more realistic model will
estimate the future drift rate, and using that plus the other
parameters, will calculate a future value for the option. The present
value is that future value discounted at the rate you choose (which is
unlikely to be the riskless rate). The Black-Scholes-Merton derivation
has confused many people on this point.

HTH,

-rex
-- 
For every complex problem there is an answer that is clear, simple, 
and wrong. --H. L. Mencken