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( [, ], axis=1 ) eps = eps - tt.dot( i, W ).squeeze() i = tt.stack( [, eps], axis=1 ) eps = eps - 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.