Skip to content

Conversation

@jessegrabowski
Copy link
Member

The purpose of the duration argument is to allow high-frequency seasonal patterns inside of low-frequency data. For example, if we have hourly observations, we might want to define a seasonal pattern over days of the week.

This is currently broken -- it simply doesn't do anything. My first cut was to get it working. To do this I made a single huge transition matrix that rotates through duration * season_length states, which are all copies of eachother. This works but was unusable.

To fix this I made the component use a time-varying transition matrix. This required some refactor of the PyMCStateSpace to look for a special n_timesteps variable. This is replaced with the length of the data or forecasts at graph building time. It should also unlock some other features, like time-varying intercept components.

Here is an example of the new component, estimating day-of-week effect on hourly data:

image

Copy link

Copilot AI left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pull request overview

This PR fixes TimeSeasonality(duration=...) by switching to time-varying state-space matrices driven by an n_timesteps placeholder that is substituted at graph-build time (fit/forecast/prior paths). It also expands test coverage around time-varying behavior and seasonal state ordering.

Changes:

  • Implement time-varying transitions for TimeSeasonality (with use_time_varying and start_state) and introduce an n_timesteps symbolic placeholder.
  • Refactor PyMCStateSpace graph building/simulation paths to replace n_timesteps with data/forecast lengths.
  • Update and add tests for time-varying matrices, start-state shifting, and multivariate seasonality composition.

Reviewed changes

Copilot reviewed 7 out of 7 changed files in this pull request and generated 8 comments.

Show a summary per file
File Description
pymc_extras/statespace/models/structural/components/seasonality.py Implements time-varying transition approach for duration>1, adds start_state and use_time_varying.
pymc_extras/statespace/core/statespace.py Adds n_timesteps placeholder handling and replaces it with data/forecast lengths in multiple build paths.
pymc_extras/statespace/models/structural/core.py Adds n_timesteps placeholder to structural components.
tests/statespace/test_utilities.py Updates symbolic-matrix unpacking to substitute n_timesteps for tests/simulation utilities.
tests/statespace/core/test_statespace.py Adds a minimal time-varying intercept model and exercises key sampling/forecast/IRF paths.
tests/statespace/models/structural/test_against_statsmodels.py Adjusts seasonal initial-state ordering to match internal vs statsmodels conventions.
tests/statespace/models/structural/components/test_seasonality.py Adds/updates tests for TimeSeasonality duration/start-state/time-varying mode and related compositions.
Comments suppressed due to low confidence (2)

pymc_extras/statespace/models/structural/core.py:549

            self.make_symbolic_graph()

pymc_extras/statespace/core/statespace.py:303

        self.make_symbolic_graph()

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

k_posdef = self._k_posdef_per_endog()

R = pt.zeros((k_states, k_posdef))[0, 0].set(1.0)
self.ssm["selection", :, :] = pt.join(0, *[R for _ in range(k_endog_effective)])
Copy link

Copilot AI Feb 8, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For multivariate (k_endog_effective > 1) this builds the selection matrix by vertically concatenating identical R blocks, which will broadcast into the full (k_states, k_posdef) shape and effectively makes each endog’s shock drive all endogs’ state blocks. Selection should be block-diagonal (one shock per endog affecting only its own seasonal states). Consider using pt.linalg.block_diag(*[R for _ in range(k_endog_effective)]) (and optionally pt.specify_shape) so R has shape (k_states_total, k_posdef_total) without cross-series coupling.

Suggested change
self.ssm["selection", :, :] = pt.join(0, *[R for _ in range(k_endog_effective)])
if k_endog_effective == 1:
self.ssm["selection", :, :] = R
else:
R_blocks = [R for _ in range(k_endog_effective)]
self.ssm["selection", :, :] = pt.linalg.block_diag(*R_blocks)

Copilot uses AI. Check for mistakes.
Comment on lines +412 to +413
decimal = s_str.split(".") if "." in (s_str := str(s)) else "0"
T = int(s * 10 ** len(decimal))
Copy link

Copilot AI Feb 8, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This computation of T is incorrect for integer season lengths: when s is an int (e.g., 5), decimal becomes the string "0" so len(decimal)==1 and T becomes 50, causing the periodicity assertion to be wrong. Compute the scale factor using the length of the decimal part only (e.g., decimal_part = s_str.split(".")[1] if "." in s_str else ""; T = int(s * 10**len(decimal_part))).

Suggested change
decimal = s_str.split(".") if "." in (s_str := str(s)) else "0"
T = int(s * 10 ** len(decimal))
s_str = str(s)
decimal_part = s_str.split(".")[1] if "." in s_str else ""
T = int(s * 10 ** len(decimal_part))

Copilot uses AI. Check for mistakes.
Comment on lines +196 to +205
x0 = np.zeros((mod.k_endog, mod.n_seasons), dtype=config.floatX)
x0[0, 0], x0[1, 0] = 1.0, 2.0
params = {"params_season": x0, "sigma_season": np.array([0.0, 0.0], dtype=config.floatX)}

for matrix, expected in zip(
[x0, T, Z, R, Q],
[expected_x0, expected_T, expected_Z, expected_R, expected_Q],
):
np.testing.assert_allclose(matrix, expected)
steps = (
s * d * 10
) # assert_pattern_repeats requires this to be divisible by pattern period s * d
x, y = simulate_from_numpy_model(mod, rng, params, steps=steps)
assert_pattern_repeats(y[:, 0], s * d, atol=ATOL, rtol=RTOL)
assert_pattern_repeats(y[:, 1], s * d, atol=ATOL, rtol=RTOL)
Copy link

Copilot AI Feb 8, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This test sets sigma_season to zeros, so it doesn’t exercise the innovations path (R/Q usage) for the time-varying multivariate case. That can hide issues where the selection matrix is shaped/built incorrectly (e.g., shocks bleeding across endog). Consider using non-zero sigma_season and asserting basic structure (e.g., R is block-diagonal / output variance differs per series) so time-varying + multivariate innovations are actually covered.

Copilot uses AI. Check for mistakes.
assert len(ss_mod.coords["state_quarterly"]) == s - 1

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"])
Copy link

Copilot AI Feb 8, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Variable _P0 is not used.

Suggested change
_P0 = pm.Deterministic("P0", pt.eye(ss_mod.k_states) * 10, dims=ss_mod.param_dims["P0"])
pm.Deterministic("P0", pt.eye(ss_mod.k_states) * 10, dims=ss_mod.param_dims["P0"])

Copilot uses AI. Check for mistakes.

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"])
_params = pm.Normal(
Copy link

Copilot AI Feb 8, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Variable _params is not used.

Suggested change
_params = pm.Normal(
pm.Normal(

Copilot uses AI. Check for mistakes.
y = rng.normal(size=(n_obs, 1)).astype(floatX)

with pm.Model(coords=ss_mod_time_varying.coords) as m:
slope = pm.Normal("slope", mu=0, sigma=1)
Copy link

Copilot AI Feb 8, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Variable slope is not used.

Suggested change
slope = pm.Normal("slope", mu=0, sigma=1)
slope = pm.Normal("slope", mu=0, sigma=1)
slope_det = pm.Deterministic("slope_det", slope)

Copilot uses AI. Check for mistakes.

with pm.Model(coords=ss_mod_time_varying.coords) as m:
slope = pm.Normal("slope", mu=0, sigma=1)
P0 = pm.Deterministic(
Copy link

Copilot AI Feb 8, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Variable P0 is not used.

Suggested change
P0 = pm.Deterministic(
pm.Deterministic(

Copilot uses AI. Check for mistakes.
P0 = pm.Deterministic(
"P0", pt.eye(ss_mod_time_varying.k_states) * 1.0, dims=["state", "state_aux"]
)
x0 = pm.Normal("x0", dims=["state"])
Copy link

Copilot AI Feb 8, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Variable x0 is not used.

Suggested change
x0 = pm.Normal("x0", dims=["state"])
x0 = pm.Normal("x0", dims=["state"])
x0_det = pm.Deterministic("x0_det", x0, dims=["state"])

Copilot uses AI. Check for mistakes.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

bug Something isn't working enhancements New feature or request statespace

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant