Heston Stochastic Volatility Stochastic Process
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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