\section{Heston Model}\label{sec:HM}
As pointed before, the BSM framework was important because it presented a close form solution to the problem of Option Pricing. Although its assumptions are unrealistic, the message is clear and calculations insightful. Somehow this model is the benchmark before which all others measure themselves.
One of those models, trying to overcome the restriction of constant volatility along the life period of the option, is the Heston model, introduced by Heston in \cite{Heston1993}. In this paper he proposes a stochastic volatility model in which the volatility, instead of being constant, follows a square-root diffusion. Mathematically, the model is stated by the system os SDE
\begin{align*} \displaystyle\frac{dS_t}{S_t} &= \mu dt +\sqrt{V_t} dZ_t^1\\ \displaystyle dV_t &= \alpha(b-V_t)dt + \sigma\sqrt{V_t}dZ_t^2, \end{align*}where $Z_t^1$ and $Z_t^2$ are correlated BM, i.e. $dZ_t^1\cdot dZ_t^2=\rho dt$, $\alpha$ is the mean reversion rate, or the rate to which the current volatility is reversing to the long run variance $b$ and $\sigma$ is the volatility of volatility. There are some stability results of the solutions of this model, but those are beyond the scope of this analysis.
The Euler scheme on which we will perform the simulations is
\begin{align*} S_{t_{i+1}} &= S_{t_{i}}(1+\mu\Delta) + S_{t_{i}}\cdot\sqrt{V_{t_{i}}}\sqrt{\Delta}\, G_{1,i+1}^S\\ V_{t_{i+1}} &= V_{t_{i}} + \alpha(b-V_{t_{i}})\Delta + \sigma\cdot\sqrt{V_{t_{i}}}\sqrt{\Delta}\, \left(\rho G_{1,i+1}^S+\sqrt{1-\rho^2} G_{2,i+1}^S\right). \end{align*}The first thing we have to do is to generate two correlated standard normal distributions using the following result.
\begin{lemma}\label{lem:corrgauss} Given two independent standard Gaussian arrays $G_1$ and $G_2$, and $\rho\in[-1,1]$. Then $$ G = \rho\, G_1 +\sqrt{1-\rho^2}\,G_2 $$ is $\rho$-correlated with $G_1$. \end{lemma}using Plots
function maxplus(x)
out = 0
if x > 0
out = x
end
return out
end
nsteps = 500
initial_t = 0
final_t = 30/365
delta = (final_t - initial_t)/nsteps
sqrdelta = sqrt(delta)
rho = 0.1
gauss1 = randn(nsteps)
gauss2 = rho * gauss1 + sqrt(1-rho^2) * randn(nsteps)
tempo = linspace(initial_t , final_t , nsteps)
vzero = 0.25
szero = 9.51
mu = 0.1 * delta
alpha = 0.1 #1.98
longrunvol = 0.01 # b parameter in Heston model
volvol = 0.15 # sigma parameter in Heston model
# the initialization of the volatility array
v_values = zeros(nsteps)
v_values[1] = vzero
for i in 2:nsteps
summand1 = v_values[i-1] + alpha *(longrunvol-v_values[i-1]) * delta
summand2 = volvol * sqrt(maxplus(v_values[i-1])) * sqrdelta * gauss2[i]
v_values[i] = summand1 + summand2
end
s_values = szero * (1 + mu * tempo + sqrdelta * map(sqrt,map(maxplus,v_values)).* gauss1 )
println("The value of volatility at the end of the period is $(v_values[nsteps])")
println("The asset price simulation at the end of period is $(s_values[nsteps])")
plot(tempo,[v_values,s_values],legend = false , layout = (2,1))
Now we can make this calculation automatic defining a function to simulate a path of the Heston model.
function hestonpath(initial_t , final_t , nsteps , vzero , szero , gausscorr , drift, alpha , longrunvol , volvol)
delta = (final_t - initial_t)/nsteps
sqrdelta = sqrt(delta)
mu = drift * delta
rho = gausscorr
gauss1 = randn(nsteps)
gauss2 = rho * gauss1 + sqrt(1-rho^2) * randn(nsteps)
tempo = linspace(initial_t , final_t , nsteps)
# the initialization of the volatility array
v_values = zeros(nsteps)
v_values[1] = vzero
for i in 2:nsteps
summand1 = v_values[i-1] + alpha *(longrunvol-v_values[i-1]) * delta
summand2 = volvol * sqrt(maxplus(v_values[i-1])) * sqrdelta * gauss2[i]
v_values[i] = summand1 + summand2
end
s_values = szero * (1 + mu * tempo + sqrdelta * map(sqrt,map(maxplus,v_values)).* gauss1 )
return v_values , s_values
end
nsteps = 500
initial_t = 0
final_t = 30/365
vzero = 0.25
szero = 9.51
gausscorr = 0.1
drift = 0.1
alpha = 0.1 #1.98
longrunvol = 0.01 # b parameter in Heston model
volvol = 0.15 # sigma parameter in Heston model
v_values , svalues = hestonpath(initial_t , final_t , nsteps , vzero , szero , gausscorr , drift, alpha , longrunvol , volvol)
tempo = linspace(initial_t , final_t , nsteps)
println("The value of volatility at the end of the period is $(v_values[nsteps])")
println("The asset price simulation at the end of period is $(s_values[nsteps])")
plot(tempo,[v_values,s_values],legend = false , layout = (2,1))
Next we compute a bunch of paths in order to simulate the values of the Call and Put Options.
initial_t = 0
final_t = 30/365
nsteps = 250
vzero = 0.15
szero = 9.51
strike = 9.0
gausscorr = 0.1
drift = 0.01
alpha = 0.01 #1.98
longrunvol = 0.01 # b parameter in Heston model
volvol = 0.25 # sigma parameter in Heston model
npaths = 100000
volpath , assetpath = hestonpath(initial_t , final_t , nsteps , vzero , szero , gausscorr , drift, alpha , longrunvol , volvol)
# The call and put values initialization
callvalue = maxplus(assetpath[nsteps]-strike)
putvalue = maxplus(strike - assetpath[nsteps])
for i in 1:npaths
volpath , assetpath = hestonpath(initial_t , final_t , nsteps , vzero , szero , gausscorr , drift, alpha , longrunvol , volvol)
# before increment i update mean_path
callvalue = (1-1/(i+1)) * callvalue + (1/(i+1)) * maxplus(assetpath[nsteps]-strike)
putvalue = (1-1/(i+1)) * putvalue + (1/(i+1)) * maxplus(strike - assetpath[nsteps])
end
callvalue = exp(-drift * (final_t - initial_t)) * callvalue
putvalue = exp(-drift * (final_t - initial_t)) * putvalue
println("The price of European Call Option is $(callvalue)")
println("The price of European Put Option is $(putvalue)")
(Heston, 1993) Heston Steven L, ``A Closed-Form Solution for Options with Stochastic Volatility with Applications to Bond and Currency Options'', Review of Financial Studies, vol. 6, number 2, pp. 327-43, 1993.