I’m trying to implement a moving average model based on this example from stan:
import theano.tensor as tt
lags = 2
nlags = Nlags(lags)
with pm.Model() as ar:
W = pm.Normal( 'W', mu=0, sd=10, shape=( lags, 1 ) )
# eps is the series of unobserved errors
# using a list for convenience for now
eps = price_train.tolist() # errors initialized with price_train
# initialize first and second error manually
i = tt.stack( [[0], [0]], axis=1 )
eps[0] = eps[0] - tt.dot( i, W ).squeeze()
i = tt.stack( [[0], eps[0]], axis=1 )
eps[1] = eps[1] - tt.dot( i, W ).squeeze()
# error is defined as the difference between actual value and prediction
for i in range(len(eps)-lags):
stride_end = i+lags
# take the next stride of 'lags' elements
# and shove them into a theano tensor
e = tt.stack( eps[i:stride_end], axis=1 )
# make a prediction
y_hat = tt.dot( e, W ).squeeze()
# update the next epsilon
# eps was initialized with observed values
eps[ stride_end ] = eps[ stride_end ] - y_hat
# transform eps back to theano tensor
eps = tt.stack( eps, axis=0 )
# transform the data to form a supervised learning problem
eps_lags = nlags.build(input=eps)
z_t = tt.dot( eps_lags, W )
sigma = pm.InverseGamma( 'sigma', alpha=3 ,beta=1 )
ar_out = pm.Normal( 'out', mu=z_t, sd=sigma, observed=nlags.trim(price_train) )
Nlags class takes a vector of time series data and maps each point [x_t] to [x_t-1, x_t-2,…,x_t-lags].
It has already been successfully used during implementation of AR(N) model.
# adds a dimension and stacks n consecutive lagged copies of X along it
# lags with respect to the first axis that encodes examples
# trims X to remove missing data
class Nlags:
def __init__( self, lags ):
self.lags=lags
def trim( self, trimmed ):
return trimmed[self.lags:, ...]
def build( self, *, input ):
results = [ tt.roll( input, i+1, axis=0 ) for i in range(self.lags) ]
result = tt.concatenate( results, axis=-1 )
return self.trim(result)
Model definition runs without errors, but python interpreter crashes as soon as I try to sample from it using NUTS.
Data is stored in price_train vector and it contains stock price which had been previously differenced (with NaNs removed).
Can somebody either point my mistake out, or direct me towards a working implementation of MA model in PYmc3? I haven’t been able to find it anywhere.