Help with Multivariate SDE Timeseries

Questions
33

I think I was able to implement the workaround, but the sampling is extremely slow.


import numpy as np
import matplotlib.pyplot as plt
import pymc3 as pm
from pymc3.distributions.timeseries import EulerMaruyama
import theano.tensor as tt
from theano import scan

returns = np.genfromtxt(pm.get_data("SP500.csv"))

def get_volatility(x,initial_vol,alpha,beta,sigma):
    x = x[:-1]

    def volatility_update(x, vol, a, b,sigma):
        return a + beta * vol + sigma * x

    vol, _ = scan(fn=volatility_update,
                  sequences=[x],
                  outputs_info=[initial_vol],
                  non_sequences=[alpha, beta,
                                 sigma])
    return tt.concatenate([[initial_vol], vol])


with pm.Model() as sp500_model:

    alpha=pm.Normal('alpha',0,sd=100)
    beta=pm.Normal('beta',0,sd=10)
    sigma_v=pm.InverseGamma('sigma_v',2.5,0.1)
    error_v=pm.Normal('error_v',mu=0,sd=1,shape=len(returns))
    rho=pm.Uniform('rho',-.9,.9)


    logV=get_volatility(error_v,.01,alpha,beta,sigma_v)

    volatility_process = pm.Deterministic('volatility_process', pm.math.exp(.5*logV))
    r=pm.Normal('r',mu=volatility_process*rho*error_v,sd=volatility_process*np.sqrt(1-rho*rho),observed=returns)



with sp500_model:
    trace = pm.sample(2000,chains=1,cores=1)

pm.traceplot(trace, [alpha,beta,sigma_v]);
    
fig, ax = plt.subplots(figsize=(14, 8))
ax.plot(returns)
ax.plot(trace['volatility_process',::5].T, 'r', alpha=.03);
ax.set(xlabel='time', ylabel='returns')
ax.legend(['S&P500', 'stoch vol']);
2 Likes