diff --git a/CONTEXT.md b/CONTEXT.md index f62ffad..bf7e53f 100644 --- a/CONTEXT.md +++ b/CONTEXT.md @@ -58,6 +58,10 @@ The dynamic response of each variable to a unit structural shock at horizons `0. The response of each endogenous variable to a unit impulse in an *exogenous* regressor at horizons `0..h`: `Psi_h = Phi_h @ B_exog`. Shares the moving-average coefficients `Phi_h` with the IRF, but needs no identification scheme — exogenous regressors are exogenous by assumption — so it lives on `FittedVAR`, not `IdentifiedVAR`, and takes no `at`. The cumulative variant is the response to a permanent unit *step* rather than a one-off impulse. All dynamics come from the endogenous lag structure; `B_exog` enters contemporaneously and carries no lags of its own. _Avoid_: "exogenous IRF" — IRF is reserved for responses to identified structural shocks. +**Deterministic design (`DeterministicDesign`)**: +A composed set of *deterministic terms* — `Trend`, `Fourier`, `SeasonalDummies`, `BreakDummy` — that generate an exogenous block from the index alone: total functions of a timestamp, never of the data. Elapsed time is counted in integer pandas period ordinals from the first timestamp of the estimation index, which buys the **continuation property**: `design.build(index[:T+h]).iloc[T:] == design.extend(index[:T], h)`. That is what makes `exog_future` safe — it extends the design and reorders columns *by name* to match `exog_names`, closing the silent-wrong-answer hole in positional `exog_future` handling. `build` rank-checks against the intercept every estimator fits; `extend` does not (short horizons are legitimately deficient). The `DeterministicTerm` Protocol is the extension point. +_Avoid_: "exogenous data" for these terms — they are generated, not observed. Real covariates (solar forcing, ENSO, CO2) are *data* and belong to the datasets spine, not here. + **Historical decomposition (HD)**: The attribution of each in-sample observation to a deterministic baseline (initial conditions, intercept, and any exogenous path) plus the *propagated* contribution of each structural shock: `c_{j,t} = P_t[:,j] ε_{j,t} + Σ_i A_i c_{j,t-i}`, with `y_t = baseline_t + Σ_j c_{j,t}` holding exactly. "Propagated" is load-bearing: a shock's contribution carries forward through the lag dynamics beyond its impact period. _Avoid_: presenting the contemporaneous split of the one-step forecast error (`u_t = Σ_j P[:,j] ε_{j,t}`) as HD — that is a residual decomposition, and was the (incorrect) behaviour of the implementation prior to the scenario-analysis stack (2026-07). diff --git a/docs/how-to/data-preparation.md b/docs/how-to/data-preparation.md index 2f1047b..2903fd7 100644 --- a/docs/how-to/data-preparation.md +++ b/docs/how-to/data-preparation.md @@ -26,6 +26,19 @@ data = VARData.from_df( ) ``` +Trends, seasonal cycles and dated breaks do not need to come from a file — +build them from the index itself with a +[deterministic design](deterministic-regressors.md), which also produces the +matching future block for forecasting: + +```python +from impulso import DeterministicDesign, Fourier, Trend + +design = DeterministicDesign(terms=[Trend(degree=1, scale=120.0), Fourier(period=12, order=2)]) +frame = pd.concat([df, design.build(df.index)], axis=1) +data = VARData.from_df(frame, endog=["gdp", "inflation", "rate"], exog=design.column_names) +``` + ## From NumPy arrays If you already have arrays, pass them directly: diff --git a/docs/how-to/deterministic-regressors.md b/docs/how-to/deterministic-regressors.md new file mode 100644 index 0000000..6b8bb87 --- /dev/null +++ b/docs/how-to/deterministic-regressors.md @@ -0,0 +1,297 @@ +# Deterministic Regressors for Climate VARs + +Climate series arrive with structure that is not dynamics. A monthly +temperature record for Berlin swings by twenty degrees between January and +July, drifts upward across four decades, and jumps the day the station moved. +None of that is the lag structure a VAR exists to estimate, and a VAR fitted on +the raw series spends its coefficients re-learning the calendar. + +This guide shows how to build that structure as *deterministic regressors* — +total functions of a timestamp, with no data in them — and feed them to the +model through the exogenous block. + +## Two routes, and when to take each + +**Transform it away.** Subtract a month-of-year climatology, standardise, and +fit the VAR on anomalies: + +```python +climatology = raw.groupby(raw.index.month).transform("mean") +anomalies = raw - climatology +anomalies = (anomalies - anomalies.mean()) / anomalies.std() +``` + +This is what the [conjugate VAR tutorial](../tutorials/conjugate-var.py) does. +It works with both estimators, keeps the model small, and gives impulse +responses that read in standard deviations. Its cost is that the seasonal +adjustment is a point estimate: its uncertainty never reaches the posterior, +and you cannot ask how large the annual cycle is or whether it changed. + +**Model it.** Put the trend, cycle and breaks in the exogenous block. The +coefficients get posteriors you can read and plot, the seasonal uncertainty +propagates into the forecast bands, and a level-shift coefficient becomes a +dynamic multiplier — the full propagated path of the regime change, not just +its impact. The cost is that only the NUTS estimator (`impulso.VAR`) consumes +exogenous regressors, and each column is a parameter per variable. + +The rest of this page is the second route. + +## Building a design + +The four term types live in `impulso.deterministic`. + +A **trend** counts periods elapsed since the start of the sample: + +```python +from impulso import Trend + +Trend(degree=1, scale=120.0) # linear, in units of decades on monthly data +Trend(degree=2, scale=120.0) # adds a squared term +``` + +`scale` is not cosmetic, though it is no longer a prior-fighting problem. The +prior on the exogenous coefficients adapts to each regressor's sample spread +(see `VAR.exog_prior_scale`), so an unscaled 540-month trend reaching 539 is +not crushed by a prior fixed in coefficient space. What `scale` still buys you +is the coefficient's units and the sampler's geometry: divide by a fixed, +interpretable constant — periods per decade — and the coefficient reads as +"change per decade" while the design column stays O(1) instead of spanning +three orders of magnitude. + +:::{admonition} Keep `scale` independent of the sample +:class: warning +`scale` must not depend on the sample length. `scale=len(index)` breaks the +continuation property below: the design you extend with is no longer the design +you fitted. +::: + +**Fourier harmonics** represent a smooth cycle of known length at two +coefficients per harmonic pair: + +```python +from impulso import Fourier + +Fourier(period=12, order=2) # annual cycle on monthly data, 4 columns +Fourier(period=4, order=1) # annual cycle on quarterly data, 2 columns +``` + +The cycle length is always explicit — nothing here infers a period from the +data. At most `period / 2` harmonic pairs are identified; more is rejected at +construction. + +**Seasonal dummies** are the non-parametric alternative — a free coefficient +per calendar unit: + +```python +from impulso import SeasonalDummies + +SeasonalDummies(season="month") # 11 columns (January dropped) +SeasonalDummies(season="quarter") # 3 columns +SeasonalDummies(season="dayofweek") # 6 columns (levels 0-6, Monday = 0) +``` + +One level is always dropped, because every Impulso estimator fits an +unconditional intercept and a full set of indicators sums to it. Use +`reference=` to choose which. + +**Break dummies** mark a known, dated discontinuity: + +```python +from impulso import BreakDummy + +BreakDummy(date="1991-06-15") # level shift from that date on +BreakDummy(date="1991-06-15", kind="pulse") # that period only +``` + +Compose them into a `DeterministicDesign`: + +```python +from impulso import BreakDummy, DeterministicDesign, Fourier, Trend + +design = DeterministicDesign( + terms=[ + Trend(degree=1, scale=120.0), + Fourier(period=12, order=2), + BreakDummy(date="1991-06-15"), + ], + freq="MS", +) +frame = design.build(anomalies.index) +``` + +## The column-name contract + +`design.column_names` is knowable before any index exists, and it is exactly +what `build` and `extend` emit, in that order: + +| Term | Columns | +| --- | --- | +| `Trend(degree=3)` | `trend`, `trend_squared`, `trend_cubed` | +| `Fourier(period=12, order=2)` | `sin(1,12)`, `cos(1,12)`, `sin(2,12)`, `cos(2,12)` | +| `SeasonalDummies(season="month")` | `month_2` … `month_12` | +| `SeasonalDummies(season="quarter")` | `quarter_2`, `quarter_3`, `quarter_4` | +| `SeasonalDummies(season="dayofweek")` | `dow_1` … `dow_6` | +| `BreakDummy(date="1991-06-15")` | `level_1991-06-15` | +| `BreakDummy(date="1991-06-15", kind="pulse")` | `pulse_1991-06-15` | + +Those names travel: they become `exog_names` on `VARData`, and PyMC labels the +`exog` coordinate of `B_exog` with them, so the posterior comes back +self-describing. + +## Composing into VARData + +There is no special constructor — concatenate and hand the column names over: + +```python +import pandas as pd +from impulso import VARData + +frame = pd.concat([anomalies, design.build(anomalies.index)], axis=1) +data = VARData.from_df( + frame, + endog=list(anomalies.columns), + exog=design.column_names, +) +``` + +### The no-NaN invariant + +Deterministic terms are total functions of a timestamp: they cannot produce a +missing value. So if `VARData` rejects your frame with + +``` +ValueError: exog contains NaN or Inf values +``` + +the design and the endogenous block are misaligned — `pd.concat` filled the +holes. The fix is the order of operations: + +1. Transform the endogenous data (differences, logs, anomalies). +2. `dropna()`. +3. Build the design **on the index that survived**. +4. Concatenate. + +Never `dropna()` after the concatenation: that silently drops observations to +paper over a misalignment. + +Note also that `VAR.fit` discards the first `p` rows of both blocks to form the +lag matrices, so a design column that is only non-zero in the first few periods +never reaches the model. + +## Fitting and forecasting + +The payoff is that one design object serves estimation and forecasting. Because +elapsed time is anchored to the calendar rather than to row position, the block +`extend` writes for the future is exactly the block `build` would have written +had the sample run longer — the *continuation property*: + +```python +design.build(index[: T + h]).iloc[T:] == design.extend(index[:T], h) +``` + +`exog_future` wraps that up and, crucially, reorders the columns **by name** to +match the fitted `exog_names`. `forecast` indexes the block positionally, so a +permuted design would otherwise be a silently wrong forecast rather than an +error. + +The following recipe is exercised end to end in +`tests/test_deterministic.py::test_deterministic_design_end_to_end`: + +```python +import numpy as np +import pandas as pd + +from impulso import VAR, DeterministicDesign, Fourier, NUTSSampler, Trend, VARData + +# A short monthly two-variable sample standing in for climate anomalies. +rng = np.random.default_rng(7) +index = pd.date_range("2000-01-01", periods=120, freq="MS") +endog = pd.DataFrame( + rng.standard_normal((len(index), 2)).cumsum(axis=0) * 0.1, + index=index, + columns=["temperature", "precipitation"], +) + +# One design, used for estimation and for forecasting. +design = DeterministicDesign( + terms=[Trend(degree=1, scale=120.0), Fourier(period=12, order=1)], + freq="MS", +) + +frame = pd.concat([endog, design.build(index)], axis=1) +data = VARData.from_df(frame, endog=list(endog.columns), exog=design.column_names) + +fitted = VAR(lags=1).fit(data, sampler=NUTSSampler(draws=50, tune=50, chains=2, cores=1, random_seed=42)) + +# The posterior labels B_exog with the design's own column names. +assert list(fitted.idata.posterior["B_exog"].coords["exog"].values) == design.column_names + +forecast = fitted.forecast(steps=12, exog_future=design.exog_future(fitted, 12)) + +assert forecast.median().shape == (12, 2) +``` + +The same array feeds `conditional_forecast` and `structural_scenario`. + +## The estimator boundary + +:::{admonition} The conjugate estimator refuses exogenous regressors +:class: warning +`ConjugateVAR` does not consume exogenous regressors. Handing it a `VARData` +that carries a deterministic design raises: + +> ConjugateVAR does not support exogenous regressors: the conjugate engine +> estimates endogenous dynamics only, and silently ignoring the exog block +> would corrupt downstream forecasts. Drop exog from VARData or use the +> PyMC/NUTS estimator (impulso.VAR), which consumes it. + +Two ways forward. Fit with `impulso.VAR` and keep the design in the model — the +coefficients get posteriors. Or residualise first: regress each endogenous +series on `design.build(index)` by OLS and fit the conjugate VAR on the +residuals, accepting that the deterministic part's uncertainty is discarded +rather than propagated. +::: + +## Collinearity rules + +`build` checks the design's rank against the intercept every estimator fits and +refuses a deficient one, because its coefficients are not identified. The +common causes, all named in the error message: + +- `SeasonalDummies(drop_first=False)` — the levels sum to the intercept. +- `Fourier(period=12, order=6)` — at exactly the Nyquist limit the top sine is + identically zero at every sampled point. +- Seasonal dummies **and** harmonics of the same cycle (`season="month"` with + `Fourier(period=12, ...)`) — the harmonics live in the span of the level + indicators. Keep one or the other. +- A level break at or before the first observation (constant in-sample) or + after the last (never occurs); a pulse break on a date the index does not + contain (zero everywhere). +- Fewer observations than columns. + +`extend` deliberately does **not** rank-check: a three-step forecast block +cannot span twelve month dummies, and nothing is being estimated from those +rows. + +## Reading the coefficients + +`B_exog` is labelled by design column, so posterior summaries name themselves. +For a level break the more interesting object is the dynamic multiplier: +`fitted.dynamic_multiplier(cumulative=True)` propagates a permanent unit step +through the lag dynamics, which for a `level_*` column *is* the full adjustment +path of the regime shift — not merely its impact-period effect. See +[`FittedVAR`](../reference/fitted.md) for the method and +[the result objects](../reference/results.md) for what it returns. + +## Scope + +This module does calendar arithmetic and nothing else. Holiday and business +calendars, interaction terms, slope breaks and inferring cycle lengths from the +data are all out of scope. So is anything that is *data* rather than +arithmetic: solar forcing, ENSO indices and CO2 concentrations are covariates +to be loaded from a dataset, not generated from a timestamp. + +One approximation is worth flagging. On a daily index, `Fourier(period=365.25)` +counts whole days, so the harmonic drifts against the calendar within a leap +cycle. It is fine over multi-year samples and exact for monthly, quarterly and +annual sampling, where a period ordinal *is* the calendar unit. diff --git a/docs/how-to/index.md b/docs/how-to/index.md index 7ff315b..64d7838 100644 --- a/docs/how-to/index.md +++ b/docs/how-to/index.md @@ -6,6 +6,7 @@ Practical recipes for solving specific problems with Impulso. :maxdepth: 1 data-preparation +deterministic-regressors custom-priors lag-selection stationarity-testing diff --git a/docs/reference/deterministic.md b/docs/reference/deterministic.md new file mode 100644 index 0000000..a6de051 --- /dev/null +++ b/docs/reference/deterministic.md @@ -0,0 +1,25 @@ +# Deterministic Regressors + +Calendar-anchored exogenous design matrices: `Trend`, `Fourier`, +`SeasonalDummies` and `BreakDummy`, composed by `DeterministicDesign` into the +exogenous block of a `VARData`. Elapsed time is counted in integer period +ordinals from the first timestamp of the estimation index, so `extend` writes +exactly the rows `build` would have written on a longer sample — which is what +lets `exog_future` hand `FittedVAR.forecast` a column-aligned block. + +See the how-to guide, [Deterministic Regressors for Climate +VARs](../how-to/deterministic-regressors.md), for the recipe. + +```{eval-rst} +.. currentmodule:: impulso.deterministic + +.. autosummary:: + :toctree: generated/ + :nosignatures: + + DeterministicDesign + Trend + Fourier + SeasonalDummies + BreakDummy +``` diff --git a/docs/reference/index.md b/docs/reference/index.md index 100e069..ab911cd 100644 --- a/docs/reference/index.md +++ b/docs/reference/index.md @@ -18,6 +18,7 @@ fitted identified identification scenario +deterministic results evidence primitives diff --git a/docs/reference/protocols.md b/docs/reference/protocols.md index 57284a5..fdc80ad 100644 --- a/docs/reference/protocols.md +++ b/docs/reference/protocols.md @@ -12,4 +12,5 @@ IdentificationScheme VolatilityProcess ErrorDistribution + DeterministicTerm ``` diff --git a/src/impulso/__init__.py b/src/impulso/__init__.py index 8fe3f84..33454c6 100644 --- a/src/impulso/__init__.py +++ b/src/impulso/__init__.py @@ -14,13 +14,20 @@ from impulso._ma import compute_ma_phi from impulso.conjugate import ConjugateVAR from impulso.conjugate_volatility import ConjugateVolatility, PandemicBreak + from impulso.deterministic import ( + BreakDummy, + DeterministicDesign, + Fourier, + SeasonalDummies, + Trend, + ) from impulso.evidence import EvidenceComparison, ModelEvidence, compare_evidence from impulso.fitted import FittedVAR from impulso.identification import Cholesky, LongRunRestriction, ProxySVAR, SignRestriction from impulso.identified import IdentifiedVAR from impulso.observation import Gaussian, StudentT from impulso.priors import MinnesotaPrior, NIWPrior - from impulso.protocols import ErrorDistribution, VolatilityProcess + from impulso.protocols import DeterministicTerm, ErrorDistribution, VolatilityProcess from impulso.results import ( CointegrationTestResult, ConditionalForecastResult, @@ -48,6 +55,7 @@ __all__ = [ "VAR", + "BreakDummy", "Cholesky", "CointegrationTestResult", "ConditionalForecastResult", @@ -55,6 +63,8 @@ "ConjugateVolatility", "Constant", "CounterfactualResult", + "DeterministicDesign", + "DeterministicTerm", "DynamicMultiplierResult", "ErrorDistribution", "EvidenceComparison", @@ -62,6 +72,7 @@ "FittedSV", "FittedVAR", "ForecastResult", + "Fourier", "Gaussian", "HDIResult", "HistoricalDecompositionResult", @@ -80,11 +91,13 @@ "SVDefaultPrior", "SVForecastResult", "ScenarioResult", + "SeasonalDummies", "ShockPath", "SignRestriction", "StationarityTestResult", "StochasticVolatility", "StudentT", + "Trend", "VARData", "VariablePath", "VolatilityProcess", @@ -113,6 +126,12 @@ "ConjugateVAR": "impulso.conjugate", "ConjugateVolatility": "impulso.conjugate_volatility", "PandemicBreak": "impulso.conjugate_volatility", + "DeterministicDesign": "impulso.deterministic", + "Trend": "impulso.deterministic", + "Fourier": "impulso.deterministic", + "SeasonalDummies": "impulso.deterministic", + "BreakDummy": "impulso.deterministic", + "DeterministicTerm": "impulso.protocols", "ModelEvidence": "impulso.evidence", "EvidenceComparison": "impulso.evidence", "compare_evidence": "impulso.evidence", diff --git a/src/impulso/deterministic.py b/src/impulso/deterministic.py new file mode 100644 index 0000000..b70804f --- /dev/null +++ b/src/impulso/deterministic.py @@ -0,0 +1,667 @@ +"""Deterministic regressors — calendar-anchored exogenous design matrices. + +Climate series arrive with structure that is not dynamics: an annual +cycle, a warming trend, a step change when an instrument was replaced. +Left in the endogenous block, a VAR spends its coefficients re-learning +the calendar. This module builds those features as *deterministic* +regressors — total functions of a timestamp — which enter the model +through `VARData`'s exogenous block and come back with posterior +coefficients you can read. + +The pieces are `Trend`, `Fourier`, `SeasonalDummies` and `BreakDummy`, +composed by `DeterministicDesign`: + +```python +design = DeterministicDesign( + terms=[Trend(degree=1, scale=120.0), Fourier(period=12, order=2)] +) +frame = pd.concat([anomalies, design.build(anomalies.index)], axis=1) +data = VARData.from_df(frame, endog=list(anomalies.columns), exog=design.column_names) +``` + +Everything is anchored to the calendar rather than to row position: +elapsed time is counted in integer period ordinals from the first +timestamp of the *estimation* index. That is what makes the +continuation property hold — + + design.build(index[: T + h]).iloc[T:] == design.extend(index[:T], h) + +— so `design.exog_future(fitted, h)` is exactly the block the estimated +coefficients were fitted against, in the column order the posterior +expects. + +Out of scope (deliberately): holiday and business calendars; interaction +terms; slope breaks; inferring cycle lengths from the data; and any +regressor that is *data* rather than calendar arithmetic — solar +forcing, ENSO indices, CO2 concentrations belong in a dataset, not here. +""" + +from __future__ import annotations + +import datetime as _datetime +from typing import TYPE_CHECKING, Literal + +import numpy as np +import pandas as pd +from pydantic import Field, field_validator, model_validator + +from impulso._base import ImpulsoBaseModel, ImpulsoModel +from impulso.data import VARData +from impulso.protocols import DeterministicTerm + +if TYPE_CHECKING: + from impulso.fitted import FittedVAR + +_TREND_NAMES = ("trend", "trend_squared", "trend_cubed") + +_SEASON_ATTR = {"month": "month", "quarter": "quarter", "dayofweek": "dayofweek"} +_SEASON_LEVELS = { + "month": tuple(range(1, 13)), + "quarter": tuple(range(1, 5)), + "dayofweek": tuple(range(7)), +} +_SEASON_PREFIX = {"month": "month", "quarter": "quarter", "dayofweek": "dow"} +_SEASON_CYCLE = {"month": 12.0, "quarter": 4.0, "dayofweek": 7.0} + + +# --------------------------------------------------------------------------- # +# Time helpers +# --------------------------------------------------------------------------- # + + +def _as_datetime_index(index: pd.DatetimeIndex, label: str = "index") -> pd.DatetimeIndex: + """Validate that `index` is a usable, strictly increasing DatetimeIndex.""" + if not isinstance(index, pd.DatetimeIndex): + raise TypeError(f"{label} must be a pandas DatetimeIndex, got {type(index).__name__}") + if len(index) == 0: + raise ValueError(f"{label} must not be empty") + if not index.is_monotonic_increasing or not index.is_unique: + raise ValueError(f"{label} must be strictly increasing with no duplicate timestamps") + return index + + +def _resolve_offset(index: pd.DatetimeIndex, freq: str | None) -> pd.offsets.BaseOffset: + """Resolve the sampling offset from an explicit alias, the index, or inference. + + Inference is *verified*: `pd.infer_freq` returns confident answers on + short irregular indices, so the candidate offset must regenerate the + index exactly before it is accepted. + """ + if freq is not None: + return pd.tseries.frequencies.to_offset(freq) + if index.freq is not None: + return index.freq + inferred = pd.infer_freq(index) if len(index) >= 3 else None + if inferred is None: + raise ValueError( + "Could not determine the sampling frequency of the index: it carries no " + "`freq` and pandas could not infer one. Pass it explicitly, e.g. " + 'DeterministicDesign(terms=[...], freq="MS").' + ) + offset = pd.tseries.frequencies.to_offset(inferred) + regenerated = pd.date_range(index[0], periods=len(index), freq=offset) + if not regenerated.equals(index): + raise ValueError( + f"pandas inferred freq={inferred!r} for this index, but regenerating the " + "index from that frequency does not reproduce it — the inference is a " + "false positive on an irregular index. Pass the true sampling frequency " + 'explicitly, e.g. DeterministicDesign(terms=[...], freq="MS").' + ) + return offset + + +def _period_alias(offset: pd.offsets.BaseOffset) -> str: + """Map a sampling offset to the pandas period alias used for ordinals. + + Period ordinals are the anchor for every elapsed-time count: they are + exact integers, unit-agnostic, and immune to the `datetime64[us]` vs + `datetime64[ns]` resolution traps of raw integer timestamp arithmetic. + """ + if isinstance(offset, pd.offsets.BusinessDay | pd.offsets.CustomBusinessDay): + raise ValueError( # noqa: TRY004 — an unusable frequency, not a wrong type + "Business-day frequencies have no period equivalent in pandas " + "(PeriodDtype[B] is deprecated), so deterministic terms cannot be " + 'anchored to them. Use freq="D" — calendar-day ordinals count ' + "business days correctly as long as the index itself is business-daily." + ) + try: + return pd.date_range(pd.Timestamp("2000-01-03"), periods=1, freq=offset).to_period().freqstr + except (ValueError, TypeError) as exc: + raise ValueError( + f"Frequency {offset.freqstr!r} has no pandas period equivalent, so " + "deterministic terms cannot be anchored to the calendar at this " + "sampling rate. Use a period-compatible frequency (MS, ME, QS, QE, " + "YS, YE, D, W, h, ...)." + ) from exc + + +def _elapsed(index: pd.DatetimeIndex, origin: pd.Timestamp, alias: str) -> np.ndarray: + """Sampling periods elapsed between `origin` and each timestamp, as float64. + + Calendar arithmetic, not row counting: a gap in the index produces a + jump in the returned values, which is the correct behaviour for a + trend (time really did pass) and for a harmonic (the phase really did + advance). + + The unit is one *sampling period*, not one period-ordinal tick. Those + differ whenever the frequency carries a multiplier: pandas stores + `15D` ordinals in days and `2h` ordinals in hours, so the raw + difference would count 15 (or 2) per observation. Dividing by the + multiplier is exact — the origin is subtracted first, so an on-grid + index yields exact integers — and it is what makes `Fourier.period` + mean what its docstring says at every sampling rate. + """ + anchor = pd.Period(origin, freq=alias) + ticks = (index.to_period(alias).asi8 - anchor.ordinal).astype(np.float64) + return ticks / anchor.freq.n + + +def _extend_index(index: pd.DatetimeIndex, offset: pd.offsets.BaseOffset, steps: int) -> pd.DatetimeIndex: + """The `steps` timestamps following the last entry of `index`. + + `pd.date_range` rolls a start that is not on the offset's anchor + forward to the next valid date, so the walk's first entry is *already* + a future period whenever the sample ends off-anchor — an irregular + index with an explicit `freq`, which this module supports. Dropping it + unconditionally would skip a period (a sample ending 2000-03-15 under + `MS` would forecast from May, silently losing April). Selecting the + entries strictly after the last observation is correct either way. + """ + walk = pd.date_range(index[-1], freq=offset, periods=steps + 1) + return walk[walk > index[-1]][:steps] + + +def _format_period(period: float) -> str: + """Render a cycle length for a column name: 12.0 -> "12", 365.25 -> "365.25".""" + return str(int(period)) if float(period).is_integer() else str(float(period)) + + +# --------------------------------------------------------------------------- # +# Terms +# --------------------------------------------------------------------------- # + + +class Trend(ImpulsoModel): + """Polynomial time trend anchored to the start of the estimation sample. + + Column `trend` is the number of periods elapsed since the first + observation, divided by `scale`; higher degrees are integer powers of + that same column. + + The prior on the exogenous coefficients adapts to each regressor's + sample spread (see `VAR.exog_prior_scale`), so an unscaled trend is + no longer fought by the prior. `scale` still matters for two other + reasons: the coefficient's units, and the sampler's geometry. Divide + by a fixed, interpretable constant — periods per decade, say — so + the coefficient reads as "change per decade", and so the design + column stays O(1) rather than reaching 539 by the end of a 540-month + sample. + + Warning: + `scale` must not depend on the sample length. A `T`-dependent + scale (`scale=len(index)`) breaks the continuation property: the + design you extend with would no longer be the design you fitted. + + Attributes: + degree: Highest power of elapsed time, 1 to 3. + scale: Divisor applied to elapsed time before exponentiation. + """ + + degree: int = Field(default=1, ge=1, le=3) + scale: float = Field(default=1.0, gt=0) + + @property + def column_names(self) -> list[str]: + """Column names contributed by this term.""" + return list(_TREND_NAMES[: self.degree]) + + def build(self, index: pd.DatetimeIndex, origin: pd.Timestamp, alias: str) -> np.ndarray: + """Evaluate the trend powers on `index`, anchored at `origin`.""" + t = _elapsed(index, origin, alias) / self.scale + return np.column_stack([t**k for k in range(1, self.degree + 1)]) + + +class Fourier(ImpulsoModel): + """Harmonic pair(s) representing a smooth cycle of known length. + + Order `k` contributes `sin(2πk·t/period)` and `cos(2πk·t/period)`, + where `t` is periods elapsed since the start of the estimation + sample. Two harmonics on a monthly index (`period=12, order=2`) buy a + seasonal shape at four coefficients per variable instead of the + eleven that month dummies cost. + + The cycle length is always explicit — nothing here infers a period + from the data. + + Note: + On a daily index, `period=365.25` is an approximation: elapsed + time is counted in whole days, so the harmonic drifts against the + calendar within a leap cycle. It is fine for multi-year samples + and exact for monthly, quarterly and annual sampling, where a + period ordinal *is* the calendar unit. + + Attributes: + period: Cycle length in sampling periods (12 for an annual cycle + on monthly data, 4 on quarterly data). + order: Number of harmonic pairs. Must satisfy `2 * order <= + period` — beyond the Nyquist limit the extra harmonics are + not identified from the sampled points. + """ + + period: float = Field(..., gt=1) + order: int = Field(..., ge=1) + + @model_validator(mode="after") + def _validate_nyquist(self) -> Fourier: + if 2 * self.order > self.period: + raise ValueError( + f"Fourier(period={self.period}, order={self.order}) exceeds the Nyquist " + f"limit: at most {int(self.period // 2)} harmonic pairs are identified " + f"from a cycle sampled {self.period} times." + ) + return self + + @property + def column_names(self) -> list[str]: + """Column names contributed by this term, sine before cosine per order.""" + label = _format_period(self.period) + names: list[str] = [] + for k in range(1, self.order + 1): + names.append(f"sin({k},{label})") + names.append(f"cos({k},{label})") + return names + + def build(self, index: pd.DatetimeIndex, origin: pd.Timestamp, alias: str) -> np.ndarray: + """Evaluate the harmonics on `index`, with phase anchored at `origin`.""" + t = _elapsed(index, origin, alias) + columns: list[np.ndarray] = [] + for k in range(1, self.order + 1): + angle = 2.0 * np.pi * k * t / self.period + columns.append(np.sin(angle)) + columns.append(np.cos(angle)) + return np.column_stack(columns) + + +class SeasonalDummies(ImpulsoModel): + """Indicator columns for a calendar season, one level dropped. + + Unlike trends and harmonics this term reads the calendar attribute + directly, so it needs no origin and is invariant to where the sample + starts. + + One level must be dropped: `VAR` and `ConjugateVAR` both fit an + unconditional intercept, and a full set of indicators sums to it. + + Attributes: + season: `"month"` (levels 1-12), `"quarter"` (1-4) or + `"dayofweek"` (0-6, Monday = 0). + drop_first: Drop one level to keep the design full rank against + the intercept. Leave `True` unless the design matrix is bound + for an estimator that fits no intercept. + reference: The level to drop. Defaults to the first level of the + season. Only meaningful when `drop_first` is `True`. + """ + + season: Literal["month", "quarter", "dayofweek"] + drop_first: bool = True + reference: int | None = None + + @model_validator(mode="after") + def _validate_reference(self) -> SeasonalDummies: + if self.reference is None: + return self + if not self.drop_first: + raise ValueError("reference is only meaningful with drop_first=True — nothing is dropped otherwise.") + levels = _SEASON_LEVELS[self.season] + if self.reference not in levels: + raise ValueError( + f"reference={self.reference} is not a valid {self.season} level; expected one of {list(levels)}." + ) + return self + + @property + def _kept_levels(self) -> tuple[int, ...]: + levels = _SEASON_LEVELS[self.season] + if not self.drop_first: + return levels + dropped = levels[0] if self.reference is None else self.reference + return tuple(level for level in levels if level != dropped) + + @property + def column_names(self) -> list[str]: + """Column names contributed by this term, in ascending level order.""" + prefix = _SEASON_PREFIX[self.season] + return [f"{prefix}_{level}" for level in self._kept_levels] + + def build(self, index: pd.DatetimeIndex, origin: pd.Timestamp, alias: str) -> np.ndarray: + """Evaluate the indicators on `index`. `origin` and `alias` are unused.""" + del origin, alias + values = np.asarray(getattr(index, _SEASON_ATTR[self.season])) + return np.column_stack([(values == level).astype(np.float64) for level in self._kept_levels]) + + +class BreakDummy(ImpulsoBaseModel): + """Indicator for a known, dated discontinuity. + + A `"level"` break is 1 from `date` onward — an instrument change, a + station move, a regime shift. Its coefficient is the permanent offset, + and `FittedVAR.dynamic_multiplier(cumulative=True)` propagates it + through the lag dynamics into the full adjustment path. + + A `"pulse"` break is 1 on `date` alone — a one-off event whose + influence you want held out of the dynamics. + + Attributes: + date: Timestamp of the break. Strings, `datetime` objects and + `numpy.datetime64` are coerced. + kind: `"level"` (step, the default) or `"pulse"` (one period). + """ + + date: pd.Timestamp + kind: Literal["level", "pulse"] = "level" + + @field_validator("date", mode="before") + @classmethod + def _coerce_timestamps(cls, value: object) -> object: + if isinstance(value, pd.Timestamp): + return value + if isinstance(value, str | _datetime.date | np.datetime64): + return pd.Timestamp(value) + return value + + @property + def column_names(self) -> list[str]: + """The single column name contributed by this term.""" + return [f"{self.kind}_{self.date.strftime('%Y-%m-%d')}"] + + def build(self, index: pd.DatetimeIndex, origin: pd.Timestamp, alias: str) -> np.ndarray: + """Evaluate the indicator on `index`. `origin` and `alias` are unused.""" + del origin, alias + hit = index >= self.date if self.kind == "level" else index == self.date + return np.asarray(hit, dtype=np.float64).reshape(-1, 1) + + def _check_in_sample(self, index: pd.DatetimeIndex) -> None: + """Raise if the break is not identified from the estimation sample.""" + if self.kind == "pulse": + if self.date not in index: + position = int(index.searchsorted(self.date)) + before = index[position - 1] if position > 0 else None + after = index[position] if position < len(index) else None + raise ValueError( + f"BreakDummy(kind='pulse', date={self.date.date()}) does not fall on " + f"the sampled index, so the column is zero everywhere and its " + f"coefficient is not identified. Nearest sampled dates: " + f"{before.date() if before is not None else None} (before) and " + f"{after.date() if after is not None else None} (after)." + ) + return + if self.date <= index[0]: + raise ValueError( + f"BreakDummy(kind='level', date={self.date.date()}) is at or before the " + f"start of the sample ({index[0].date()}), so the column is 1 everywhere " + f"and is collinear with the intercept every estimator fits." + ) + if self.date > index[-1]: + raise ValueError( + f"BreakDummy(kind='level', date={self.date.date()}) is after the end of " + f"the sample ({index[-1].date()}), so the shift never occurs in-sample " + f"and its coefficient is not identified." + ) + + +# --------------------------------------------------------------------------- # +# Design +# --------------------------------------------------------------------------- # + + +class DeterministicDesign(ImpulsoBaseModel): + """A composed set of deterministic terms, built and extended as one block. + + The design owns the calendar: it resolves a single `(origin, alias)` + pair from the estimation index and passes it to every term, in-sample + and out. That is what makes the continuation property hold — + + design.build(index[: T + h]).iloc[T:] == design.extend(index[:T], h) + + — and it is why `exog_future` can hand `FittedVAR.forecast` a block + that matches, column for column, what the model was fitted against. + + Attributes: + terms: The `DeterministicTerm` instances to concatenate, in + column order. At least one is required. + freq: Explicit pandas frequency alias for the sampling rate (e.g. + `"MS"`, `"QS"`, `"D"`). Optional: the design falls back to the + index's own `freq`, then to *verified* inference. Pass it + whenever the index has gaps — calendar anchoring handles gaps + correctly, but pandas cannot infer a frequency through them. + """ + + terms: tuple[DeterministicTerm, ...] + freq: str | None = None + + @model_validator(mode="after") + def _validate_terms(self) -> DeterministicDesign: + if not self.terms: + raise ValueError("DeterministicDesign requires at least one term.") + seen: dict[str, int] = {} + for position, term in enumerate(self.terms): + for name in term.column_names: + if name in seen: + raise ValueError( + f"Duplicate column name {name!r}: terms {seen[name]} " + f"({type(self.terms[seen[name]]).__name__}) and {position} " + f"({type(term).__name__}) both emit it. Deterministic column " + f"names are the contract that keeps forecast blocks aligned, so " + f"they must be unique." + ) + seen[name] = position + return self + + @property + def column_names(self) -> list[str]: + """Every column the design emits, in build order.""" + return [name for term in self.terms for name in term.column_names] + + # -- calendar resolution ------------------------------------------------ # + + def _resolve_calendar(self, index: pd.DatetimeIndex) -> tuple[pd.Timestamp, str, pd.offsets.BaseOffset]: + offset = _resolve_offset(index, self.freq) + return index[0], _period_alias(offset), offset + + def _assemble(self, index: pd.DatetimeIndex, origin: pd.Timestamp, alias: str) -> pd.DataFrame: + blocks: list[np.ndarray] = [] + for term in self.terms: + block = np.asarray(term.build(index, origin, alias), dtype=np.float64) + expected = (len(index), len(term.column_names)) + if block.shape != expected: + raise ValueError( + f"{type(term).__name__}.build returned shape {block.shape}, but its column_names imply {expected}." + ) + blocks.append(block) + return pd.DataFrame(np.hstack(blocks), index=index, columns=self.column_names) + + # -- public surface ----------------------------------------------------- # + + def build(self, index: pd.DatetimeIndex) -> pd.DataFrame: + """Construct the in-sample design matrix. + + The result is checked for rank deficiency against an intercept + column, because every Impulso estimator fits one and a deficient + exogenous block is not identified. + + Args: + index: The estimation index — strictly increasing, no + duplicates. Its first timestamp becomes the origin for + every elapsed-time count. + + Returns: + A `float64` DataFrame indexed by `index`, with columns + `self.column_names` in that order and no missing values. + + Raises: + ValueError: If the frequency cannot be resolved, a break + dummy is not identified in-sample, there are fewer + observations than columns, or the design is collinear. + """ + index = _as_datetime_index(index) + origin, alias, _ = self._resolve_calendar(index) + for term in self.terms: + if isinstance(term, BreakDummy): + term._check_in_sample(index) + frame = self._assemble(index, origin, alias) + self._check_rank(frame) + return frame + + def extend( + self, + index: pd.DatetimeIndex, + steps: int, + future_index: pd.DatetimeIndex | None = None, + ) -> pd.DataFrame: + """Construct the design matrix for `steps` periods past `index`. + + Takes the *estimation* index, not the future one, so the origin + and frequency resolve exactly as they did in `build` (the + statsmodels `out_of_sample` contract). No rank check is applied: + short horizons are legitimately rank-deficient — a 3-step block + cannot span twelve month dummies — and nothing is being estimated + from these rows. + + Args: + index: The estimation index the design was built on. + steps: Number of future periods, at least 1. + future_index: Optional explicit future timestamps, overriding + the frequency walk. Must have length `steps`. + + Returns: + A `float64` DataFrame of `steps` rows with the same columns, + in the same order, as `build`. + """ + index = _as_datetime_index(index) + if steps < 1: + raise ValueError(f"steps must be >= 1, got {steps}") + origin, alias, offset = self._resolve_calendar(index) + if future_index is None: + resolved = _extend_index(index, offset, steps) + else: + resolved = _as_datetime_index(future_index, label="future_index") + if len(resolved) != steps: + raise ValueError(f"future_index has length {len(resolved)}, but steps={steps}.") + return self._assemble(resolved, origin, alias) + + def future_index(self, index: pd.DatetimeIndex, steps: int) -> pd.DatetimeIndex: + """The timestamps `extend` would use for `steps` periods past `index`. + + Args: + index: The estimation index. + steps: Number of future periods, at least 1. + + Returns: + A `DatetimeIndex` of length `steps`. + """ + index = _as_datetime_index(index) + if steps < 1: + raise ValueError(f"steps must be >= 1, got {steps}") + _, _, offset = self._resolve_calendar(index) + return _extend_index(index, offset, steps) + + def exog_future( + self, + fitted: FittedVAR | VARData, + steps: int, + future_index: pd.DatetimeIndex | None = None, + ) -> np.ndarray: + """The future exogenous block for `FittedVAR.forecast`, column-aligned. + + `forecast`, `conditional_forecast` and `structural_scenario` all + index `exog_future` positionally, so a permuted column order is a + silently wrong forecast rather than an error. This method removes + that hazard: it extends the design, then reorders the columns by + *name* to match `exog_names` on the data the model was fitted + with. + + Args: + fitted: A `FittedVAR`, or the `VARData` it was fitted on. + steps: Forecast horizon, at least 1. + future_index: Optional explicit future timestamps, of length + `steps`. + + Returns: + A `(steps, k)` float64 array whose columns are ordered as + `exog_names`. + + Raises: + ValueError: If the data carries no exogenous block, or if the + design's columns and `exog_names` are not the same set. + """ + data = fitted if isinstance(fitted, VARData) else fitted.data + if data.exog_names is None: + raise ValueError( + "This model was fitted without exogenous regressors, so there is no " + "future exogenous block to build. Refit with the design included in " + "VARData before forecasting with it." + ) + frame = self.extend(data.index, steps, future_index=future_index) + missing = [name for name in data.exog_names if name not in frame.columns] + extra = [name for name in frame.columns if name not in data.exog_names] + if missing or extra: + raise ValueError( + f"The design does not match the fitted exogenous block. Fitted " + f"exog_names: {list(data.exog_names)}; design columns: " + f"{list(frame.columns)}. Missing from the design: {missing}; not fitted: " + f"{extra}. Build the forecast block with the same design the model was " + f"fitted with." + ) + return np.ascontiguousarray(frame[list(data.exog_names)].to_numpy(dtype=np.float64)) + + # -- rank diagnostics --------------------------------------------------- # + + def _check_rank(self, frame: pd.DataFrame) -> None: + design = frame.to_numpy(dtype=np.float64) + augmented = np.column_stack([np.ones(len(frame)), design]) + n_rows, n_cols = augmented.shape + if n_rows < n_cols: + raise ValueError( + f"Too few observations for this design: {len(frame)} rows against " + f"{design.shape[1]} deterministic columns plus the intercept every " + f"estimator fits. Shorten the design or lengthen the sample." + ) + if int(np.linalg.matrix_rank(augmented)) == n_cols: + return + hints = self._collinearity_hints() + detail = ("\n - " + "\n - ".join(hints)) if hints else "" + raise ValueError( + "This deterministic design has collinear columns once the intercept every " + f"estimator fits is included ({n_cols} columns, rank " + f"{int(np.linalg.matrix_rank(augmented))}), so its coefficients are not " + f"identified.{detail}" + ) + + def _collinearity_hints(self) -> list[str]: + hints: list[str] = [] + seasons = [t for t in self.terms if isinstance(t, SeasonalDummies)] + fouriers = [t for t in self.terms if isinstance(t, Fourier)] + for term in seasons: + if not term.drop_first: + hints.append( + f"SeasonalDummies(season={term.season!r}, drop_first=False) emits every " + f"level, and they sum to the intercept column. Set drop_first=True." + ) + for term in fouriers: + if 2 * term.order == term.period: + hints.append( + f"Fourier(period={term.period}, order={term.order}) sits exactly at the " + f"Nyquist limit: sin(2*pi*{term.order}*t/{term.period}) is identically " + f"zero at integer t. Drop the last harmonic pair (order=" + f"{term.order - 1})." + ) + for season in seasons: + for fourier in fouriers: + if fourier.period == _SEASON_CYCLE[season.season]: + hints.append( + f"SeasonalDummies(season={season.season!r}) and Fourier(period=" + f"{fourier.period}) describe the same cycle; harmonics of a cycle " + f"sampled that many times are linear combinations of its level " + f"indicators. Keep one or the other." + ) + return hints diff --git a/src/impulso/fitted.py b/src/impulso/fitted.py index f329a61..ed8aeef 100644 --- a/src/impulso/fitted.py +++ b/src/impulso/fitted.py @@ -153,6 +153,32 @@ def innovation_covariance(self) -> np.ndarray: extra_dims = sigma.ndim - inflation.ndim return sigma * inflation.reshape(inflation.shape + (1,) * extra_dims) + def _resolve_exog_future(self, exog_future: np.ndarray | None, steps: int) -> np.ndarray | None: + """Validate a future exogenous block against the fitted posterior. + + Checked here rather than left to the propagation `einsum`: the block + is routinely generated now (`DeterministicDesign.exog_future`), and a + mis-shaped or mis-ordered one must not surface as an opaque einsum + error. `conditional_forecast` enforces the same contract. + """ + if self.has_exog and exog_future is None: + raise ValueError("exog_future is required when model includes exogenous variables") + if not self.has_exog and exog_future is not None: + raise ValueError("exog_future provided but model has no exogenous variables") + if exog_future is None: + return None + if "B_exog" not in self.idata.posterior: + raise ValueError( + "This FittedVAR's data carries exogenous regressors the estimator " + "never consumed (no B_exog in the posterior); refit with an " + "estimator that supports them before forecasting." + ) + block = np.asarray(exog_future, dtype=float) + n_exog = self.idata.posterior["B_exog"].shape[-1] + if block.shape != (steps, n_exog): + raise ValueError(f"exog_future must have shape ({steps}, {n_exog}), got {block.shape}.") + return block + def forecast( self, steps: int, @@ -192,10 +218,7 @@ def forecast( from impulso.results import ForecastResult - if self.has_exog and exog_future is None: - raise ValueError("exog_future is required when model includes exogenous variables") - if not self.has_exog and exog_future is not None: - raise ValueError("exog_future provided but model has no exogenous variables") + exog_future = self._resolve_exog_future(exog_future, steps) rng = np.random.default_rng(seed) if not isinstance(seed, np.random.Generator) else seed diff --git a/src/impulso/protocols.py b/src/impulso/protocols.py index e6a99dc..cb7983b 100644 --- a/src/impulso/protocols.py +++ b/src/impulso/protocols.py @@ -6,6 +6,7 @@ if TYPE_CHECKING: import arviz as az + import pandas as pd import pymc as pm import pytensor.tensor as pt import xarray as xr @@ -27,6 +28,60 @@ class Sampler(Protocol): def sample(self, model: "pm.Model") -> "az.InferenceData": ... +@runtime_checkable +class DeterministicTerm(Protocol): + """Contract for deterministic regressor terms. + + A term is a total function of a timestamp: given an index it returns + real-valued columns with no missing entries, and it does so from the + calendar alone — never from the endogenous data. `DeterministicDesign` + composes terms into an exogenous design matrix. + + The `(origin, alias)` pair is resolved **once** from the estimation + index and handed to every term, for both in-sample construction and + out-of-sample extension. Terms that count elapsed time (trends, + harmonics) must anchor on it rather than on the position of a row + inside the index they are handed; that is what makes + `design.extend(index, h)` reproduce the rows `design.build` would + have written had the index been `h` periods longer. + + Concrete implementations: `Trend`, `Fourier`, `SeasonalDummies`, + `BreakDummy` (all in `impulso.deterministic`). + """ + + @property + def column_names(self) -> list[str]: + """Names of the columns this term contributes, in build order. + + Must be knowable without an index — the design's column contract + is static — and must have the same length as `build`'s second + axis. + """ + ... + + def build(self, index: "pd.DatetimeIndex", origin: "pd.Timestamp", alias: str) -> np.ndarray: + """Evaluate the term on `index`. + + Args: + index: Timestamps to evaluate at. Not necessarily the + estimation index — `DeterministicDesign.extend` passes a + future index while keeping `origin` and `alias` fixed. + origin: First timestamp of the estimation index; the zero + point for elapsed-time counts. + alias: pandas period alias for the sampling frequency (e.g. + `"M"`, `"Q-DEC"`, `"D"`, `"15D"`), used to convert + timestamps to integer period ordinals. Note that a + multiplied alias stores ordinals in its *base* unit, so + elapsed time must be divided by the multiplier to be + counted in sampling periods. + + Returns: + Float array of shape `(len(index), len(self.column_names))` + containing only finite values. + """ + ... + + @runtime_checkable class IdentificationScheme(Protocol): """Contract for structural identification schemes. diff --git a/tests/test_deterministic.py b/tests/test_deterministic.py new file mode 100644 index 0000000..289a2bb --- /dev/null +++ b/tests/test_deterministic.py @@ -0,0 +1,807 @@ +"""Tests for deterministic regressors (`impulso.deterministic`). + +The headline property is *continuation*: the rows `extend` writes for the +future must be exactly the rows `build` would have written had the sample +run that much longer. Everything else — trend anchoring, harmonic phase, +dummy calendars, column order, dtypes, the future index — falls out of it, +so `TestContinuationProperty` carries most of the weight here. +""" + +import numpy as np +import pandas as pd +import pandas.testing as pdt +import pytest +from numpy.testing import assert_allclose + +from impulso.data import VARData +from impulso.deterministic import ( + BreakDummy, + DeterministicDesign, + Fourier, + SeasonalDummies, + Trend, + _format_period, +) + +# --------------------------------------------------------------------------- # +# Index fixtures +# --------------------------------------------------------------------------- # + +_FREQ_SPECS = { + "MS": ("1980-01-01", 540), + "QS": ("1980-01-01", 240), + "D": ("2010-01-01", 1200), + "YS": ("1900-01-01", 150), +} + +#: Reserved tail length: continuation cases build on `index[: T + h]`. +_TAIL = 25 + + +def _index_for(freq_key: str) -> pd.DatetimeIndex: + start, periods = _FREQ_SPECS[freq_key] + return pd.date_range(start, periods=periods, freq=freq_key) + + +#: Per-frequency designs that are full rank in-sample, keyed by a test id. +_DESIGNS_BY_FREQ: dict[str, dict[str, list]] = { + "MS": { + "trend1": [Trend(degree=1, scale=120.0)], + "trend3": [Trend(degree=3, scale=120.0)], + "fourier": [Fourier(period=12, order=2)], + "month": [SeasonalDummies(season="month")], + "level": [BreakDummy(date="2000-01-01")], + "pulse": [BreakDummy(date="2000-01-01", kind="pulse")], + "composed": [ + Trend(degree=1, scale=120.0), + Fourier(period=12, order=2), + SeasonalDummies(season="quarter"), + BreakDummy(date="2000-01-01"), + ], + }, + "QS": { + "trend2": [Trend(degree=2, scale=40.0)], + "fourier": [Fourier(period=4, order=1)], + "quarter": [SeasonalDummies(season="quarter")], + "pulse": [BreakDummy(date="2000-01-01", kind="pulse")], + "composed": [ + Trend(degree=1, scale=40.0), + SeasonalDummies(season="quarter"), + BreakDummy(date="2000-01-01"), + ], + }, + "D": { + "trend1": [Trend(degree=1, scale=365.0)], + "fourier": [Fourier(period=365.25, order=2)], + "dayofweek": [SeasonalDummies(season="dayofweek")], + "level": [BreakDummy(date="2012-01-01")], + "composed": [ + Trend(degree=1, scale=365.0), + Fourier(period=365.25, order=2), + SeasonalDummies(season="dayofweek"), + BreakDummy(date="2012-01-01"), + ], + }, + "YS": { + "trend1": [Trend(degree=1, scale=10.0)], + "fourier": [Fourier(period=11, order=2)], + "level": [BreakDummy(date="1950-01-01")], + "composed": [Trend(degree=1, scale=10.0), BreakDummy(date="1950-01-01")], + }, +} + + +def _continuation_cases(): + for freq_key, designs in _DESIGNS_BY_FREQ.items(): + for design_id, terms in designs.items(): + for h in (1, 3, 13, 25): + yield pytest.param(freq_key, terms, h, id=f"{freq_key}-{design_id}-h{h}") + + +@pytest.fixture +def monthly_index(): + return _index_for("MS") + + +@pytest.fixture +def quarterly_index(): + return _index_for("QS") + + +@pytest.fixture +def daily_index(): + return _index_for("D") + + +@pytest.fixture +def annual_index(): + return _index_for("YS") + + +@pytest.fixture +def gappy_monthly_index(monthly_index): + """Monthly index with a six-month hole and therefore no inferable freq.""" + return monthly_index.delete(range(36, 42)) + + +@pytest.fixture +def simple_design(): + """All four term types, jointly full rank on a monthly index.""" + return DeterministicDesign( + terms=[ + Trend(degree=1, scale=120.0), + Fourier(period=12, order=2), + SeasonalDummies(season="quarter"), + BreakDummy(date="2000-01-01"), + ] + ) + + +# --------------------------------------------------------------------------- # +# A. The continuation property +# --------------------------------------------------------------------------- # + + +class TestContinuationProperty: + """`extend` must reproduce the rows `build` writes on a longer index.""" + + @pytest.mark.parametrize(("freq_key", "terms", "h"), _continuation_cases()) + def test_extend_matches_build_on_longer_index(self, freq_key, terms, h): + index = _index_for(freq_key) + design = DeterministicDesign(terms=terms) + cut = len(index) - _TAIL + + expected = design.build(index[: cut + h]).iloc[cut:] + actual = design.extend(index[:cut], h) + + pdt.assert_frame_equal(expected, actual) + + @pytest.mark.parametrize(("freq_key", "terms", "h"), _continuation_cases()) + def test_future_index_matches_extend_index(self, freq_key, terms, h): + index = _index_for(freq_key) + design = DeterministicDesign(terms=terms) + cut = len(index) - _TAIL + + assert design.extend(index[:cut], h).index.equals(design.future_index(index[:cut], h)) + + def test_columns_and_dtypes_are_stable(self, monthly_index, simple_design): + built = simple_design.build(monthly_index[:400]) + extended = simple_design.extend(monthly_index[:400], 12) + + assert list(built.columns) == simple_design.column_names + assert list(extended.columns) == simple_design.column_names + assert (built.dtypes == np.float64).all() + assert (extended.dtypes == np.float64).all() + + def test_extend_rejects_nonpositive_steps(self, monthly_index, simple_design): + with pytest.raises(ValueError, match="steps must be >= 1"): + simple_design.extend(monthly_index, 0) + with pytest.raises(ValueError, match="steps must be >= 1"): + simple_design.future_index(monthly_index, -3) + + +# --------------------------------------------------------------------------- # +# B. Slicing invariance (and the deliberate exception) +# --------------------------------------------------------------------------- # + + +class TestSlicingInvariance: + """Calendar terms ignore where the sample starts; trends deliberately do not.""" + + @pytest.mark.parametrize( + ("term_id", "term"), + [ + ("fourier", Fourier(period=12, order=2)), + ("month", SeasonalDummies(season="month")), + ("level", BreakDummy(date="2000-01-01")), + ("pulse", BreakDummy(date="2000-01-01", kind="pulse")), + ], + ) + def test_calendar_terms_are_slice_invariant(self, monthly_index, term_id, term): + design = DeterministicDesign(terms=[term]) + cut = 60 + + pdt.assert_frame_equal(design.build(monthly_index).iloc[cut:], design.build(monthly_index[cut:])) + + def test_trend_shifts_by_the_exact_offset(self, monthly_index): + design = DeterministicDesign(terms=[Trend(degree=1, scale=1.0)]) + cut = 60 + + on_full = design.build(monthly_index).iloc[cut:]["trend"].to_numpy() + on_slice = design.build(monthly_index[cut:])["trend"].to_numpy() + + # The origin moved with the sample start, so the two differ by exactly + # the number of periods dropped — an affine shift, not a bug. + assert not np.array_equal(on_full, on_slice) + assert_allclose(on_full - on_slice, float(cut)) + + +# --------------------------------------------------------------------------- # +# C. Column-name contract +# --------------------------------------------------------------------------- # + + +class TestColumnNames: + def test_trend_names(self): + assert Trend(degree=1).column_names == ["trend"] + assert Trend(degree=3).column_names == ["trend", "trend_squared", "trend_cubed"] + + def test_fourier_names(self): + assert Fourier(period=12, order=2).column_names == [ + "sin(1,12)", + "cos(1,12)", + "sin(2,12)", + "cos(2,12)", + ] + assert Fourier(period=365.25, order=1).column_names == ["sin(1,365.25)", "cos(1,365.25)"] + + @pytest.mark.parametrize( + ("period", "expected"), + [(12, "12"), (12.0, "12"), (365.25, "365.25"), (4, "4"), (52.18, "52.18")], + ) + def test_format_period(self, period, expected): + assert _format_period(period) == expected + + def test_seasonal_names_drop_the_reference_level(self): + assert SeasonalDummies(season="quarter").column_names == ["quarter_2", "quarter_3", "quarter_4"] + assert SeasonalDummies(season="quarter", reference=3).column_names == [ + "quarter_1", + "quarter_2", + "quarter_4", + ] + assert SeasonalDummies(season="dayofweek").column_names == [f"dow_{i}" for i in range(1, 7)] + assert len(SeasonalDummies(season="month").column_names) == 11 + assert SeasonalDummies(season="month", drop_first=False).column_names == [f"month_{i}" for i in range(1, 13)] + + def test_break_names_use_the_resolved_timestamp(self): + assert BreakDummy(date="2000-01-01").column_names == ["level_2000-01-01"] + assert BreakDummy(date=pd.Timestamp("2000-03-15"), kind="pulse").column_names == ["pulse_2000-03-15"] + + def test_design_column_names_concatenate_in_term_order(self, simple_design): + assert simple_design.column_names == [ + "trend", + "sin(1,12)", + "cos(1,12)", + "sin(2,12)", + "cos(2,12)", + "quarter_2", + "quarter_3", + "quarter_4", + "level_2000-01-01", + ] + + def test_duplicate_column_names_are_rejected(self): + with pytest.raises(ValueError, match="Duplicate column name 'trend'"): + DeterministicDesign(terms=[Trend(degree=1), Trend(degree=1, scale=12.0)]) + + def test_empty_design_is_rejected(self): + with pytest.raises(ValueError, match="requires at least one term"): + DeterministicDesign(terms=[]) + + +# --------------------------------------------------------------------------- # +# D. Validation errors +# --------------------------------------------------------------------------- # + + +class TestConstructionValidation: + @pytest.mark.parametrize("degree", [0, 4]) + def test_trend_degree_bounds(self, degree): + with pytest.raises(ValueError, match="degree"): + Trend(degree=degree) + + def test_trend_scale_must_be_positive(self): + with pytest.raises(ValueError, match="scale"): + Trend(degree=1, scale=0.0) + + def test_fourier_nyquist_limit(self): + with pytest.raises(ValueError, match="Nyquist"): + Fourier(period=12, order=7) + + def test_fourier_period_must_exceed_one(self): + with pytest.raises(ValueError, match="period"): + Fourier(period=1, order=1) + + def test_seasonal_reference_must_be_a_level(self): + with pytest.raises(ValueError, match="not a valid month level"): + SeasonalDummies(season="month", reference=13) + + def test_seasonal_reference_requires_drop_first(self): + with pytest.raises(ValueError, match="only meaningful with drop_first=True"): + SeasonalDummies(season="quarter", drop_first=False, reference=2) + + def test_unknown_season_rejected(self): + with pytest.raises(ValueError, match="season"): + SeasonalDummies(season="dayofyear") + + +class TestBuildValidation: + def test_pulse_not_on_the_index_names_its_neighbours(self, monthly_index): + design = DeterministicDesign(terms=[BreakDummy(date="1990-01-15", kind="pulse")]) + with pytest.raises(ValueError, match=r"1990-01-01 \(before\) and 1990-02-01 \(after\)"): + design.build(monthly_index) + + def test_level_break_at_sample_start_is_rejected(self, monthly_index): + design = DeterministicDesign(terms=[BreakDummy(date="1980-01-01")]) + with pytest.raises(ValueError, match="collinear with the intercept"): + design.build(monthly_index) + + def test_level_break_after_sample_end_is_rejected(self, monthly_index): + design = DeterministicDesign(terms=[BreakDummy(date="2050-01-01")]) + with pytest.raises(ValueError, match="never occurs in-sample"): + design.build(monthly_index) + + def test_full_dummy_set_is_collinear_with_the_intercept(self, monthly_index): + design = DeterministicDesign(terms=[SeasonalDummies(season="month", drop_first=False)]) + with pytest.raises(ValueError, match="drop_first=False"): + design.build(monthly_index) + + def test_degenerate_top_harmonic_is_rejected(self, monthly_index): + design = DeterministicDesign(terms=[Fourier(period=12, order=6)]) + with pytest.raises(ValueError, match="is identically"): + design.build(monthly_index) + + def test_dummies_and_harmonics_of_the_same_cycle_clash(self, monthly_index): + design = DeterministicDesign(terms=[SeasonalDummies(season="month"), Fourier(period=12, order=2)]) + with pytest.raises(ValueError, match="describe the same cycle"): + design.build(monthly_index) + + def test_too_few_observations(self, monthly_index): + design = DeterministicDesign(terms=[SeasonalDummies(season="month")]) + with pytest.raises(ValueError, match="Too few observations"): + design.build(monthly_index[:8]) + + def test_index_must_be_a_datetime_index(self, simple_design): + with pytest.raises(TypeError, match="must be a pandas DatetimeIndex"): + simple_design.build(pd.RangeIndex(10)) + + def test_index_must_be_strictly_increasing(self, simple_design, monthly_index): + shuffled = monthly_index[::-1] + with pytest.raises(ValueError, match="strictly increasing"): + simple_design.build(shuffled) + + def test_empty_index_rejected(self, simple_design, monthly_index): + with pytest.raises(ValueError, match="must not be empty"): + simple_design.build(monthly_index[:0]) + + def test_break_date_of_an_unsupported_type_is_rejected(self): + with pytest.raises(ValueError, match="date"): + BreakDummy(date=12345) + + def test_a_term_whose_width_contradicts_its_names_is_caught(self, monthly_index): + class BadTerm: + """A custom term that promises two columns and delivers one.""" + + @property + def column_names(self): + return ["a", "b"] + + def build(self, index, origin, alias): + return np.zeros((len(index), 1)) + + design = DeterministicDesign(terms=[BadTerm()]) + with pytest.raises(ValueError, match=r"BadTerm.build returned shape"): + design.build(monthly_index) + + +# --------------------------------------------------------------------------- # +# E. Frequency resolution +# --------------------------------------------------------------------------- # + + +class TestFrequencyResolution: + def test_explicit_freq_wins_over_a_gappy_index(self, gappy_monthly_index): + design = DeterministicDesign(terms=[Trend(degree=1)], freq="MS") + assert design.build(gappy_monthly_index).shape == (len(gappy_monthly_index), 1) + + def test_index_freq_is_used_when_present(self, monthly_index): + assert monthly_index.freq is not None + design = DeterministicDesign(terms=[Trend(degree=1)]) + assert_allclose(design.build(monthly_index)["trend"].to_numpy()[:4], [0, 1, 2, 3]) + + def test_inference_accepted_when_it_regenerates_the_index(self, monthly_index): + stripped = pd.DatetimeIndex(list(monthly_index)) + assert stripped.freq is None + design = DeterministicDesign(terms=[Trend(degree=1)]) + assert_allclose(design.build(stripped)["trend"].to_numpy()[:4], [0, 1, 2, 3]) + + def test_false_positive_inference_is_rejected(self): + # pandas confidently infers WOM-1SAT here; it has no period equivalent, + # so the design refuses rather than silently anchoring to nonsense. + irregular = pd.DatetimeIndex(["2020-01-04", "2020-02-01", "2020-03-07"]) + assert pd.infer_freq(irregular) == "WOM-1SAT" + design = DeterministicDesign(terms=[Trend(degree=1)]) + with pytest.raises(ValueError, match="no pandas period equivalent"): + design.build(irregular) + + def test_inference_that_does_not_regenerate_the_index_is_rejected(self, monkeypatch, monthly_index): + # `pd.infer_freq` is confident on short irregular indices, so the + # candidate must reproduce the index before it is trusted. Forcing a + # wrong-but-valid answer exercises that guard directly. + stripped = pd.DatetimeIndex(list(monthly_index)) + monkeypatch.setattr(pd, "infer_freq", lambda index: "QS") + design = DeterministicDesign(terms=[Trend(degree=1)]) + + with pytest.raises(ValueError, match="does not reproduce it"): + design.build(stripped) + + def test_unresolvable_frequency_errors(self): + irregular = pd.DatetimeIndex(["2020-01-01", "2020-01-05", "2020-03-17", "2020-09-02"]) + design = DeterministicDesign(terms=[Trend(degree=1)]) + with pytest.raises(ValueError, match="Could not determine the sampling frequency"): + design.build(irregular) + + def test_business_day_frequency_names_the_alternative(self): + business = pd.date_range("2020-01-01", periods=60, freq="B") + design = DeterministicDesign(terms=[Trend(degree=1)]) + with pytest.raises(ValueError, match=r'Business-day frequencies.*freq="D"'): + design.build(business) + + @pytest.mark.parametrize("freq", ["MS", "ME", "QS", "QE", "YS", "YE", "D", "W", "h", "15D", "2h", "15min"]) + def test_extend_walks_the_sampling_offset(self, freq): + index = pd.date_range("2000-01-03", periods=60, freq=freq) + design = DeterministicDesign(terms=[Trend(degree=1)]) + + expected_index = pd.date_range(index[-1], periods=6, freq=freq)[1:] + extended = design.extend(index, 5) + + assert extended.index.equals(expected_index) + assert design.future_index(index, 5).equals(expected_index) + last = design.build(index)["trend"].to_numpy()[-1] + assert_allclose(extended["trend"].to_numpy(), last + np.arange(1.0, 6.0)) + + @pytest.mark.parametrize("freq", ["15D", "2h", "15min"]) + def test_multiplied_offsets_count_in_sampling_periods(self, freq): + # pandas stores 15D ordinals in days and 2h ordinals in hours. Elapsed + # time must still advance by one per observation, or `Fourier.period` + # would silently mean something other than "cycle length in sampling + # periods" — a 12-observation cycle on 2-hourly data is 24 hours. + index = pd.date_range("2000-01-03", periods=60, freq=freq) + design = DeterministicDesign(terms=[Trend(degree=1)]) + + assert_allclose(design.build(index)["trend"].to_numpy(), np.arange(60.0)) + + @pytest.mark.parametrize("freq", ["15D", "2h", "15min"]) + @pytest.mark.parametrize("h", [1, 3, 13]) + def test_continuation_holds_for_multiplied_offsets(self, freq, h): + index = pd.date_range("2000-01-03", periods=60, freq=freq) + design = DeterministicDesign( + terms=[Trend(degree=1, scale=12.0), Fourier(period=12, order=2)], + ) + cut = len(index) - 13 + + pdt.assert_frame_equal(design.build(index[: cut + h]).iloc[cut:], design.extend(index[:cut], h)) + + def test_extend_does_not_skip_a_period_off_anchor(self): + # `pd.date_range` rolls an off-anchor start forward, so the walk's + # first entry is already April here. Dropping it would forecast from + # May and silently lose a month. + index = pd.DatetimeIndex(["2000-01-01", "2000-02-01", "2000-03-15"]) + design = DeterministicDesign(terms=[Trend(degree=1)], freq="MS") + + future = design.future_index(index, 3) + + assert list(future.strftime("%Y-%m-%d")) == ["2000-04-01", "2000-05-01", "2000-06-01"] + assert design.extend(index, 3).index.equals(future) + # The trend keeps counting calendar months from the origin. + assert_allclose(design.extend(index, 3)["trend"].to_numpy(), [3.0, 4.0, 5.0]) + + def test_on_anchor_extend_is_unchanged(self): + index = pd.date_range("2000-01-01", periods=3, freq="MS") + design = DeterministicDesign(terms=[Trend(degree=1)]) + + assert list(design.future_index(index, 3).strftime("%Y-%m-%d")) == [ + "2000-04-01", + "2000-05-01", + "2000-06-01", + ] + + def test_future_index_override_is_honoured(self, monthly_index): + design = DeterministicDesign(terms=[Trend(degree=1)]) + override = pd.DatetimeIndex(["2030-01-01", "2030-06-01", "2031-01-01"]) + + extended = design.extend(monthly_index, 3, future_index=override) + + assert extended.index.equals(override) + # Anchored to the estimation origin, so elapsed months are absolute. + assert_allclose(extended["trend"].to_numpy(), [600.0, 605.0, 612.0]) + + def test_future_index_length_must_match_steps(self, monthly_index): + design = DeterministicDesign(terms=[Trend(degree=1)]) + override = pd.DatetimeIndex(["2030-01-01", "2030-02-01"]) + with pytest.raises(ValueError, match="future_index has length 2, but steps=3"): + design.extend(monthly_index, 3, future_index=override) + + +# --------------------------------------------------------------------------- # +# F. Gaps +# --------------------------------------------------------------------------- # + + +class TestGaps: + def test_trend_jumps_across_a_gap(self, monthly_index): + gappy = monthly_index.delete(range(3, 6)) + design = DeterministicDesign(terms=[Trend(degree=1)], freq="MS") + + trend = design.build(gappy)["trend"].to_numpy() + + assert_allclose(trend[:8], [0, 1, 2, 6, 7, 8, 9, 10]) + + def test_dummies_stay_calendar_correct_across_a_gap(self, gappy_monthly_index): + design = DeterministicDesign(terms=[SeasonalDummies(season="month")], freq="MS") + + frame = design.build(gappy_monthly_index) + + months = gappy_monthly_index.month + for position, month in enumerate(months[:60]): + row = frame.iloc[position] + if month == 1: # the dropped reference level + assert row.sum() == 0.0 + else: + assert row[f"month_{month}"] == 1.0 + assert row.sum() == 1.0 + + def test_continuation_holds_across_a_gap(self, gappy_monthly_index): + design = DeterministicDesign( + terms=[Trend(degree=1, scale=12.0), Fourier(period=12, order=2)], + freq="MS", + ) + cut = len(gappy_monthly_index) - 12 + + expected = design.build(gappy_monthly_index).iloc[cut:] + actual = design.extend(gappy_monthly_index[:cut], 12) + + assert_allclose(actual.to_numpy(), expected.to_numpy()) + + +# --------------------------------------------------------------------------- # +# G. Finiteness +# --------------------------------------------------------------------------- # + + +class TestFiniteness: + """Terms are total functions of a timestamp — they cannot emit NaN.""" + + @pytest.mark.parametrize(("freq_key", "terms", "h"), _continuation_cases()) + def test_no_missing_values_anywhere(self, freq_key, terms, h): + index = _index_for(freq_key) + design = DeterministicDesign(terms=terms) + cut = len(index) - _TAIL + + assert np.isfinite(design.build(index[:cut]).to_numpy()).all() + assert np.isfinite(design.extend(index[:cut], h).to_numpy()).all() + + +# --------------------------------------------------------------------------- # +# H. VARData round trip +# --------------------------------------------------------------------------- # + + +def _endog_frame(index: pd.DatetimeIndex, n_vars: int = 2) -> pd.DataFrame: + rng = np.random.default_rng(0) + return pd.DataFrame( + rng.standard_normal((len(index), n_vars)), + index=index, + columns=[f"y{i + 1}" for i in range(n_vars)], + ) + + +class TestVARDataRoundTrip: + def test_design_columns_become_exog_names(self, monthly_index, simple_design): + endog = _endog_frame(monthly_index) + frame = pd.concat([endog, simple_design.build(monthly_index)], axis=1) + + data = VARData.from_df(frame, endog=list(endog.columns), exog=simple_design.column_names) + + assert data.exog_names == simple_design.column_names + assert data.exog.shape == (len(monthly_index), len(simple_design.column_names)) + assert data.exog.flags.writeable is False + + def test_misaligned_index_surfaces_as_the_nan_invariant(self, monthly_index, simple_design): + endog = _endog_frame(monthly_index) + # The documented failure mode: the design was built on the index that + # survived a transform, but concatenated against the untrimmed endog. + # Deterministic terms cannot themselves emit NaN, so any NaN in the + # exog block is misalignment — and VARData catches it at construction. + frame = pd.concat([endog, simple_design.build(monthly_index[3:])], axis=1) + + with pytest.raises(ValueError, match="exog contains NaN or Inf values"): + VARData.from_df(frame, endog=list(endog.columns), exog=simple_design.column_names) + + +# --------------------------------------------------------------------------- # +# I. exog_future +# --------------------------------------------------------------------------- # + + +def _data_with_design(index, design, order=None): + endog = _endog_frame(index) + frame = pd.concat([endog, design.build(index)], axis=1) + columns = design.column_names if order is None else order + return VARData.from_df(frame, endog=list(endog.columns), exog=columns) + + +class TestExogFuture: + def test_shape_and_dtype(self, monthly_index, simple_design): + data = _data_with_design(monthly_index, simple_design) + + block = simple_design.exog_future(data, 6) + + assert block.shape == (6, len(simple_design.column_names)) + assert block.dtype == np.float64 + assert_allclose(block, simple_design.extend(monthly_index, 6).to_numpy()) + + def test_columns_are_reordered_to_match_exog_names(self, monthly_index): + design = DeterministicDesign(terms=[Trend(degree=1, scale=120.0), Fourier(period=12, order=2)]) + permuted = ["cos(1,12)", "trend", "sin(2,12)", "sin(1,12)", "cos(2,12)"] + data = _data_with_design(monthly_index, design, order=permuted) + assert data.exog_names == permuted + + block = design.exog_future(data, 4) + + assert_allclose(block, design.extend(monthly_index, 4)[permuted].to_numpy()) + # Positional forecasting would otherwise silently use the wrong column. + assert not np.allclose(block, design.extend(monthly_index, 4).to_numpy()) + + def test_name_mismatch_names_both_sets(self, monthly_index, simple_design): + data = _data_with_design(monthly_index, simple_design) + other = DeterministicDesign(terms=[Trend(degree=1, scale=120.0), Fourier(period=12, order=1)]) + + with pytest.raises(ValueError, match="does not match the fitted exogenous block"): + other.exog_future(data, 4) + + def test_data_without_exog_errors(self, monthly_index, simple_design): + endog = _endog_frame(monthly_index) + data = VARData.from_df(endog, endog=list(endog.columns)) + + with pytest.raises(ValueError, match="fitted without exogenous regressors"): + simple_design.exog_future(data, 4) + + def test_accepts_a_fitted_var(self, monthly_index, simple_design): + fitted = _fitted_with_exog(monthly_index, simple_design) + + from_fitted = simple_design.exog_future(fitted, 5) + from_data = simple_design.exog_future(fitted.data, 5) + + assert_allclose(from_fitted, from_data) + + def test_future_index_override_flows_through(self, monthly_index, simple_design): + data = _data_with_design(monthly_index, simple_design) + override = simple_design.future_index(monthly_index, 3) + + assert_allclose( + simple_design.exog_future(data, 3, future_index=override), + simple_design.exog_future(data, 3), + ) + + +# --------------------------------------------------------------------------- # +# J. Fast end-to-end (synthetic posterior, no MCMC) +# --------------------------------------------------------------------------- # + + +def _fitted_with_exog(index, design, n_lags: int = 1): + """A FittedVAR over a synthetic posterior that carries `B_exog`.""" + import arviz as az + import xarray as xr + + from impulso.fitted import FittedVAR + from impulso.volatility import Constant + + data = _data_with_design(index, design) + n_chains, n_draws, n_vars = 2, 20, 2 + n_exog = len(design.column_names) + rng = np.random.default_rng(11) + + L = np.broadcast_to(np.eye(n_vars) * 0.2, (n_chains, n_draws, n_vars, n_vars)).copy() + posterior = xr.Dataset({ + "B": xr.DataArray( + rng.standard_normal((n_chains, n_draws, n_vars, n_vars * n_lags)) * 0.2, + dims=["chain", "draw", "var", "coeff"], + ), + "intercept": xr.DataArray( + rng.standard_normal((n_chains, n_draws, n_vars)) * 0.01, + dims=["chain", "draw", "var"], + ), + "B_exog": xr.DataArray( + rng.standard_normal((n_chains, n_draws, n_vars, n_exog)), + dims=["chain", "draw", "var", "exog"], + coords={"exog": design.column_names}, + ), + "L": xr.DataArray(L, dims=["chain", "draw", "var1", "var2"]), + }) + return FittedVAR.model_construct( + idata=az.InferenceData(posterior=posterior), + n_lags=n_lags, + data=data, + var_names=data.endog_names, + volatility=Constant(), + pymc_model=None, + ) + + +class TestForecastIntegrationFast: + def test_forecast_consumes_the_generated_block(self, monthly_index, simple_design): + fitted = _fitted_with_exog(monthly_index, simple_design) + + block = simple_design.exog_future(fitted, 6) + forecast = fitted.forecast(steps=6, exog_future=block, include_shock_uncertainty=False) + + assert forecast.median().shape == (6, 2) + + zeroed = fitted.forecast( + steps=6, + exog_future=np.zeros_like(block), + include_shock_uncertainty=False, + ) + # If B_exog were ignored the two would coincide; the design must bite. + assert not np.allclose(forecast.median().to_numpy(), zeroed.median().to_numpy()) + + def test_posterior_exog_coord_matches_the_design(self, monthly_index, simple_design): + fitted = _fitted_with_exog(monthly_index, simple_design) + + coord = list(fitted.idata.posterior["B_exog"].coords["exog"].values) + + assert coord == simple_design.column_names + + +# --------------------------------------------------------------------------- # +# K. Estimator boundary +# --------------------------------------------------------------------------- # + + +class TestConjugateEstimatorBoundary: + def test_conjugate_var_rejects_a_deterministic_design(self, monthly_index, simple_design): + from impulso.conjugate import ConjugateVAR + from impulso.priors import NIWPrior + + data = _data_with_design(monthly_index, simple_design) + estimator = ConjugateVAR(lags=1, prior=NIWPrior(), draws=2, tune=0, seed=0) + + with pytest.raises(ValueError, match="does not support exogenous regressors"): + estimator.fit(data) + + +# --------------------------------------------------------------------------- # +# L. Slow integration — the onboarding recipe, end to end +# --------------------------------------------------------------------------- # + + +@pytest.mark.slow +def test_deterministic_design_end_to_end(): + """The documented recipe, exercised against real MCMC. + + The body of this test is reproduced verbatim in + `docs/how-to/deterministic-regressors.md`. + """ + import numpy as np + import pandas as pd + + from impulso import VAR, DeterministicDesign, Fourier, NUTSSampler, Trend, VARData + + # A short monthly two-variable sample standing in for climate anomalies. + rng = np.random.default_rng(7) + index = pd.date_range("2000-01-01", periods=120, freq="MS") + endog = pd.DataFrame( + rng.standard_normal((len(index), 2)).cumsum(axis=0) * 0.1, + index=index, + columns=["temperature", "precipitation"], + ) + + # One design, used for estimation and for forecasting. + design = DeterministicDesign( + terms=[Trend(degree=1, scale=120.0), Fourier(period=12, order=1)], + freq="MS", + ) + + frame = pd.concat([endog, design.build(index)], axis=1) + data = VARData.from_df(frame, endog=list(endog.columns), exog=design.column_names) + + fitted = VAR(lags=1).fit(data, sampler=NUTSSampler(draws=50, tune=50, chains=2, cores=1, random_seed=42)) + + # The posterior labels B_exog with the design's own column names. + assert list(fitted.idata.posterior["B_exog"].coords["exog"].values) == design.column_names + + forecast = fitted.forecast(steps=12, exog_future=design.exog_future(fitted, 12)) + + assert forecast.median().shape == (12, 2) diff --git a/tests/test_fitted.py b/tests/test_fitted.py index ccd90e4..a8645b6 100644 --- a/tests/test_fitted.py +++ b/tests/test_fitted.py @@ -500,3 +500,72 @@ def test_time_varying_sigma_broadcasts_over_the_time_axis(self): ratios = np.diagonal(actual, axis1=-2, axis2=-1) / np.diagonal(sigma, axis1=-2, axis2=-1) np.testing.assert_allclose(ratios[0, 0], 2.0) np.testing.assert_allclose(ratios[0, 1], 1.5) + + +class TestForecastExogValidation: + """`forecast` validates the future exogenous block before propagating it. + + Users generate this block programmatically now + (`DeterministicDesign.exog_future`), so a mis-shaped one must fail with a + named contract rather than an opaque einsum error deep in the loop. + """ + + @pytest.fixture + def fitted_with_exog(self, synthetic_idata_2v, var_data_2v): + import xarray as xr + + rng = np.random.default_rng(3) + exog = rng.standard_normal((var_data_2v.endog.shape[0], 2)) + data = VARData( + endog=var_data_2v.endog, + endog_names=list(var_data_2v.endog_names), + exog=exog, + exog_names=["x1", "x2"], + index=var_data_2v.index, + ) + idata = synthetic_idata_2v.copy() + idata.posterior["B_exog"] = xr.DataArray( + rng.standard_normal((2, 50, 2, 2)), + dims=["chain", "draw", "var", "exog"], + coords={"exog": ["x1", "x2"]}, + ) + return FittedVAR.model_construct( + idata=idata, + n_lags=1, + data=data, + var_names=["y1", "y2"], + ) + + def test_correct_shape_forecasts(self, fitted_with_exog): + result = fitted_with_exog.forecast( + steps=4, + exog_future=np.zeros((4, 2)), + include_shock_uncertainty=False, + ) + assert result.median().shape == (4, 2) + + def test_wrong_column_count_raises(self, fitted_with_exog): + with pytest.raises(ValueError, match=r"exog_future must have shape \(4, 2\), got \(4, 1\)"): + fitted_with_exog.forecast(steps=4, exog_future=np.zeros((4, 1))) + + def test_wrong_step_count_raises(self, fitted_with_exog): + with pytest.raises(ValueError, match=r"exog_future must have shape \(4, 2\), got \(6, 2\)"): + fitted_with_exog.forecast(steps=4, exog_future=np.zeros((6, 2))) + + def test_posterior_without_b_exog_raises(self, synthetic_idata_2v, var_data_2v): + rng = np.random.default_rng(4) + data = VARData( + endog=var_data_2v.endog, + endog_names=list(var_data_2v.endog_names), + exog=rng.standard_normal((var_data_2v.endog.shape[0], 1)), + exog_names=["x1"], + index=var_data_2v.index, + ) + fitted = FittedVAR.model_construct( + idata=synthetic_idata_2v, + n_lags=1, + data=data, + var_names=["y1", "y2"], + ) + with pytest.raises(ValueError, match="never consumed"): + fitted.forecast(steps=3, exog_future=np.zeros((3, 1))) diff --git a/tests/test_public_api.py b/tests/test_public_api.py index 7a5af57..8b4b447 100644 --- a/tests/test_public_api.py +++ b/tests/test_public_api.py @@ -197,3 +197,52 @@ def test_long_run_restriction_in_all(self): import impulso assert "LongRunRestriction" in impulso.__all__ + + +class TestDeterministicPublicAPI: + def test_design_and_terms_importable_from_impulso(self): + import impulso + from impulso.deterministic import ( + BreakDummy, + DeterministicDesign, + Fourier, + SeasonalDummies, + Trend, + ) + + assert impulso.DeterministicDesign is DeterministicDesign + assert impulso.Trend is Trend + assert impulso.Fourier is Fourier + assert impulso.SeasonalDummies is SeasonalDummies + assert impulso.BreakDummy is BreakDummy + + def test_deterministic_term_protocol_importable_from_impulso(self): + import impulso + from impulso.protocols import DeterministicTerm + + assert impulso.DeterministicTerm is DeterministicTerm + + def test_deterministic_names_in_all(self): + import impulso + + for name in ( + "BreakDummy", + "DeterministicDesign", + "DeterministicTerm", + "Fourier", + "SeasonalDummies", + "Trend", + ): + assert name in impulso.__all__ + + def test_terms_satisfy_the_protocol(self): + from impulso.deterministic import BreakDummy, Fourier, SeasonalDummies, Trend + from impulso.protocols import DeterministicTerm + + for term in ( + Trend(degree=1), + Fourier(period=12, order=1), + SeasonalDummies(season="month"), + BreakDummy(date="2000-01-01"), + ): + assert isinstance(term, DeterministicTerm)