# Approximate asian geometric option with Heston

I am trying to implement Theorem 1 from this Journal in RStudio.

The journal says the it is possible to find a approximate price of a geometric asian option in a Heston setup this way: $$X_{1cGAO}=e^{0.5\hat{\sigma}-rT+a}N(d_+)-Ke^{-rT}N(d_-)$$ where $$a=log(S_0)+0.5rT- \frac{1}{2T}\int_0^T[\int_0^t\sigma^2(s)ds]dt$$ $$\hat{\sigma}^2=\frac{1}{n^2}\sum_{j=1}^n[2(n-j)+1]\int_0^{j\frac{T}{n}}\sigma^2(s)ds$$ where $n$ goes to $\infty$ and $d_-=\frac{a-log(K)}{\hat{\sigma}}$, $d_+=d_-+\hat{\sigma}$

The trick is, according to the journal to substitute the stochastic volatility from the Heston process with its deterministic expected value: $$E[\sigma^2_t]=v_0e^{-\kappa*t}+\theta(1-e^{-\kappa t})$$

And here comes my problem:

When i am trying to implement the code in RStuido (see code below) I do not get a plausibly result. I get the price $8.592864$ and when I am using Monte Carlo simulation I get about $5.430642$

Can someone see what I am doing wrong or what I am missing?

kappa=2
theta=0.04
v0=0.04
S=100
rf=0.05
K=100
T=1

MC_Heston_Std_Asian_GEO(S,K,rf,T,v0,theta,10^6,10,kappa,0,0)

#Expected value of CIR-Prices
ExpectedVarians <- function(t) {
exp(-kappa*t)*v0+(theta)*(1-exp(-kappa*t))
}

sqrt(ExpectedVarians(1))

#My function
AsianGEO = function(S,K,rf,T,n){

a=log(S)+0.5*rf*T-(1/(2*T))*integrate(function(y) {
sapply(y, function(y) {
integrate(ExpectedVarians, 0, y)$value }) }, 0, T)$value

sum=0
for (j in 1:n){
sum=sum+(2*(n-j)+1)*integrate(ExpectedVarians, 0, j*(T/n))\$value
}
sigma_hat=sqrt(sum/n^2)
d_minus=(a-log(K))/(sigma_hat)
d_plus=d_minus+(sigma_hat)

(exp(0.5*sigma_hat-rf*T+a)*pnorm(d_plus)-K*exp(-rf*T)*pnorm(d_minus))
}

AsianGEO(S,K,rf,T,1000)


MC_Heston_Std_Asian_GEO <- function(S,K,rf,T,V0,theta,n,b,kappa,sigma_vol,rho) {

values <- c(4)
# Start the clock!
ptm <- proc.time()

S_pris=0
Price_sum_i_2=0
Price_sum=0
delta=(T)/(b)

Z1=matrix(rnorm(n*b),n, b)
Z2=matrix(rnorm(n*b),n, b)

Z_vol=Z1
Z=rho * Z_vol + sqrt(1-rho^2)*Z2

for (i in 1:n) {
x_t=log(S)
vol_i=V0
S_produkt = 1*S
S_mean_GEO=0

for (j in 1:b) {

vol_i=V0
x_t=x_t + (rf - 0.5*max(vol_i,0))*delta+(max(vol_i,0))**(0.5)*delta**(0.5)*Z[i,j]
vol_i = vol_i + kappa * (theta - max(vol_i,0)) * delta + sigma_vol * (max(vol_i,0))**(0.5)*(delta)**(0.5)*Z_vol[i,j]
S_produkt = S_produkt*exp(x_t)
}

S_mean_GEO=S_produkt^(1/(b+1))

C = exp(-1*rf*T)*max(S_mean_GEO-K,0)
Price_sum=Price_sum+C
Price_sum_i_2=Price_sum_i_2+C^2

if (i==n) {
print("*STANDARD HESTON MONTE CARLO MODELLEN ASIAN GEO*")
values[1]=Price_sum/n #PRIS
print(paste0("Option Price: ", values[1]))
values[2]=(Price_sum_i_2/n)-(Price_sum/n)^2 #VARIANS
values[3]=sqrt(values[2])/sqrt(n) #STD
print(paste0("Std on estimat: ", values[3]))
}
}
# Stop the clock
values[4]=(proc.time() - ptm)[3]
print(paste0("TIme: ", values[4]))

return(list(Price=values[1],Varians=values[2],Std.dev=values[3],Time=values[4]))
}

• Can you post your Monte Carlo code, too? If you are using Euler discretization that could be a problem. – Brian B Aug 31 '15 at 14:57
• I have added the Monte Carlo code. I am using Euler discretization. Why is that a problem? – Rasmus Aug 31 '15 at 15:07
• The skew of stochastic vol increments is quite, high, making Euler scheme prices inaccurate...though you have 6 digits of precision on your MC price estimate, it is likely that the true price is significantly different. – Brian B Aug 31 '15 at 18:22
• your delta may be too large. Euler discretization is good enough with daily time slices. Try by setting b to 500, (i.e., delta = 1/500). – Gordon Aug 31 '15 at 18:53
• I dont think that is the problem. Because this setup converge to the closed-form solution when I price a EU call option. – Rasmus Sep 2 '15 at 15:40