I try to valuate option by using GARCH simulation, but I have call-put
disparity of estimated option premiums.

 

Maybe the reason is non-usage of risk neutral measure, but I have option on
futures contract with interest rate = 0, and short time interval = 24/250.

 

The R code:

 

# GARCH modeling to option valuation

# oleg.mubarakshin@gmail.com

 

require(rugarch)

require(fOptions)

 

# donwload file with market data

file_name_http = 'http://quant-lab.com/wp-content/uploads/2013/04/S_ret.txt'

file_name_here = 'prices_returns.txt'

download.file(file_name_http, destfile = file_name_here, method =
'internal', mode = "wb")

 

market_data <- read.table(file_name_here, sep = ' ', header = TRUE,
stringsAsFactors=FALSE)

tail(market_data)

 

S = market_data$S

ret = market_data$ret

rm(market_data)

 

# uGARCHspec

spec = ugarchspec(variance.model=list(model="apARCH", garchOrder=c(1,1)),

               mean.model=list(armaOrder=c(0,0), include.mean=FALSE),

               distribution.model="sstd", fixed.pars=list(omega=0))

print(spec)

 

# uGARCHfit

fit = ugarchfit(data=ret, spec=spec, solver="solnp",

               fit.control = list(stationarity=1, fixed.se=1, scale=0,
rec.init="all"))              

print(fit)

 

# uGARCHsim

nSim = 24 # the simulation horizon.

mSim = 100000 # the number of simulations

sim = ugarchsim(fit, n.sim=nSim, n.start=1, m.sim=mSim,
startMethod="sample")

 

# Density of S(T) distribution

St = NULL

for(i in 1:mSim)

               St[i] = S[length(S)] * exp(sum(sim@simulation$seriesSim[,i]))

d = density(St)

min_x = min(d$x); max_x = max(d$x)

f_density = splinefun(d$x, d$y)

err = as.numeric(integrate(f_density, min(d$x), max(d$x), subdivisions=1000,
abs.tol = 1e-8)$value)

 

# Option pricing

f <- function(x) {

               ret_val = NULL

               for (i in 1:length(x)) {

                              if (x[i] >= min_x & x[i] <= max_x) {

                                            ret_val[i] = f_density(x[i])

                              } else {

                                            ret_val[i] = 0

                              }

               }

               return(ret_val)

}

               

d_payoff <- function(x) {

               ret_val = NULL

               if (put) {

                              for (i in 1:length(x)) {

                                            if (strike > x[i]) {

                                                           ret_val[i] =
(strike - x[i]) * f(x[i])

                                            } else {

                                                           ret_val[i] = 0

                                            }

                              }

               } else {

                              for (i in 1:length(x)) {

                                            if (strike < x[i]) {

                                                           ret_val[i] =
(x[i] - strike) * f(x[i])

                                            } else {

                                                           ret_val[i] = 0

                                            }

                              }

               }

               return(ret_val)

}

 

dK = 5000; kK = c(0.7, 1.3)

K = NULL; i = 1; K[i] = round(S[length(S)] * kK[1]/dK, 0) * dK

while (K[i] < S[length(S)] * kK[2]) {

               i = i + 1

               K[i] = K[i - 1] + dK             

}

               

C = P = NULL

for (i in 1:length(K)) {

               strike = K[i]

               put = FALSE; C[i] = as.numeric(integrate(d_payoff, K[i],
max_x, subdivisions=1000, abs.tol = 1e-8)$value)/err

               put = TRUE; P[i] = as.numeric(integrate(d_payoff, min_x,
K[i], subdivisions=1000, abs.tol = 1e-8)$value)/err

}

 

# Put-Call Parity prices

Cp = P - K + S[length(S)]

Pc = C + K - S[length(S)]

 

dev.new()

plot(K, C, type='l', col='blue', ylim=range(C,P,Cp,Pc), xlab='Strike',
ylab='Level', main='Option prices', 

               sub='Call, Put - blue lines; Sintetic Call, Put - red lines')

lines(K, Cp, type='l', col='red')

lines(K, P, type='l', col='blue')

lines(K, Pc, type='l', col='red')

grid(col='lightgray', lty='dotted')

 

dev.new()

plot(K, C - Cp, type='l', col='blue', ylim=range(C-Cp,P-Pc), xlab='Strike',
ylab='Level', main='Option prices', 

               sub='Call - Sintetic Call (blue line), Put - Sintetic (red
line)')

lines(K, P - Pc, type='l', col='red')

grid(col='lightgray', lty='dotted')

 

# Implied Volatility

T = nSim/250; tol = 1e-6; IV = IVc = IVp = NULL

for(i in 1:length(K)) {

               if(K[i] > S[length(S)]) {

                              if(C[i] > tol) {

                                            IV[i] =
max(GBSVolatility(price=C[i], TypeFlag='c', S=S[length(S)], 

                                                           X=K[i], Time=T,
r=0, b=0, tol=1e-8, maxiter=1e8), 0)

                              } else {

                                            IV[i] = NA

                              }

               } else {

                              if(P[i] > tol) {

                                            IV[i] =
max(GBSVolatility(price=P[i], TypeFlag='p', S=S[length(S)], 

                                                           X=K[i], Time=T,
r=0, b=0, tol=1e-8, maxiter=1e8), 0)

                              } else {

                                            IV[i] = NA

                              }

               }

               if(C[i] > tol) {

                              IVc[i] = max(GBSVolatility(price=C[i],
TypeFlag='c', S=S[length(S)], 

                                                           X=K[i], Time=T,
r=0, b=0, tol=1e-8, maxiter=1e8), 0)

               } else {

                              IVc[i] = NA

               }

               if(P[i] > tol) {

                              IVp[i] = max(GBSVolatility(price=P[i],
TypeFlag='p', S=S[length(S)], 

                                                           X=K[i], Time=T,
r=0, b=0, tol=1e-8, maxiter=1e8), 0)

               } else {

                              IVp[i] = NA

               }

}

 

dev.new()

plot(K, IV, type='l', col='black', ylim=range(IV,IVc,IVp), xlab='Strike',
ylab='IV', main='Implied Volatility', 

               sub='IV by OTM (black), IV by Call (red), IV by Put (blue)')

lines(K, IVc, type='l', col='red')

lines(K, IVp, type='l', col='blue')

grid(col='lightgray', lty='dotted')

 

 

 


Sincerely,

Oleg Mubarakshin

 


	[[alternative HTML version deleted]]

