Option valuation for arbitrary distribution using monte carlo simulation
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.
On 23/12/2011 14:59, Joachim Breit wrote:
I gave it a try. See the (ugly) code below. Comments are welcome. There
are still a few issues to solve:
- How can one do the Monte Carlo simulation faster/more efficient?
- Can somebody explain Winston's daily risk free rate?
- Is there a better optimization algorithm for the utility function than
genoud?
----------------------------------------------
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
ov=0 #### monte carlo simulation
for (i in 1:100000) {
ov[i] = max(80 - 90 * prod(1+ sample(ret, size=63, replace = TRUE, prob
= rnprob)), 0 )
}
mov = 1/1.08^(63/252) * mean(ov) ### 2.512722 = option value
Am 25.11.2011 09:48, schrieb Charles Ward:
Luenberger discusses the task of creating risk-neutral probabilities from any set of observed returns. I cannot find the specific reference in his book (Investment Science) but there is a discussion paper by Winston (google) that applies his method using Excel and @Risk to a sample of Microsoft returns. I think it would be relatively simple to program it in R. Charles Ward
_______________________________________________ R-SIG-Finance at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-finance -- Subscriber-posting only. If you want to post, subscribe first. -- Also note that this is not the r-help list where general R questions should go.
Patrick Burns patrick at burns-stat.com http://www.burns-stat.com http://www.portfolioprobe.com/blog twitter: @portfolioprobe