A non-Gaussian autoregressive model for peat accumulation

Dear community,

I am implementing an age-depth model which which uses autoregressive gamma process.

A non-Gaussian autoregressive model is proposed for the sediment accumulation rates xj = w*xj+1 + zj. where zj ~ Gamma(alpha, beta / (1 - w)) are iid innovations. And w ~ Beta( aw, bw).

I have implemented this in pyMC as follows

w = pm.Beta('w', aw, bw)
alpha = pm.Gamma('alpha', a_alpha, (1 - w) / b_alpha, shape=K)

x = pm.Data('x', np.zeros(K))

for j in range(K - 2, -1, -1):
    x[j].set(w * x[j + 1] + (1 - w) * alpha[j])

Is this the correct way to do this?

Thank you.

Have you had a look at the time series examples in the example gallery? The custom time series model one might be of particular interest. There are also plenty of resources and discussions about similar models on this forum if you search.

Forgive my naiveness, I am very new to Bayesian statistics and PyMC. I just wanted to know that all AR models that I have seen has dependence on previous (e.g. xj = rho*xj-1) but the paper I am following does this the other way around.

So what can I do about that?

Thank you.