MA model for time series prediction

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.

You probably need to rewrite the for loop using a tt.scan (otherwise even if you manage to make the code compile you will have serious performance issue.)

I’m still very new to this.
The intention of the for loop was to build the expression tree for creating epsilon values that are sequentially dependent. Are expressions like these slower than tt.scan?

much slower. See here for some more information https://docs.pymc.io/notebooks/PyMC3_tips_and_heuristic.html`

I tried porting the code into scan, here’s my best attempt:

import theano.tensor as tt
lags = 2

nlags = nn.Nlags(lags)
    
with pm.Model() as ar:
    W = pm.Normal( 'W', mu=0, sd=10, shape=( lags, 1 ), dtype='float32' )
    
    y = tt.dmatrix('y')
    y.tag.test_value = price_train
    eps0 = tt.dmatrix('eps0')
    eps0.tag.test_value = np.array([[0],[0]], dtype='float64')
    
    def compute_eps( price, eps_m1, eps_m2, W ):
        # eps(n) = price_train(n) - dot( [ eps(n-1), eps(n-2), ..., eps(n-lags) ], W )
        lagged = tt.stack( [eps_m1, eps_m2], axis=1 )
        return price - tt.dot( lagged, W ).squeeze()
    
    
    eps, updates = theano.scan( fn=compute_eps,
                                 sequences=y,
                                 outputs_info=dict(initial=eps0, taps=[-1,-2]),
                                 non_sequences=W )
    
    calc_eps = theano.function([y, eps0, W], outputs=eps,updates=updates)
    
    # transform the data to form a supervised learning problem
    eps_lags = nlags.build(input=calc_eps(price_train, np.array([[0],[0]], dtype='float64'), W))
    
    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) )

Unfortunately, the loop code is dependent on W, which is a parameter matrix that I’m trying to estimate with PYMC. Theano reports the following error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-49-7a89896c9cba> in <module>
     29 
     30     # transform the data to form a supervised learning problem
---> 31     eps_lags = nlags.build(input=calc_eps(price_train, np.array([[0],[0]], dtype='float64'), W))
     32 
     33     z_t = tt.dot( eps_lags, W )

C:\ProgramData\Anaconda3\envs\PYmc3\lib\site-packages\theano\compile\function_module.py in __call__(self, *args, **kwargs)
    811                         s.storage[0] = s.type.filter(
    812                             arg, strict=s.strict,
--> 813                             allow_downcast=s.allow_downcast)
    814 
    815                     except Exception as e:

C:\ProgramData\Anaconda3\envs\PYmc3\lib\site-packages\theano\tensor\type.py in filter(self, data, strict, allow_downcast)
     85         if isinstance(data, Variable):
     86             raise TypeError(
---> 87                 'Expected an array-like object, but found a Variable: '
     88                 'maybe you are trying to call a function on a (possibly '
     89                 'shared) variable instead of a numeric array?')

TypeError: Bad input argument with name "W" to theano function with name 
"<ipython-input-49-7a89896c9cba>:28" at index 2 (0-based). Expected an array-like object, 
but found a Variable: maybe you are trying to call a function on a (possibly shared) variable 
instead of a numeric array?

Do I have to wrap the W in some theano construct?

I think you are complicating the problem, you dont need to wrap it into a theano function, just the tensor is sufficient. It might help if you have a look at the similar posts on this discourse (search theano.scan)

I think I finally got it to work, thanks for encouragement. I haven’t evaluated the model yet, but the coefficients seem reasonable. Here’s the code:

import theano.tensor as tt
lags = 2

nlags = nn.Nlags(lags)
    
with pm.Model() as ar:
    W = pm.Normal( 'W', mu=0, sd=10, shape=( lags, 1 ) )
    
    e = tt.TensorConstant( tt.TensorType('float64', (False,True)), np.array([[0.],[0.]]))
    
    def compute_eps( price, eps_m1, eps_m2, W ):
        # eps(n) = price_train(n) - dot( [ eps(n-1), eps(n-2), ..., eps(n-lags) ], W )
        lagged = tt.stack( [eps_m1, eps_m2], axis=1 )
        return tt.unbroadcast(price - tt.dot( lagged, W )[0,...])
    
    
    eps_, _ = theano.scan( fn=compute_eps,
                          sequences=price_train,
                          outputs_info=dict(initial=e, taps=[-1,-2]),
                          non_sequences=W )
    
    eps = pm.Deterministic( 'eps', eps_ )
    
    # 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) )

The quality of sampling seems okay, but speed is very slow (~10 samples/s) compared to simple autoregressive model (~800 samples/s). I have a relevant error message that says:

C:\Users\Tomek\Anaconda3\envs\pymc3\lib\site-packages\theano\scan_module\scan_perform_ext.py:76: UserWarning: The file scan_perform.c is not available. This donot happen normally. You are probably in a strangesetup. This mean Theano can not use the cython code for scan. If youwant to remove this warning, use the Theano flag'cxx=' (set to an empty string) to disable all ccode generation.

I installed PYMC3 via anaconda ( conda install pymc3 ) and all libraries from this site are up to date or newer.
I’m working on windows 10.
Have you encountered such an error before?

2 Likes

So, I was able to fix that issue as well. It turns out that the anaconda distribution of pymc does not include the file scan_perform.c in scan_module/c_code directory. Manually copying the file from github solved the issue and the code is 4 times faster (~40 draws/s).

Too bad that I’m not able to use multiple cores while sampling on windows 10 yet. Has there been any progress on the issue?

2 Likes