Stochastic Volatility Model

Reference: Python for Finance by Yves Hilpisch, O'Reilly Media, Chapter 10, Stochastics

Multidimensional Stochastic Processes

import numpy as np
import numpy.random as npr
import matplotlib.pyplot as plt
import pandas as pd
import scipy as sp
import scipy.linalg 
import scipy.stats

Features of the implied volatility surface such as volatility smile and skew indicate that implied volatility varies with respect to strike price and expiry. Is it also a function of the underlying price?

The Heston SV model assumes that the variance is a random process that exhibits:

• a tendency to revert towards a long-term mean $\omega$ at a rate $\theta$
• a volatility proportional to the square root of its level
• and whose source of randomness is negatively correlated with the randomness of the underlying's price process, $\rho<0$.

Geometric Brownian motion SDE:

$dS_t = rS_tdt + \sqrt{v_t} S_t dZ_t^1$, where $v_t = \sigma_t^2 $

$ S_t = S_{t-dt} exp [ (r- v_t/2)dt + \sqrt{v_t} \sqrt{dt} z_t^1 ] $

AR1 process for the variance of the process:

$dv_t = k_v(\theta_v-v_t)dt + \sigma_v \sqrt{v_t} dZ_t^2$, where $\theta$ is the long-term mean of the variance process.

$dZ_t^1dZ_t^2 = \rho dt$

I have read but not verified that the condition $ k_v \theta_v > \sigma^2_v/2 $ ensures positivity of the variance process.

#Simulate N process paths with arbitrary parameters: 
I=100
v0=0.1
sigmav = 0.1 # std deviation factor of the variance process 
k = 3.0
th = 0.25 #theta
r = 0.05
T = 1
M=50
dt = T/M
S0=100.
rho = -0.6

To account for the correlation between two stochastic processes, we need to perform the Choleski decomposition of the correlation matrix. A hermitic(...) positive-definite matrix $A = LL^*$, the product of two triangular matrices. If $A$ is symmetric (which is true iff $A$ is real), the factorization may be written $A = LL^T$.

corr_mat = np.zeros((2,2))
corr_mat[0,:] = [1.0, rho]
corr_mat[1,:] = [rho, 1.0]
#L = scipy.linalg.cholesky(corr_mat, lower=True)
#U = scipy.linalg.cholesky(corr_mat, lower=False)
chol_mat = np.linalg.cholesky(corr_mat); chol_mat

array([[ 1. ,  0. ],
       [-0.6,  0.8]])

We generate two sets of correlated random numbers, then we use the Euler scheme to take into account the correlation.

rand = npr.standard_normal((2,M+1, I))
v = np.zeros_like(rand[0]) # We need to impose v positive, so we use an auxiliary vh.
vh = np.zeros_like(v)
v[0] = v0 
vh[0] = v0

for t in range(1, M+1):
    z = np.dot(chol_mat, rand[:,t,:]) #2-dim process
    clamped = np.maximum(vh[t-1],0)
    vh[t] = vh[t-1] + k*(th-clamped)*dt + sigmav*np.sqrt(clamped)*np.sqrt(dt)*z[1]
    v = np.maximum(vh,0)    

S = np.zeros_like(rand[0]) 
S[0] = S0 
for t in range(1, M+1):
    z = np.dot(chol_mat, rand[:,t,:])
    S[t] = S[t-1]*np.exp((r-0.5*v[t])*dt + np.sqrt(v[t])*np.sqrt(dt)*z[0])      

plt.plot(S[:,:20], lw = 1)
plt.xlabel('# periods'); plt.ylabel('Index Level')
plt.grid(True)
#To do: plot on log scale

plt.plot(v[:,:25], lw = 1)
plt.xlabel('# periods'); plt.ylabel('Index Level')
plt.grid(True)

#State variable
# Terminal Value: 
ST = S[-1]
#The state variable follows a Lognormal distribution.
fig, ax1 = plt.subplots()#
ax1.hist(ST, bins=50); ax1.set_title('Process Values - Lognormal Distribution')

Text(0.5, 1.0, 'Process Values - Lognormal Distribution')

#Variance of the process:
# Terminal Value: 
vT = v[-1]
#The state variable follows a Lognormal distribution.
fig, ax1 = plt.subplots()#
ax1.hist(vT, bins=50); ax1.set_title('Process Variance - Normal Distribution \n with mean theta')

Text(0.5, 1.0, 'Process Variance - Normal Distribution \n with mean theta')

scipy.stats.describe((S[-1]))

DescribeResult(nobs=100, minmax=(30.529011677930516, 245.42697805069005), mean=104.35388861375988, variance=2375.9483710690743, skewness=0.8361075659188333, kurtosis=0.07602336542600119)

scipy.stats.describe((v[-1]))

DescribeResult(nobs=100, minmax=(0.19283207326610985, 0.2958076774498704), mean=0.24321525147797518, variance=0.0004992162361855502, skewness=0.11176406215837527, kurtosis=-0.19539346470800556)

Expected compounded average return

$CAGR = (S_T/S_0)^{1/T}-1 \approx \frac{log(S_T/S_0)}{T}$

LOG = np.log(np.mean(ST)/S0)/M
CAGR = (np.mean(ST)/S0)**(1/M)-1
print(round(LOG,5), round(CAGR,5))

0.00085 0.00085

#Predicted returns are lower than the returns predicted by the GBM model.
# In this case, average returns are negative even though the drift term is positive.
log_Rt = np.log(S[1:,:]/S[:-1,:])
print(round(np.mean((log_Rt)),5), round(np.std((log_Rt)), 5))

-0.00127 0.06411

Appendix: Real Values

import pandas_datareader as pr
from pandas_datareader import data

start_date='6/30/2010'; end_date='08/31/2020'
google = data.DataReader('^GSPC', 'yahoo', start_date, end_date)

plt.plot(google['Close'])
plt.ylabel('Google Close Price')
plt.grid(True)

google['logReturn']=np.log(google['Close']/google['Close'].shift(1))
google['logReturn'].describe()

count    2560.000000
mean        0.000478
std         0.010942
min        -0.127652
25%        -0.003314
50%         0.000653
75%         0.005234
max         0.089683
Name: logReturn, dtype: float64

goog = np.array(google['logReturn'][1:])

The asymptotic properties of the SV model, necessary for the Bird's eye view calibration, are less straight-forward than for the previous models.

Is there a way to obtain $\rho$ directly from the original price series, instead of the rolling series? I need to brush up the theory. For now:

rollVol = google["logReturn"].rolling(252).var() #Variance, not standard deviation!
rollmean = google["logReturn"].rolling(252).mean()

Y = rollVol[~np.isnan(rollVol)]
X = rollmean[~np.isnan(rollmean)]

rho = np.corrcoef(X, Y)[1,0]; rho

-0.4330162712248781

MeanV = (np.mean(Y**2)); round(MeanV,10) #Mean rolling variance
print(round(np.var(goog), 6), round(MeanV,6))
#Two different measures of stock price variance are not very far from eachother.

0.00012 0.0

However, $\theta$ is NOT simply the variance of the price process. I am not even sure if there is a closed-form formula for it (there probably is). Can we find both theta and k from the coefficients of the AR(1) equation for stochastic volatility?

$v_t = k\theta + (1-k)v_{t-1}+\epsilon_t$

Unfortuantely, the instantaneous variance of the price process is not obesrvable. I have no better idea than to aproximate it by rolling variance. Another caveat for this estimation is that the error term $\epsilon$ is heteroscedastic, it rises proportionately with the instantaneous standard deviation of the process.

YY = Y[2:] #dependent vt
l=len(YY); l
YX = np.vstack((np.ones(l), Y[1:-1])).T #independent variable v(t-1)

a,b = np.linalg.lstsq(YX, YY, rcond=None)[0]
#pd.DataFrame(result)
print("alpha:   beta: ")
print(round(a,5),round(b,5))

alpha:   beta: 
-0.0 1.00278

#alpha = k*theta
#beta = 1-k
k = 1-b
th = a/k
print("k:      theta: ")
print(round(k,5),round(th,5))

k:      theta: 
-0.00278 4e-05

v0 = Y[0]
sigmav = np.std(Y) #std of rolling variance...
r = np.mean(goog)+th/2 #the drift term for the price process
print(round(r,5), round(sigmav,5))

0.0005 9e-05

S0 = google['Close'][0];S0

1030.7099609375

corr_mat = np.zeros((2,2))
corr_mat[0,:] = [1.0, rho]
corr_mat[1,:] = [rho, 1.0]
chol_mat = np.linalg.cholesky(corr_mat); chol_mat

array([[ 1.        ,  0.        ],
       [-0.43301627,  0.9013861 ]])

rand = npr.standard_normal((2,M+1, I))
v = np.zeros_like(rand[0]) # We need to impose v positive, so we use an auxiliary vh.
vh = np.zeros_like(v)
v[0] = v0 
vh[0] = v0

for t in range(1, M+1):
    z = np.dot(chol_mat, rand[:,t,:]) #2-dim process
    clamped = np.maximum(vh[t-1],0)
    vh[t] = vh[t-1] + k*(th-clamped)*dt + sigmav*np.sqrt(clamped)*np.sqrt(dt)*z[1]
    v = np.maximum(vh,0)    

S = np.zeros_like(rand[0]) 
S[0] = S0 
for t in range(1, M+1):
    z = np.dot(chol_mat, rand[:,t,:])
    S[t] = S[t-1]*np.exp((r-0.5*v[t])*dt + np.sqrt(v[t])*np.sqrt(dt)*z[0])      

plt.plot(S)
plt.xlabel('# periods'); plt.title('Bird Eye Calibration'); plt.ylabel('Price Process is GBM with SV')
plt.grid(True)

plt.plot(v)
plt.xlabel('# periods'); plt.title('Bird Eye Calibration'); plt.ylabel('Stochastic Volatility Process')
plt.grid(True)

#The results are all over the place. 
#Next I will attempt to obtain more exact parameters by minimizing the squared difference between observed index values and predicted values, for the two-dimensional process.

import scipy as sp
from scipy import stats

$\theta^\star = \min_\theta \sum_i \|y_i-f(x_i,\theta)\|^2 = \min_\theta g(\theta)$

$f(\theta)= Realized S_T$

$g =\|y - f(\theta)\|^2$

GPrice = np.array(google['Close'])

S0 = GPrice[0]

theta = [rho, r, k, th, sigmav]#SV parameters
I = 1
M=100

def f(theta):
    #M = len(GPrice)
    corr_mat = np.zeros((2,2))
    corr_mat[0,:] = [1.0, theta[0]]
    corr_mat[1,:] = [theta[0], 1.0]
    chol_mat = np.linalg.cholesky(corr_mat); chol_mat

    rand = npr.standard_normal((2,M+1, I))
    v = np.zeros_like(rand[0]) 
    vh = np.zeros_like(v)
    v[0] = v0 
    vh[0] = v0
    for t in range(1, M+1):
        z = np.dot(chol_mat, rand[:,t,:]) #2-dim process
        clamped = np.maximum(vh[t-1],0)
        vh[t] = vh[t-1] + theta[2]*(theta[3]-clamped)*dt + theta[4]*np.sqrt(clamped)*np.sqrt(dt)*z[1]
        v = np.maximum(vh,0)    
    
    S = np.zeros_like(rand[0]) 
    S[0] = S0 
    for t in range(1, M+1):
        z = np.dot(chol_mat, rand[:,t,:])
        S[t] = S[t-1]*np.exp((theta[1]-0.5*v[t])*dt + np.sqrt(v[t])*np.sqrt(dt)*z[0])   
        
    return S

ordd = 2

g = lambda theta: np.linalg.norm(GPrice[:99]-f(theta), ord=ordd)**ordd

import scipy.optimize
thetaopt = scipy.optimize.fmin(func=g, x0=[rho, r, k, th, sigmav],ftol=1e-7,xtol=1e-7)
print('Parameter calibration for market prices: \n', thetaopt)

Warning: Maximum number of function evaluations has been exceeded.
Parameter calibration for market prices: 
 [-4.53318954e-01  5.15164249e-04 -2.78723551e-03  4.21371711e-05
  9.20919314e-05]

print(g([rho, r, k, th, sigmav]))
print(g(thetaopt))

118672664.72905757
123315738.35774077

plt.plot(GPrice[:101], f(thetaopt), 'r*')

[<matplotlib.lines.Line2D at 0x1a25470f28>]

#Simulate the calibrated GBM path: 
I=100
M = len(GPrice)
Sopt = f(thetaopt)
plt.plot(Sopt)
plt.ylabel('Fmin Calibration')
plt.grid(True)

The Fmin minimization leads to poor results (all the parameters are close to zero). I need to use the "earth moving" minimization of the distance between two stochastic processes (TO DO).