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

[Paul Ringseth]

Hi:
    I think this could work.  

    Now, I always find it useful to identify the hard parts up front.  Let's see if I can sketch this in a couple minutes.

1)  You need to estimate variance (vol^2) from the historical prices.  Now 50 days closing prices is a little low and when close to expiration you'd want some intraday data too.  But the model is a little crude anyhow, so perhaps it doesn't matter.  So assume the stock price St follows a log-normal process:  dlog(St) = (r - vol^2/2) dt + vol dBt and take first differences of log(St) for the returns process Rt ~ log(St/S(t-1)).  Rt is normally distributed so vol^2 ~ Sum{(Rt - R(t-1))^2}/48 - (1/49 Sum{Rt})^2}.  
2)  You need the risk-free rate, so set r=treasusry-rate for the same maturity as the stock option time-to-expiration -- this is in calendar days.
3)  Let t = today and T=expiration. 
    The sample paths of the stock price can be assumed to follow: Ss -> exp{r s + vol N(0, 1) sqrt(s)}  // N(0,1) =normal distribution.  
                 Let c_(ns) = exp{-r (T-t)} max{ exp{r (T-t) + vol ns sqrt(T-t)} - X, 0} -- where ns is a random sample from N(0,1).
                 european call(X,S,r,vol,T) = { c_(ns1) + c_(ns2) + ... + c_(ns(num_samples)) } / num_samples.
   
My math could be off a little (e.g. maybe a ito factor), but at least for me, it really helps to get some brief clean calculations up front, and then code.;-)


________________________________________
From: r-sig-finance-bounces at r-project.org [r-sig-finance-bounces at r-project.org] on behalf of msalese [massimo.salese at gmail.com]
Sent: Tuesday, November 22, 2011 7:30 AM
To: r-sig-finance at r-project.org
Subject: Re: [R-SIG-Finance] Option valuation for arbitrary distribution using monte carlo simulation

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)




--
View this message in context: http://r.789695.n4.nabble.com/Option-valuation-for-arbitrary-distribution-using-monte-carlo-simulation-tp4095718p4096122.html
Sent from the Rmetrics mailing list archive at Nabble.com.

_______________________________________________
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.