Hi everyone! Best wishes for 2016!
In this post, I’ll show you how to use ESGtoolkit, for the simulation of Heston stochastic volatility model for stock prices.
If you’re interested in seeing other examples of use of ESGtoolkit, you can read these two posts: the Hull and White short rate model and the 2-factor Hull and White short rate model (G2++).
The Heston model was introduced by Steven Heston’s A closed-form solution for options with stochastic volatility with applications to bonds an currency options, 1993. For a fixed risk-free interest rate 
where 
In this model, under a certain probability, the stock price’s returns on very short periods of time of length 
On the other hand, 
By using this model, one can derive prices for European call options, as described in Calibrating Option Pricing Models with Heuristics. The authors provide a useful function called ‘callHestoncf’, which calculates these prices in R and Matlab.
Here’s the function’s description. I won’t reproduce the function here, please refer to the paper for details:
callHestoncf(S, X, tau, r, v0, vT, rho, k, sigma){
# S = Spot, X = Strike, tau = time to maturity
# r = risk-free rate, q = dividend yield
# v0 = initial variance, vT = long run variance (theta)
# rho = correlation, k = speed of mean reversion (kappa)
# sigma = volatility of volatility
}
Now, it’s time to use ESGtoolkit for Monte Carlo pricing. We are going to price 3 European Call options on 
rm(list=ls()) library(ESGtoolkit) #Initial stock price S0 <- 100 # Number of simulations (feel free to reduce this) n <- 100000 # Sampling frequency freq <- "monthly" # volatility mean-reversion speed kappa <- 0.003 # volatility of volatility volvol <- 0.009 # Correlation between stoch. vol and spot prices rho <- -0.5 # Initial variance V0 <- 0.04 # long-term variance theta <- 0.04 #Initial short rate r0 <- 0.015 # Options maturities horizon <- 15 # Options' exercise prices strikes <- c(140, 100, 60)
For the simulation of the Heston model with ESGtoolkit, we first need to define how to make simulations of the terms
# Simulation of shocks with given correlation set.seed(5) # reproducibility seed shocks <- simshocks(n = n, horizon = horizon, frequency = freq, method = "anti", family = 1, par = rho)
This function provides a list with 2 components, each containing simulated random gaussian increments. Both of these components will be useful for the simulation of 
# Stochastic volatility simulation sim.vol <- simdiff(n = n, horizon = horizon, frequency = freq, model = "CIR", x0 = V0, theta1 = kappa*theta, theta2 = kappa, theta3 = volvol, eps = shocks[[1]]) # Stock prices simulation sim.price <- simdiff(n = n, horizon = horizon, frequency = freq, model = "GBM", x0 = S0, theta1 = r0, theta2 = sqrt(sim.vol), eps = shocks[[2]])
We are now able to calculate options prices, with the 3 different
exercise prices. You’ll need to have imported ‘callHestoncf’, for this piece of code to work.
# Stock price at maturity (15 years)
S_T <- sim.price[nrow(sim.price), ]
### Monte Carlo prices
#### Estimated Monte Carlo price
discounted.payoff <- function(x)
{
(S_T - x)*(S_T - x > 0)*exp(-r0*horizon)
}
mcprices <- sapply(strikes, function(x)
mean(discounted.payoff(x)))
#### 95% Confidence interval around the estimation
mcprices95 <- sapply(strikes, function(x)
t.test(discounted.payoff(x),
conf.level = 0.95)$conf.int)
#### 'Analytical' prices given by 'callHestoncf'
pricesAnalytic <- sapply(strikes, function(x)
callHestoncf(S = S0, X = x, tau = horizon,
r = r0, q = 0, v0 = V0, vT = theta,
rho = rho, k = kappa, sigma = volvol))
results <- data.frame(cbind(strikes, mcprices,
t(mcprices95), pricesAnalytic))
colnames(results) <- c("strikes", "mcprices", "lower95",
"upper95", "pricesAnalytic")
print(results)
strikes mcprices lower95 upper95 pricesAnalytic
1 140 25.59181 25.18569 25.99793 25.96174
2 100 37.78455 37.32418 38.24493 38.17851
3 60 56.53187 56.02380 57.03995 56.91809
From these results, we see that the Monte Carlo prices for the 3 options are
fairly close to the price calculated by using the function ‘callHestoncf’
(uses directly formulas for the prices calculation). The 95% confidence
interval contains the theoretical price. Do not hesitate to change the seed,
and re-run the previous code.
Below are the option prices, as functions of the number of simulations.
The theoretical price calculated by ‘callHestoncf’ is drawn in blue,
the average Monte Carlo price in red, and the shaded region represents
the 95% confidence interval around the mean (the Monte Carlo price).
and 
is a standard brownian motion under a risk-neutral probability. 
, which is constant in Vasicek’s model, is a function constructed so as to correctly match the initial term structure of interest rates.
are market-implied instantaneous forward rates for maturities 
.
Thierry Moudiki's blog 