# Autocorrelation in a model

**URL:** <https://discourse.pymc.io/t/autocorrelation-in-a-model/1195>\
**Category:** Questions\
**Created:** [May 9, 2018, 3:55pm UTC](https://discourse.pymc.io/t/autocorrelation-in-a-model/1195 "2018-05-09T15:55:12Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![hajasw](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/hajasw/32/1525_2.png) [@hajasw](https://discourse.pymc.io/u/hajasw)\
**Post date:** [May 9, 2018, 3:55pm UTC](https://discourse.pymc.io/t/autocorrelation-in-a-model/1195/1 "2018-05-09T15:55:13Z")

</div>

Is there a practical way to implement autocorrelation in a pymc3 model?

What I want to implement is:

x\_0~normal(0,s^2)  
x\_i~normal( ro\*x\_i-1 , (1-ro^2)\*s^2 ) for i=1,2,3,…n-1

In my current application, n=106

I can’t figure out how to implement this type of model in pymc3. I am hoping there is a simple way.

Thanks,  
Wayne Hajas

---

<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:** [May 9, 2018, 7:20pm UTC](https://discourse.pymc.io/t/autocorrelation-in-a-model/1195/2 "2018-05-09T19:20:58Z")

</div>

There is AR(n) distribution for this purpose, see eg [http://docs.pymc.io/notebooks/AR.html](http://docs.pymc.io/notebooks/AR.html)

---

<div class="post-metadata">

**Author:** ![hajasw](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/hajasw/32/1525_2.png) [@hajasw](https://discourse.pymc.io/u/hajasw)\
**Post date:** [May 9, 2018, 9:25pm UTC](https://discourse.pymc.io/t/autocorrelation-in-a-model/1195/3 "2018-05-09T21:25:40Z")

</div>

Thanks Junpeng. I took a look at the AR function.

Unfortunately, the autocorrelation is just part of my model. I need to estimate all the x-values as well as the ro and sigma-values. As far as I can tell, the AR function will not help me estimate the x-values.

This model has been implemented in pymc2. I suspect the implementation is also possible in pymc3 - but I am just not seeing it.

Any assistance is appreciated.  
Wayne

---

<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:** [May 9, 2018, 9:47pm UTC](https://discourse.pymc.io/t/autocorrelation-in-a-model/1195/4 "2018-05-09T21:47:44Z")

</div>

could you post the pymc2 code?

---

<div class="post-metadata">

**Author:** ![hajasw](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/hajasw/32/1525_2.png) [@hajasw](https://discourse.pymc.io/u/hajasw)\
**Post date:** [May 10, 2018, 3:38pm UTC](https://discourse.pymc.io/t/autocorrelation-in-a-model/1195/5 "2018-05-10T15:38:35Z")

</div>

I am in the process of cleaning up the pymc2 code. I haven’t fully tested this version but this is what I have right now. I am generating auto-correlated values of LogRecruit. Thanks for taking a look at this.

```python
def Recruitment(self):
    '''Log of Relative annual Recruitment'''

    #Will be subjected to antiLogit transformation
    self.T = Normal('T', 7.78, 0.667)
    #Coefficient for autocorrelation
    self.Ro = Lambda('Ro', lambda t=self.T: expit(t))
    self.Omega = Lambda('Omega', lambda t=self.Ro: np.sqrt(1 - t * t))

    #Standard deviation of recruitment rate
    self.lnSigmaRecruit = Normal('lnSigmaRecruit', 2.5, 1)
    self.tauRecruit = Lambda(
        'tauRecruit', lambda t=self.lnSigmaRecruit: np.exp(-2 * t))

    self.LogRecruit = numpy.empty(self.nyear, dtype=object)
    self.LogRecruit[0] = Normal('LogRecruit_' + str(self.MinYear), 0,
                                self.tauRecruit)
    for y in range(1, self.nyear):
        self.LogRecruit[y] = Normal('LogRecruit_' + str(self.MinYear + y),
                                    self.Ro * self.LogRecruit[y - 1],
                                    self.tauRecruit / self.Omega / self.Omega)

```

---

<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:** [May 10, 2018, 4:35pm UTC](https://discourse.pymc.io/t/autocorrelation-in-a-model/1195/6 "2018-05-10T16:35:35Z")

</div>

What you are descripting is exactly a AR model. For example, the PyMC3 code similar to about would be something like:

```python
with pm.Model():
    T = pm.Normal('T', 7.78, tau=0.667)
    Ro = expit(T)
    Omega = tt.sqrt(1.-Ro**2)
    lnSigmaRecruit = pm.Normal('lnSigmaRecruit', 2.5, tau=1)
    tauRecruit = tt.exp(-2. * lnSigmaRecruit)
    data = pm.AR('y', 
                 Ro, 
                 tau=tauRecruit / Omega / Omega, 
                 init=pm.Normal.dist(0., tau=tauRecruit), 
                 shape=nyear)

```

---

<div class="post-metadata">

**Author:** ![hajasw](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/hajasw/32/1525_2.png) [@hajasw](https://discourse.pymc.io/u/hajasw)\
**Post date:** [May 10, 2018, 5:24pm UTC](https://discourse.pymc.io/t/autocorrelation-in-a-model/1195/7 "2018-05-10T17:24:03Z")

</div>

Maybe I am misunderstanding documentation and examples like [https://github.com/pymc-devs/pymc3/blob/master/docs/source/notebooks/AR.ipynb](https://github.com/pymc-devs/pymc3/blob/master/docs/source/notebooks/AR.ipynb)

As I understand, the code you generated will let me estimate T and lnSigmaRecruit. From those estimates, I can also estimate Ro, Omega and tauRecruit.

I also need to generate LogRecruit (from my pymc2 code). LogRecruit represents one process. I have another array of nodes representing another process. Then I multiply the two arrays together to represent the combined processes - and then this product is compared against data.

Does pm.AR generate an array of LogRecruit-nodes?

I hope that bit of explanation helps. Let me know if I am misunderstanding pm.AR

Thanks for looking at this,  
Wayne

---

<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:** [May 10, 2018, 6:30pm UTC](https://discourse.pymc.io/t/autocorrelation-in-a-model/1195/8 "2018-05-10T18:30:14Z")

</div>

Just to be clear:

```python
data = pm.AR('y', 
                 Ro, 
                 tau=tauRecruit / Omega / Omega, 
                 init=pm.Normal.dist(0., tau=tauRecruit), 
                 shape=nyear)

```

is the `LogRecruit` as in the pymc2 code.

So yes, you can generate two random arrays similarly with `pm.AR` and multiply them.

The example above is conditioned on an observed, but `pm.AR` works similar when you are treating it as a latent Variable following some distribution. See the docstring: [http://docs.pymc.io/api/distributions/timeseries.html#pymc3.distributions.timeseries.AR](http://docs.pymc.io/api/distributions/timeseries.html#pymc3.distributions.timeseries.AR)

---

<div class="post-metadata">

**Author:** ![hajasw](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/hajasw/32/1525_2.png) [@hajasw](https://discourse.pymc.io/u/hajasw)\
**Post date:** [May 10, 2018, 9:41pm UTC](https://discourse.pymc.io/t/autocorrelation-in-a-model/1195/9 "2018-05-10T21:41:24Z")

</div>

> [@hajasw](#):
>
> LogRecruit

Thanks very much Junpeng!  
pm.AR works as you said it would.  
Now all my problems are my own fault.

Cheers,  
Wayne