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

[Joachim Breit]

Dear Massimo,

thanks for the hint to Iacus' book. I was aware of this title already 
but could not get hold of it by now. Thank you for the code snippet as 
well. Unfortunately, it does not seem to help with my main problem: How 
to correct the original data such that Monte Carlo sampling from it 
gives me samples from a risk-neutral but NONNORMAL process. Maybe my 
description is a bit awkward. Also I may lack the math and statistics 
skills. I try do say it another way:

Let
r : original log return series
     ( with e.g. mean(r)  = -0.0001574849 and sd(r)  = 0.01479017 )
x:  risk-free rate of, say, 5%.
Snow : stock price now
time to maturity = 1 year = 252 trading days
St : stock price in 1 year
X: call option strike
E[] : expected value
C : the option price

I am looking for a corrected return series rcorr = r + z such that
E[ Snow * exp(  sum(  sample(rcorr, 252)     ) ) ] = 100 * exp(0.05)

because then I can find the call option price with
for (i in 1:bignumber) {
   St[i] = Snow * exp(  sum(   sample(rcorr, 252)     ) )
   op[i] = max(0, exp(-.05) * (St[i] - X)  )
}
C = mean(op)

I could of course use some search algorithm to find z. In fact thats 
what I did. But I thought this is so common a problem for finance folks 
that there must be some elegant R procedure for th task.

Your code seems to assume that we are dealing with normal log returns 
but this is not the case for my sample.

Do you see what I mean? Is this stuff dealt with in the Iacus book?

Greetings
Joachim

>
> Hi Joachim , I'm not a prof but I recommend you the following book:
> Option Pricing and Estimation of Financial Models with R
> Stafano M. Iacus
> http://www.amazon.com/Option-Pricing-Estimation-Financial-Models/dp/0470745843/ref=sr_1_48?ie=UTF8&qid=1321974343&sr=8-48
>
> it'll explain a lot of thinks with code samples.
> Some very base simple code, I started with this one and while was reading
> the book I started to improve it
>
> #define functions
> require(package='fBasics',quietly=FALSE)
>
> #This function to create innovation
> #We keep innovation function autside of process because
> #once computed leave it stored for next usage
> innovation.norm<-function(pTrials=5000){
>    # remember to add the seed
>    #generate standard normal innovation
>    z<- rnorm(n=pTrials,mean=0,sd=1)
>    #genereta anthitetic innovation
>    ant.z<--z
>    #combine in the new one:
>    z<- c(z,ant.z)
>    #return innovation vector
>    return(z)
> }
>
> #this function to create the price distribution
> st.proc.norm<-function(pInnVect=z,mu=0,pSnow=15,pVol=0.15,pTradDaysToExp=20){
>    #set the number of trading days in one year
>    yearTradDays<-252
>    #compute the interval time
>    dt<-pTradDaysToExp/yearTradDays
>    s.final<- rep(0,length(pInnVect))
>    s.final<-pSnow*exp(pVol*sqrt(dt)*pInnVect)
>    out<-c(s.final,pSnow)
>    return(s.final)
> }
>
> #this function to build the option chain
> price.chain<-function(pSFinal,pIntRate=0.02,pStrikeMin=10,pStrikeMax=15,pDStrike=0.5,pCalDaysToExp=30){
>    #build strike chain
>    K<-seq(from=pStrikeMin,to=pStrikeMax,by=pDStrike)
>    #compute dicount factor
>    calDays<-365
>    dt<-pCalDaysToExp/calDays
>    #prepare storage for chain
>    prCall<-rep(x=0,length(K))
>    prPut<-rep(x=0,length(K))
>    se.call<-rep(x=0,length(K))
>    se.put<-rep(x=0,length(K))
>
>    for (i in (1:length(K))){
>      #print(i)
>      #print(K[i])
>      #formula for pricing
>      c<- ifelse(pSFinal>K[i],exp(-pIntRate*dt)*(pSFinal-K[i]), 0)
>      p<- ifelse(pSFinal<K[i],exp(-pIntRate*dt)*(K[i]-pSFinal), 0)
>      #price as average
>      prCall[i]<-mean(c)
>      prPut[i]<-mean(p)
>      #standard error
>      se.call[i]<-sd(c)/sqrt(length(pSFinal))
>      se.put[i]<-sd(p)/sqrt(length(pSFinal))
>
>      ## make a confidence interval estimate
>      alpha<-0.05
>      t.star<-qt(p=(1-alpha/2),df=length(pSFinal)-1,lower.tail=TRUE)
>      prCall[i]+c(-1,1)*t.star*se.call[i]
>      prPut[i]+c(-1,1)*t.star*se.put[i]
>    }
>
>    #build the data frame for option's chain
>
> pr.chain<-data.frame(CALL=round(prCall,digits=3),Strike=K,PUT=round(prPut,digits=3))
>
>    return(pr.chain)
> }
>
>
> ##################################################################################
> ##  PRICING WITH NORMAL
> inn.n<-innovation.norm(pTrials=500)
> hist(inn.n)
> #you must estimate volatility from past data
> s.final.n<-st.proc.norm(pInnVect=inn.n,pSnow=13.40,pVol=0.25,pTradDaysToExp=20)
> hist(s.final.n)
> #you must estimate interest rate
> price.chain(pSFinal=s.final.n,pIntRate=0.02,pStrikeMin=10,pStrikeMax=15,pDStrike=0.5,pCalDaysToExp=20)