[R-SIG-Finance] Option valuation for arbitrary distribution using monte carlo simulation

[Joachim Breit]

Thank you, Patrick. That was indeed a huge speed-up. I can now easily do 
the monte carlo simulation with 2e6 trials and the calculated option 
value seems quite stable:

1e5 trials --> option values 2.481841, 2.464703, 2.505218
1e6 trials --> 2.491286, 2.502324, 2.491899, 2.495442
2e6 trials --> 2.493493, 2.495717

--------------------------------
library("fImport")
library('fOptions')
library("rgenoud")
msft = 
yahooImport(query="MSFT&a=00&b=4&c=1999&d=03&e=9&f=1999&ignore=.csv",file = 
"msft", source = "http://ichart.finance.yahoo.com/table.csv?", save = TRUE)

cl = msft at data[,"Close"]   ### unadjusted close prices
cl[10:67] = cl[10:67]/2    ### a 2:1 split on 1999-03-29

ret = cl[1:66] /cl[2:67] - 1   ### the returns
logret = log(1+ret)            ### log returns
vola = sd(logret) * sqrt(252)    ###  0.4091751
bsm = GBSOption(TypeFlag = "p", S = 90, X = 80, Time = 63/252, r = 
log(1+0.08), b = log(1+0.08), sigma = vola)@price  #### 2.577398 (just 
to compare)

rf = .000317   ## don't know where this comes from in Winston's paper;
                ## should be 1.08^(1/252)-1 = 0.0003054476, right?

meanutil <- function (alpha) {
x = log(alpha * (1 + ret) + (1 - alpha) * (1 + rf))  ### portf. utility
x[which(is.nan(x))] = -100000    ### correct for ln(x<0) undefined
x = mean(x)
return(x)
}

optalpha = genoud(meanutil, nvars = 1, max = TRUE, pop.size = 1000)$par
### the optimal portfolio under utility function ln(x)

portret = optalpha * (1 + ret) + (1 - optalpha) * (1 + rf)
## the optimal portfolio returns (in the sense of: daily end wealth)

num = (1 / length(ret)) * (1 / portret)
rnprob = num/sum(num)  ### the risk neutral probabilities

sampmat <- array(sample(logret,    size=63 * 2e6, replace=TRUE, 
prob=rnprob), c(63, 2e6))   ### array with 63 rows ans 2000000 columns
samplret <- colSums(sampmat)   #### calc sum for each column
ov = pmax(0, 80 - 90 * exp(samplret))
mov = 1/1.08^(63/252) * mean(ov) ### gecheckt, dass da dasselbe 
rauskommt wie bei der Schleife!!!
rm(sampmat, samplret, ov)



Am 23.12.2011 16:35, schrieb Patrick Burns:
> Joachim,
>
> I'm guessing that you could get
> some better speed by doing the
> 'sample' call all at once, do it
> with 'logret' instead of 'ret'
> and sum them efficiently. Something
> like:
>
> sampmat <- array(sample(logret,
> size=63 * 1e5, replace=TRUE,
> prob=rnprob), c(63, 1e5))
>
> samplret <- colSums(sampmat)
>
> # then transform to simple returns
> # and use 'pmax' instead of 'max'
>
> If the matrix is too large for your
> memory, you could do it in a few
> blocks.
>