The meaning of the single variables and parameters can now be inferred easily from the
discussion of the geometric Brownian motion and the square-root diffusion. The parameter
represents the instantaneous correlation between the two standard Brownian motions
. This allows us to account for a stylized fact called the leverage effect, which in
essence states that volatility goes up in times of stress (declining markets) and goes down
in times of a bull market (rising markets).
Consider the following parameterization of the model:
In [ 32 ]: S0 = 100.
r = 0.05
v0 = 0.1
kappa = 3.0
theta = 0.25
sigma = 0.1
rho = 0.6
T = 1.0
To account for the correlation between the two stochastic processes, we need to determine
the Cholesky decomposition of the correlation matrix:
In [ 33 ]: corr_mat = np.zeros(( 2 , 2 ))
corr_mat[ 0 , :] = [1.0, rho]
corr_mat[ 1 , :] = [rho, 1.0]
cho_mat = np.linalg.cholesky(corr_mat)
In [ 34 ]: cho_mat
Out[34]: array([[ 1. , 0. ],
[ 0.6, 0.8]])
Before we start simulating the stochastic processes, we generate the whole set of random
numbers for both processes, looking to use set 0 for the index process and set 1 for the
volatility process:
In [ 35 ]: M = 50
I = 10000
ran_num = npr.standard_normal(( 2 , M + 1 , I))
For the volatility process modeled by the square-root diffusion process type, we use the
Euler scheme, taking into account the correlation parameter:
In [ 36 ]: dt = T / M
v = np.zeros_like(ran_num[ 0 ])
vh = np.zeros_like(v)
v[ 0 ] = v0
vh[ 0 ] = v0
for t in range( 1 , M + 1 ):
ran = np.dot(cho_mat, ran_num[:, t, :])
vh[t] = (vh[t - 1 ] + kappa * (theta - np.maximum(vh[t - 1 ], 0 )) * dt
+ sigma * np.sqrt(np.maximum(vh[t - 1 ], 0 )) * np.sqrt(dt)
* ran[ 1 ])
v = np.maximum(vh, 0 )
For the index level process, we also take into account the correlation and use the exact
Euler scheme for the geometric Brownian motion:
In [ 37 ]: S = np.zeros_like(ran_num[ 0 ])
S[ 0 ] = S0
for t in range( 1 , M + 1 ):
ran = np.dot(cho_mat, ran_num[:, t, :])
S[t] = S[t - 1 ] * np.exp((r - 0.5 * v[t]) * dt +
np.sqrt(v[t]) * ran[ 0 ] * np.sqrt(dt))
This illustrates another advantage of working with the Euler scheme for the square-root
diffusion: correlation is easily and consistently accounted for since we only draw standard
normally distributed random numbers. There is no simple way of achieving the same with
