Time Series: State Space Model with stochastic level and stochastic seasonal components

Questions
1

Hi,

I’m new to PyMC3 and I’m implementing different state space models - in this case a local level model with stochastic seasonal component. Something like this:

\begin{align*} y_t &= \mu_t + \gamma_t + \epsilon_t, & \epsilon_t \sim NID(0, \sigma^2_\epsilon) \\ \mu_{t+1} &= \mu_t + \xi_t, & \xi_t \sim NID(0, \sigma^2_\xi) \\ \gamma_{1,t+1} &= -\gamma_{1,t} - \gamma_{2,t} - ... - \gamma_{12,t} + w_t, & w_t \sim NID(0, \sigma^2_w) \\ \gamma_{2,t+1} &= \gamma_{1,t} \\ \gamma_{3,t+1} &= \gamma_{2,t} \\ [...] \end{align*}

As an example we could use the data found here (use the log of the observations). It is an example from State Space Time Series Analysis. This model should yield good statistical tests both in terms of independence and normality of the residuals as it is a good fit for the data.

This was my attempt.


w = T.dvector("w")

γ = T.dvector("γ")

def fn_γ_t(w_t, γ_pre1, γ_pre2, γ_pre3, γ_pre4, γ_pre5, γ_pre6, γ_pre7, γ_pre8, γ_pre9, γ_pre10, γ_pre11, γ_pre12):
    γ = -γ_pre1 - γ_pre2 - γ_pre3 - γ_pre4 - γ_pre5 - γ_pre6 - γ_pre7 - γ_pre8 - γ_pre9 - γ_pre10 - γ_pre11 - γ_pre12 + w_t
    return γ

with pm.Model() as model:
    σ_ϵ = pm.HalfCauchy('σ_ϵ', 1)
    σ_ξ = pm.HalfCauchy('σ_ξ', 1)
    σ_w = pm.HalfCauchy('σ_w', 1)
    
    w = pm.Normal('w', mu=0, sigma=σ_w, shape=len(y))
    
    i = theano.shared(np.zeros(12))
    γ_, updates = theano.scan(fn=fn_γ_t, 
                                sequences = dict(input=w), 
                                outputs_info = dict(initial = i, 
                                                    taps=[-1, -2, -3, -4, -5, -6, -7, -8, -9, -10, -11, -12]), 
                                n_steps=len(y))
    γ = pm.Deterministic('γ', γ_ )

    μ = pm.GaussianRandomWalk('μ', mu=0, sd=σ_ξ, shape=len(y))

    y_pred = pm.Normal('y_pred', mu=μ + γ, sd=σ_ϵ, observed=y)

    trace = pm.sample(2000, tune=2000)

I’m clearly not defining it well. What would be the best/correct way? This is just a simple component of a more complex model, I also would need it to be as performant as possible. Feel free to use this example or, if it is simpler, simulated data.

2

Hi and welcome

I didn’t quite follow your seasonality definition. I don’t have much experience with time series but here is a model with random walk and additive monthly adjustments.

with pm.Model() as model3:
    σ_ϵ = pm.HalfCauchy('σ_ϵ', 1)
    σ_ξ = pm.HalfCauchy('σ_ξ', 1)
    
    μ = pm.GaussianRandomWalk('μ', mu=0, sd=σ_ξ, init=pm.Normal.dist(0, 5), shape=len(y))
    
    alpha_pre = pm.Normal('alpha_pre', 0, 1, shape=12)
    alpha = pm.Deterministic('alpha', alpha_pre - alpha_pre.mean())
    seasonality = pm.Deterministic('seasonality', alpha[np.arange(len(y))%12])
    
    loc = pm.Deterministic('loc', seasonality + μ)
    obs = pm.Normal('obs', mu=loc, sigma=σ_ϵ, observed=y)
    idata3 = pm.sample(return_inferencedata=True)
1 Like
3

That actually defines a deterministic seasonal component, but it is easy to change to a stochastic one. Thanks, it was very helpful!

4

Great, happy to help.

Are you able to post your final model for other’s benefit?

2 Likes