TimeSeasonality#

class pymc_experimental.statespace.models.structural.TimeSeasonality(season_length: int, innovations: bool = True, name: str | None = None, state_names: list | None = None, pop_state: bool = True)[source]#

Seasonal component, modeled in the time domain

Parameters:

  • season_length (int) – The number of periods in a single seasonal cycle, e.g. 12 for monthly data with annual seasonal pattern, 7 for daily data with weekly seasonal pattern, etc.

  • innovations (bool, default True) – Whether to include stochastic innovations in the strength of the seasonal effect

  • name (str, default None) – A name for this seasonal component. Used to label dimensions and coordinates. Useful when multiple seasonal components are included in the same model. Default is f"Seasonal[s={season_length}]"

  • state_names (list of str, default None) –

    List of strings for seasonal effect labels. If provided, it must be of length season_length. An example would be state_names = ['Mon', 'Tue', 'Wed', 'Thur', 'Fri', 'Sat', 'Sun'] when data is daily with a weekly seasonal pattern (season_length = 7).

    If None, states will be numbered [State_0, ..., State_s]

Notes

A seasonal effect is any pattern that repeats every fixed interval. Although there are many possible ways to model seasonal effects, the implementation used here is the one described by [1] as the “canonical” time domain representation. The seasonal component can be expressed:

\[\gamma_t = -\sum_{i=1}^{s-1} \gamma_{t-i} + \omega_t, \quad \omega_t \sim N(0, \sigma_\gamma)\]

Where \(s\) is the seasonal_length parameter and \(\omega_t\) is the (optional) stochastic innovation. To give interpretation to the \(\gamma\) terms, it is helpful to work through the algebra for a simple example. Let \(s=4\), and omit the shock term. Define initial conditions \(\gamma_0, \gamma_{-1}, \gamma_{-2}\). The value of the seasonal component for the first 5 timesteps will be:

\[\begin{split}\begin{align} \gamma_1 &= -\gamma_0 - \gamma_{-1} - \gamma_{-2} \\ \gamma_2 &= -\gamma_1 - \gamma_0 - \gamma_{-1} \\ &= -(-\gamma_0 - \gamma_{-1} - \gamma_{-2}) - \gamma_0 - \gamma_{-1} \\ &= (\gamma_0 - \gamma_0 )+ (\gamma_{-1} - \gamma_{-1}) + \gamma_{-2} \\ &= \gamma_{-2} \\ \gamma_3 &= -\gamma_2 - \gamma_1 - \gamma_0 \\ &= -\gamma_{-2} - (-\gamma_0 - \gamma_{-1} - \gamma_{-2}) - \gamma_0 \\ &= (\gamma_{-2} - \gamma_{-2}) + \gamma_{-1} + (\gamma_0 - \gamma_0) \\ &= \gamma_{-1} \\ \gamma_4 &= -\gamma_3 - \gamma_2 - \gamma_1 \\ &= -\gamma_{-1} - \gamma_{-2} -(-\gamma_0 - \gamma_{-1} - \gamma_{-2}) \\ &= (\gamma_{-2} - \gamma_{-2}) + (\gamma_{-1} - \gamma_{-1}) + \gamma_0 \\ &= \gamma_0 \\ \gamma_5 &= -\gamma_4 - \gamma_3 - \gamma_2 \\ &= -\gamma_0 - \gamma_{-1} - \gamma_{-2} \\ &= \gamma_1 \end{align}\end{split}\]

This exercise shows that, given a list initial_conditions of length s-1, the effects of this model will be:

  • Period 1: -sum(initial_conditions)

  • Period 2: initial_conditions[-1]

  • Period 3: initial_conditions[-2]

  • …

  • Period s: initial_conditions[0]

  • Period s+1: -sum(initial_condition)

And so on. So for interpretation, the season_length - 1 initial states are, when reversed, the coefficients associated with state_names[1:].

Warning

Although the state_names argument expects a list of length season_length, only state_names[1:] will be saved as model dimensions, since the 1st coefficient is not identified (it is defined as \(-\sum_{i=1}^{s} \gamma_{t-i}\)).

Examples

Estimate monthly with a model with a gaussian random walk trend and monthly seasonality:

from pymc_experimental.statespace import structural as st
import pymc as pm
import pytensor.tensor as pt
import pandas as pd

# Get month names
state_names = pd.date_range('1900-01-01', '1900-12-31', freq='MS').month_name().tolist()

# Build the structural model
grw = st.LevelTrendComponent(order=1, innovations_order=1)
annual_season = st.TimeSeasonality(season_length=12, name='annual', state_names=state_names, innovations=False)
ss_mod = (grw + annual_season).build()

# Estimate with PyMC
with pm.Model(coords=ss_mod.coords) as model:
    P0 = pm.Deterministic('P0', pt.eye(ss_mod.k_states) * 10, dims=ss_mod.param_dims['P0'])
    intitial_trend = pm.Deterministic('initial_trend', pt.zeros(1), dims=ss_mod.param_dims['initial_trend'])
    annual_coefs = pm.Normal('annual_coefs', sigma=1e-2, dims=ss_mod.param_dims['annual_coefs'])
    trend_sigmas = pm.HalfNormal('trend_sigmas', sigma=1e-6, dims=ss_mod.param_dims['trend_sigmas'])
    ss_mod.build_statespace_graph(data, mode='JAX')
    idata = pm.sample(nuts_sampler='numpyro')

References

__init__(season_length: int, innovations: bool = True, name: str | None = None, state_names: list | None = None, pop_state: bool = True)[source]#

Methods

__init__(season_length[, innovations, name, ...])

build([name, filter_type, verbose])

Build a StructuralTimeSeries statespace model from the current component(s)

make_and_register_data(name, shape[, dtype])

Helper function to create a pytensor symbolic variable and register it in the _name_to_data dictionary

make_and_register_variable(name, shape[, dtype])

Helper function to create a pytensor symbolic variable and register it in the _name_to_variable dictionary

make_symbolic_graph()

populate_component_properties()