# How to incorporate a lagged value?

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?

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 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["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

# 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): 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 =[ ].