Simple linear timeseries model failed to sample

Questions version agnostic
1

I am trying to do this time series model where the error measured at time t depends on the error at time t-1:

$$ Y_{t} = a + b*cos(\pi/12(t-c)) + d * activity_t + \epsilon_t$$

$$ \epsilon_t \sim \mathcal{N} (\theta* \epsilon_{t-1}, \sigma^2)$$

I couldnt get my model to sample correctly. Could you please help me figure out where I got it wrong?

with pm.Model() as ar1_cos_steps:
    # priors 
    a = pm.Normal('a', mu=0, sigma = 10.0) #Hr_base
    b = pm.Normal('b', mu=1, sigma = 1.0)  #amplitude
    c = pm.Normal('c', mu=1, sigma = 1.0)  #phase param
    d = pm.Normal('d', mu=0, sigma = 10.0) #activity coefficient
    # sigma = pm.HalfNormal("sigma", sigma=1)
    
    theta = pm.Normal("theta", 0.0, 1.0) #error correlation terms 
    tau = pm.Exponential("tau", 0.5) # precision of the innovation term tau=1/sigma
    
    # Gaussian noise
    err_center = pm.Normal("center", mu=0.0, sigma=1.0)
    
    # Expected value of outcome
    mu_y = a + b * pm.math.cos(np.pi*(t-c)/(12*60/bin_size)) + d*activity
    
    # noise term 
    err = y - mu_y 
    
    likelihood = pm.AR1("err", k=theta, tau_e=tau, observed= err - err_center)

    trace_ar1 = pm.sample(3000, tune=2000, step = pm.Metropolis()) # metropolis sampler is a lot easier to work with
    idata_ar1 = az.from_pymc3(trace_ar1)

I am coding up this model as an exercise from a Cell Reports paper that I found. Could you please help me understand why the authors wanted the errors to correlate (as they provided no explanation in their method section)? To me this is counterintuitive as why one would want their model to always underestimate or overestimate. Please forgive my little understanding of the topic. Thanks so much!

2

Correlated errors like this are called a Moving Average (MA) process. It’s not that the authors “want” the errors to correlate. They are positing a particular stochastic structure to the time series being studied.

Consider quarterly oil well production. An i.i.d. normally distributed error might be a freak well fire. As a result, the well’s production at time t will be reduced. In addition, at t+1, clean-up and repair might not be done, so we might expect additional disturbance, proportional to the amount of damage done by the original fire shock (\theta \in (0, 1)). After those two quarters, though, the effect of the fire is totally gone.

In contrast, if we modeled the oil well production as an autoregressive process, like y_t = \rho y_{t-1} + \epsilon_t, the fire shock will last far into the future. To see this, repeatedly substitute the AR process into itself. After \tau quarters, the effect of the fire will be \rho^\tau \epsilon_t, which might stay quite large for quite some time, depending on the magnitude of \rho.

As for the model, what specific errors are you getting? I strongly recommend you use NUTS instead of Metropolis. Metropolis is “easier to work with” in the sense that it complains less, but when NUTS complains it is telling you something is wrong with your model. Metropolis, on the other hand, is happy to silently return incorrect results (biased/incomplete posterior exploration).

1 Like