Does someone know why the code below doesn’t terminate at:
model = make_stochastic_volatility_model(sp500_df)
This line of code is from the following in which I try to estimate the parameters of a stochastic volatility model:
import os as os
import sys as sys
path = os.path.dirname(os.path.realpath('__file__'))
#import pandas_market_calendars as mcal
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc3 as pm
from pymc3.distributions.timeseries import GaussianRandomWalk
import pandas_datareader as pdr
import datetime
import pprint
import seaborn as sns
sns.set(palette = 'bright',color_codes=True)
# obtain sp500 log returns (scaled by 100)
start_date = datetime.datetime(1987, 12, 31)
end_date = datetime.datetime(2000, 12, 31)
sp500_df = pdr.DataReader('^GSPC', 'yahoo',start_date,end_date)
sp500_df["returns"] = sp500_df["Close"]/sp500_df["Close"].shift(1)
sp500_df.dropna(inplace=True)
sp500_df["log_returns"] = 100*np.log(sp500_df["returns"])
# plot sp500 log returns
fig, ax = plt.subplots(figsize=(16, 8))
sp500_df['log_returns'].plot(y="change", label='S&P 500', ax=ax)
ax.set(xlabel='time', ylabel='returns')
def make_stochastic_volatility_model(data):
with pm.Model() as model:
obs = data.shape[0]
# prior distribution of parameters
mu = pm.Normal(name='mu',mu=0,sigma=1)
epsilon = pm.Normal(name='epsilon',mu=0,sigma=1,shape=obs)
eta = pm.Normal(name='eta',mu=0,sigma=1,shape=obs)
beta_help = pm.Beta('beta_help',alpha=20,beta=1.5) #testval=0.5125
beta = 2*beta_help-1
alpha_help = pm.Normal('alpha_help',mu=0,sigma=5)
alpha = alpha_help*(1-beta) #HalfNormal?
varH = pm.InverseGamma('varH',alpha=2.5,beta=0.025) #testval=0.025
rho = pm.Uniform('rho',-1,1) #testval=-0.4
r_list = []
ln_h_list = []
for t in range(obs):
if t == 0:
h_current = pm.HalfNormal(f'h_{t}',sigma=1.0)
ln_h_current = pm.Deterministic(f'ln_h_{t}',np.log(h_current))
else:
ln_h_current = pm.Normal(f'ln_h_{t}',
mu=alpha+(beta)*ln_h_current,
sigma=np.sqrt(varH))
h_current = pm.Deterministic(f'h_{t}',np.exp(ln_h_current))
ln_h_list.append(ln_h_current)
r_current = pm.Normal(f'r_{t}',mu=mu,
sigma=np.sqrt(h_current),
observed=data['log_returns'][t])
r_list.append(r_current)
return model
model = make_stochastic_volatility_model(sp500_df)
# Check model
# pm.model_to_graphviz(model)
with model:
# Use Maximum A Posteriori (MAP) optimisation
# as initial value for MCMC
# start = pm.find_MAP()
# Use the No-U-Turn Sampler
step = pm.NUTS()
# Calculate the trace
trace = pm.sample(2000, step,
random_seed=42, progressbar=True, cores=1, tune=1000)
