Heston Stochastic Volatility Stochastic Process
import math | |
import numpy | |
import random | |
import decimal | |
import scipy.linalg | |
import numpy.random as nrand | |
import matplotlib.pyplot as plt | |
""" | |
Note that this Gist uses the Model Parameters class found here - https://gist.github.com/StuartGordonReid/f01f479c783dd40cc21e | |
""" | |
def cox_ingersoll_ross_heston(param): | |
""" | |
This method returns the rate levels of a mean-reverting cox ingersoll ross process. It is used to model interest | |
rates as well as stochastic volatility in the Heston model. Because the returns between the underlying and the | |
stochastic volatility should be correlated we pass a correlated Brownian motion process into the method from which | |
the interest rate levels are constructed. The other correlated process is used in the Heston model | |
:param param: the model parameters objects | |
:return: the interest rate levels for the CIR process | |
""" | |
# We don't multiply by sigma here because we do that in heston | |
sqrt_delta_sigma = math.sqrt(param.all_delta) * param.all_sigma | |
brownian_motion_volatility = nrand.normal(loc=0, scale=sqrt_delta_sigma, size=param.all_time) | |
a, mu, zero = param.heston_a, param.heston_mu, param.heston_vol0 | |
volatilities = [zero] | |
for i in range(1, param.all_time): | |
drift = a * (mu - volatilities[i-1]) * param.all_delta | |
randomness = math.sqrt(volatilities[i - 1]) * brownian_motion_volatility[i - 1] | |
volatilities.append(volatilities[i - 1] + drift + randomness) | |
return numpy.array(brownian_motion_volatility), numpy.array(volatilities) | |
def heston_model_levels(param): | |
""" | |
NOTE - this method is dodgy! Need to debug! | |
The Heston model is the geometric brownian motion model with stochastic volatility. This stochastic volatility is | |
given by the cox ingersoll ross process. Step one on this method is to construct two correlated GBM processes. One | |
is used for the underlying asset prices and the other is used for the stochastic volatility levels | |
:param param: model parameters object | |
:return: the prices for an underlying following a Heston process | |
""" | |
assert isinstance(param, ModelParameters) | |
# Get two correlated brownian motion sequences for the volatility parameter and the underlying asset | |
# brownian_motion_market, brownian_motion_vol = get_correlated_paths_simple(param) | |
brownian, cir_process = cox_ingersoll_ross_heston(param) | |
brownian, brownian_motion_market = heston_construct_correlated_path(param, brownian) | |
heston_market_price_levels = [param.all_s0] | |
for i in range(1, param.all_time): | |
drift = param.gbm_mu * heston_market_price_levels[i - 1] * param.all_delta | |
vol = cir_process[i - 1] * heston_market_price_levels[i - 1] * brownian_motion_market[i - 1] | |
heston_market_price_levels.append(heston_market_price_levels[i - 1] + drift + vol) | |
return numpy.array(heston_market_price_levels), numpy.array(cir_process) | |
def heston_construct_correlated_path(param, brownian_motion_one): | |
""" | |
This method is a simplified version of the Cholesky decomposition method for just two assets. It does not make use | |
of matrix algebra and is therefore quite easy to implement. | |
:param param: model parameters object | |
:return: a correlated brownian motion path | |
""" | |
# We do not multiply by sigma here, we do that in the Heston model | |
sqrt_delta = math.sqrt(param.all_delta) | |
# Construct a path correlated to the first path | |
brownian_motion_two = [] | |
for i in range(param.all_time - 1): | |
term_one = param.cir_rho * brownian_motion_one[i] | |
term_two = math.sqrt(1 - math.pow(param.cir_rho, 2.0)) * random.normalvariate(0, sqrt_delta) | |
brownian_motion_two.append(term_one + term_two) | |
return numpy.array(brownian_motion_one), numpy.array(brownian_motion_two) | |
def get_correlated_geometric_brownian_motions(param, correlation_matrix, n): | |
""" | |
This method can construct a basket of correlated asset paths using the Cholesky decomposition method | |
:param param: model parameters object | |
:param correlation_matrix: nxn correlation matrix | |
:param n: the number of assets i.e. the number of paths to return | |
:return: n correlated log return geometric brownian motion processes | |
""" | |
assert isinstance(param, ModelParameters) | |
decomposition = scipy.linalg.cholesky(correlation_matrix, lower=False) | |
uncorrelated_paths = [] | |
sqrt_delta_sigma = math.sqrt(param.all_delta) * param.all_sigma | |
# Construct uncorrelated paths to convert into correlated paths | |
for i in range(param.all_time): | |
uncorrelated_random_numbers = [] | |
for j in range(n): | |
uncorrelated_random_numbers.append(random.normalvariate(0, sqrt_delta_sigma)) | |
uncorrelated_paths.append(numpy.array(uncorrelated_random_numbers)) | |
uncorrelated_matrix = numpy.matrix(uncorrelated_paths) | |
correlated_matrix = uncorrelated_matrix * decomposition | |
assert isinstance(correlated_matrix, numpy.matrix) | |
# The rest of this method just extracts paths from the matrix | |
extracted_paths = [] | |
for i in range(1, n + 1): | |
extracted_paths.append([]) | |
for j in range(0, len(correlated_matrix)*n - n, n): | |
for i in range(n): | |
extracted_paths[i].append(correlated_matrix.item(j + i)) | |
return extracted_paths |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment