# MA model for time series prediction

**URL:** <https://discourse.pymc.io/t/ma-model-for-time-series-prediction/2617>\
**Category:** Questions\
**Created:** [January 31, 2019, 1:05pm UTC](https://discourse.pymc.io/t/ma-model-for-time-series-prediction/2617 "2019-01-31T13:05:20Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Tomasz\_Gac](https://avatars.discourse-cdn.com/v4/letter/t/439d5e/32.png) [@Tomasz\_Gac](https://discourse.pymc.io/u/Tomasz_Gac)\
**Post date:** [January 31, 2019, 1:05pm UTC](https://discourse.pymc.io/t/ma-model-for-time-series-prediction/2617/1 "2019-01-31T13:05:20Z")

</div>

I’m trying to implement a moving average model based on [this](https://github.com/stan-dev/example-models/blob/master/misc/moving-avg/ma2.stan) 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.

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [January 31, 2019, 4:06pm UTC](https://discourse.pymc.io/t/ma-model-for-time-series-prediction/2617/2 "2019-01-31T16:06:40Z")

</div>

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.)

---

<div class="post-metadata">

**Author:** ![Tomasz\_Gac](https://avatars.discourse-cdn.com/v4/letter/t/439d5e/32.png) [@Tomasz\_Gac](https://discourse.pymc.io/u/Tomasz_Gac)\
**Post date:** [January 31, 2019, 4:25pm UTC](https://discourse.pymc.io/t/ma-model-for-time-series-prediction/2617/3 "2019-01-31T16:25:05Z")

</div>

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`?

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [January 31, 2019, 9:14pm UTC](https://discourse.pymc.io/t/ma-model-for-time-series-prediction/2617/4 "2019-01-31T21:14:53Z")

</div>

_much_ slower. See here for some more information [https://docs.pymc.io/notebooks/PyMC3\_tips\_and\_heuristic.html`](https://docs.pymc.io/notebooks/PyMC3_tips_and_heuristic.html%60)

---

<div class="post-metadata">

**Author:** ![Tomasz\_Gac](https://avatars.discourse-cdn.com/v4/letter/t/439d5e/32.png) [@Tomasz\_Gac](https://discourse.pymc.io/u/Tomasz_Gac)\
**Post date:** [February 1, 2019, 1:18pm UTC](https://discourse.pymc.io/t/ma-model-for-time-series-prediction/2617/5 "2019-02-01T13:18:30Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [February 1, 2019, 2:56pm UTC](https://discourse.pymc.io/t/ma-model-for-time-series-prediction/2617/6 "2019-02-01T14:56:02Z")

</div>

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`)

---

<div class="post-metadata">

**Author:** ![Tomasz\_Gac](https://avatars.discourse-cdn.com/v4/letter/t/439d5e/32.png) [@Tomasz\_Gac](https://discourse.pymc.io/u/Tomasz_Gac)\
**Post date:** [February 1, 2019, 7:58pm UTC](https://discourse.pymc.io/t/ma-model-for-time-series-prediction/2617/7 "2019-02-01T19:58:52Z")

</div>

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](http://deeplearning.net/software/theano/install_windows.html) are up to date or newer.  
I’m working on windows 10.  
Have you encountered such an error before?

---

<div class="post-metadata">

**Author:** ![Tomasz\_Gac](https://avatars.discourse-cdn.com/v4/letter/t/439d5e/32.png) [@Tomasz\_Gac](https://discourse.pymc.io/u/Tomasz_Gac)\
**Post date:** [February 1, 2019, 8:34pm UTC](https://discourse.pymc.io/t/ma-model-for-time-series-prediction/2617/8 "2019-02-01T20:34:50Z")

</div>

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](https://github.com/Theano/Theano/blob/master/theano/scan_module/c_code/scan_perform.c) 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?