Feeds:
Posts

## Option pricing with long-memory stochastic volatility

Christodoulou and Viens have produced an algorithm for option pricing with long memory stochastic volatility using particle filters and their empirical distributions.  Their paper is ChronopoulosViensEstimationPricingLMSV.  I should have a usable interpretation of their algorithm in the next days in R.

Since I am fond of simple algorithms, here is one that I will be playing with.  Suppose we have returns $x_1, \dots x_N$ where $N$ is today and then we want to price a call option 30 days in the future.  First we estimate the parameters of an ARFIMA(1,d,1) using $y_t = log(x_t^2+c)$ where $c=0.000001$ and then we simulate $n$ series using this model.  Then we set $\phi(x) = e^{-2|x|}$ and the contraction $\phi_N = N^{1/3}\phi(N^{1/3}x)$ to calculate a probability distribution using the $N$-th component of the predicted $x$ from the simulated $y_N^j$ with $j=1,\dots,n$.  The probability distribution will be normalization of $\phi_N(X_N^j - x_N)$ If we simulated series of length $N+30$ then this probability distribution can weigh the endpoint distribution when we calculate $E( x_{N+30} - S_K )_+$

We can do the implementation in two steps.  First, a function that calculates the call price given the underlying volatility; then the function that simulates $MC$ arfima processes and then provides probabilities for each based on the observation at the last known point.

> callPrice
function(svol,P0,Strike){
P0<-log(P0)
Strike<-log(Strike)
n<-length(svol)
callP<-0
for (k in 1:100){
z<-rnorm(n,sd=1.0)
e<-sum( exp(svol/2) %*% z )
p<-0
if (e+P0-Strike>0){
p<-e+P0-Strike
}
callP<-callP+(1/100)*p
}
callP
}

> callPriceLM
function(rets, T, P0, Strike){
z<-as.numeric(log(rets^2+0.000001))
z[is.na(z)]<-log(0.000001)
fit<-arfima(z,order=c(1,0,1))
phi<-fit$modes[[1]]$phi
theta<-fit$modes[[1]]$theta
sigma2<-fit$modes[[1]]$sigma2
df<-fit$modes[[1]]$dfrac
N<-length(rets)
N1<-N+T
MC<-1000
initY<-rep(0,MC)
probY<-rep(0,MC)
callP<-rep(0,MC)
for (k in 1:MC){
qq<-arfima.sim(N1,model=list(phi=phi,theta=theta,frac=df))
initY[k]<-qq[N]
h<-rnorm(1)
x<-h*exp(initY[k]/2)
probY[k]<-exp(-2*abs(x-rets[N]))
callP[k]<-callPrice(qq[N:N1],P0,Strike)
}
for (k in 1:MC){
probY[k]<-probY[k]/sum(probY)
}
p<-0
p<-sum(probY %*% callP)
p
}
>