How to incorporate a lagged value?

Questions
1

Hi all,

Currently, I am estimating a stochastic model. This model requires a lagged value, for example,

  1. y_{t} = mu + x_{t}*epsilon_{1,t}
  2. x_{t} = alpha +beta_0x_{t-1} + beta_1 epsilon{1, t-1} + epsilon_{2,t}

With beta_1 == 0, I used the pm.ar(.) function but this is not possible when beta_1 is, for example, uniform(0,1).

How is it possible to incorporate this epislon{1, t-1}? Is it possible to remember the lagged values in pymc3?

Thank you in advance!

2

Hi there,

I am working on almost the same problem as yours. You could have a look here: Stochastic volatility model estimation with MCMC doesnt terminate

Basically, I created a list and a for loop to simulate the time series and store them such that I could reuse them later on (as it was more flexible than pm.AR()). You could perceive alpha +beta_0 x_{t-1} + beta_1 epsilon{1, t-1} + epsilon_{2,t} as two separate problems, that is, alpha +beta_0 *x_{t-1} + epsilon_{2,t}, and beta_1 epsilon{1, t-1}. The former could be performed with what I’m already trying to do, and the latter by creating an epsilon list and then as follows:

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 = []
        epsilon_list = []
        eta_list = []
        epsilon_list.append(epsilon)
        eta_list.append(eta)
        
        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.Deterministic(f'ln_h_{t}',
                                                alpha + 
                                                beta*ln_h_current +
                                                np.sqrt(varH)*rho*epsilon_list[t-1] +
                                                np.sqrt(varH)*np.sqrt(1-rho**2)*eta_list[t])
                h_current = pm.Deterministic(f'h_{t}',np.exp(ln_h_current))

            ln_h_list.append(ln_h_current)
            r_current = pm.Deterministic(f'r_{t}',mu + np.sqrt(h_current)*epsilon_list[t],
                                         observed=data['log_returns'][t])
            r_list.append(r_current)
    return model

The model I’m trying to replicate (without delta):
image

Take note though, my code (see link) doesn’t run yet (also I havent tested the current one yet). Furthermore, I assume it is possible to remember lagged values from what I saw in the containers section 4.4:
https://pymc-devs.github.io/pymc/modelbuilding.html#containers

Finally, regarding readability, you could replace the current values, e.g. ln_h_current, by ln_h_list[t] after specifying ln_h_list = list(np.zeros(obs)) instead of ln_h_list =[ ].