# Autocorrelation in a model

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

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

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

could you post the pymc2 code?

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.

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

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

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

Maybe I am misunderstanding documentation and examples like 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

Just to be clear:

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

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

Cheers,
Wayne