Question2

A) Heston Monte Carlo

Function to Implement the Heston Monte Carlo for different models in Euler Discretization

HestonMC<-function(case,nIterations,nSteps)
{
    
  S_0=100.0
  K = 100.0
  r = 0.0319
  v_0 = 0.010201
  T = 1.00
  
  dt = T/nSteps
  
  rho = -0.7
  kappa = 6.21
  theta = 0.019
  xi = 0.61
  
  rmse=0
  payoff_sum=0
  sumbias=0
  
  ptm <- proc.time()
  for( j in 1:nIterations) {
    W1=W2=0
    vol_path=v_0
    spot_path=S_0
    vmax=0
    for( i in 1:nSteps){
      switch(case,
        Full={# Full truncation
          vmax= max(vol_path,0)
          vol_path=vol_path+kappa*dt*(theta-vmax)+xi*sqrt(vmax*dt)*W1
          spot_path=spot_path* exp((r - 0.5*vmax)*dt +sqrt(vmax*dt)*W2)
        },
        Partial={# Partial truncation
          vmax= max(vol_path,0)
          vol_path=vol_path+kappa*dt*(theta-vol_path)+xi*sqrt(vmax*dt)*W1
          spot_path=spot_path* exp((r - 0.5*vmax)*dt +sqrt(vmax*dt)*W2)
        },
        HM={# Higham and Mao
          vmax= abs(vol_path)
          vol_path=vol_path+kappa*dt*(theta-vol_path)+xi*sqrt(vmax*dt)*W1
          spot_path=spot_path* exp((r - 0.5*vmax)*dt +sqrt(vmax*dt)*W2)
        },
        Refl={# Reflection
          vmax= abs(vol_path)
          vol_path=abs(vol_path)+kappa*dt*(theta-vmax)+xi*sqrt(vmax*dt)*W1
          spot_path=spot_path* exp((r - 0.5*vmax)*dt +sqrt(vmax*dt)*W2)
        },
        Abso={# Absorption
          vmax= max(vol_path,0)
          vol_path=vmax+kappa*dt*(theta-vmax)+xi*sqrt(vmax*dt)*W1
          spot_path=spot_path* exp((r - 0.5*vmax)*dt +sqrt(vmax*dt)*W2)
        }
        )
      
      X1 = rnorm(1,0,1) 
      X2 = rnorm(1,0,1) 
      W1=X1
      W2= rho*X1+sqrt(1-rho^2)*X2
    }
    payoff=max(spot_path-K,0)
    payoff_sum=payoff_sum+payoff
    rmse=rmse+(payoff-6.8061)^2
    sumbias=sumbias+(max(spot_path-K,0)-6.8061)
    # print(max(spot_path-K,0))
  }
  option_price=(payoff_sum/nIterations)*exp(-r*dt)
  time=as.double((proc.time() - ptm)[3])
  rmse=sqrt(rmse)/nIterations
  sumbias=sumbias/nIterations
  
  list(price=option_price,time_taken=time,rmse=rmse,bias=sumbias)
  # print((proc.time() - ptm)[3])
  # print(option_price)
}

Completed Table

full=    HestonMC(nIterations = 10000,nSteps = 100,case = "Full")
partial= HestonMC(nIterations = 10000,nSteps = 100,case = "Partial")
hm=      HestonMC(nIterations = 10000,nSteps = 100,case = "HM")
refl=    HestonMC(nIterations = 10000,nSteps = 100,case = "Refl")
abso=    HestonMC(nIterations = 10000,nSteps = 100,case = "Abso")

methods=c("Full","Partial","Higham and Mao","Reflection","Absorption")
df=data.frame()
df=rbind(df,full,partial,hm,refl,abso)
df$Methods=methods
print(df)

##      price time_taken       rmse      bias        Methods
## 1 6.971993      17.78 0.07496057 0.1681174           Full
## 2 7.054487      21.22 0.07739652 0.2506375        Partial
## 3 6.977992      15.94 0.07833858 0.1741181 Higham and Mao
## 4 7.311140      15.75 0.08456486 0.5073730     Reflection
## 5 7.112863      15.92 0.08008125 0.3090323     Absorption

2) Bonus: Heston Analytical

source("HestonPrice-1.R")
source("HestonProb.R")

HestonPriceExample = function(
  S = 100,         # Spot price
  K = 100,         # Strike price
  tau = 1,        # Maturity
  r = 0.0319,        # Risk free rate
  div = 0.0,         # Dividend yield
  kappa = 6.21,       # Heston parameter : reversion speed
  sigma = 0.61,     # Heston parameter : volatility of variance
  rho   = -0.7,    # Heston parameter : correlation
  theta = 0.019,    # Heston parameter : reversion level
  v0    = 0.010201,    # Heston parameter : initial variance
  lambda = 0       # Heston parameter : risk preference
)
{
  # Expression for the characteristic function ----
  Trap = 0;   # 0 = Original Heston formulation
              # 1 = Albrecher et al formulation
  
  # rmse=0
  # payoff_sum=0
  # Integration range                
  Lphi = 0.000001;  # Lower limit
  dphi = 0.01;      # Increment
  Uphi = 50;        # Upper limit
  
  # Obtain the Heston put and call
  HPut  = HestonPrice('P',kappa,theta,lambda,rho,sigma,tau,K,S,r,div,v0,Trap,Lphi,Uphi,dphi);
  HCall = HestonPrice('C',kappa,theta,lambda,rho,sigma,tau,K,S,r,div,v0,Trap,Lphi,Uphi,dphi);
  
  # Output the result
  list(Call = HCall, Put = HPut)
}


ptm <- proc.time()
call_price=as.numeric(HestonPriceExample()[1])
time=as.double((proc.time() - ptm)[3])

K=100
rmse=(call_price-6.8061)^2
bias=(max(call_price,0)-6.8061)
rmse=sqrt(rmse)
print(paste("Call Option Price=",round(call_price,2)," Time taken=",round(time,2)," seconds"
            ,"RMSE=",round(rmse,2),"Bias=",round(bias,2)))

## [1] "Call Option Price= 8.48  Time taken= 0.95  seconds RMSE= 1.68 Bias= 1.68"

We can clearly see that the analytical price is faster. Among the different monte carlo methods, we can see that accuracy and computation times are inversly related and looking at the above result Higham and Mao gives the fastest result while the accurate price is from the Partial truncation method.