API reference¶
The complete callable surface of tsecon, generated from the type stub (bindings/python/python/tsecon/__init__.pyi). Array arguments are float64 NumPy arrays (_ArrayLike = npt.NDArray[np.float64]; strided views are fine, plain lists and other dtypes are rejected at the boundary). Every function returns a plain dictionary, a NumPy array, or a Python scalar — no framework objects. Vector-valued keys are float64 NumPy arrays; matrix- and higher-rank-valued keys in the VAR/SVAR, Bayesian, multivariate-GARCH, panel and term-structure families (and the top-level results of var_irf, var_fevd and bvar_irf_draws) are nested Python lists — np.asarray(...) converts them; the docstring says which. For the why and when of each method, see the model cards and the guide.
192 functions.
diagnostics¶
acf¶
Autocorrelation function with Bartlett standard errors.
Returned keys: `acf`, `bartlett_se`.
Further arguments, with defaults: `nlags` (20), `adjusted` (False).
pacf¶
Partial autocorrelation function; method is "yw" or "ols".
Further arguments, with defaults: `nlags` (20).
ljung_box¶
Ljung-Box and Box-Pierce portmanteau tests for lags 1..=nlags.
Returned keys: `bp_pvalue`, `bp_stat`, `lags`, `lb_pvalue`, `lb_stat`.
jarque_bera¶
Jarque-Bera normality test (statistic, p_value, skewness, kurtosis, n).
Returned keys: `kurtosis`, `n`, `p_value`, `skewness`, `statistic`.
arch_lm¶
Engle's ARCH-LM test for conditional heteroskedasticity.
Returned keys: `df`, `nobs`, `p_value`, `statistic`.
Further arguments, with defaults: `nlags` (4).
unit roots / workflow¶
adf¶
def adf(
y: _ArrayLike,
regression: str = ...,
autolag: str | None = ...,
maxlag: int | None = ...,
) -> dict[str, Any]:
Augmented Dickey-Fuller test with MacKinnon p-values.
Returned keys: `crit`, `nobs`, `p_value`, `statistic`, `used_lag`.
Further arguments, with defaults: `regression` ("c"), `maxlag` (None).
kpss¶
KPSS stationarity test (null: stationary).
Returned keys: `lags`, `p_value`, `statistic`.
Further arguments, with defaults: `regression` ("c"), `nlags` (None).
check_stationarity¶
The ADF + KPSS confirmatory-quadrant workflow with a recommendation.
Returned keys: `adf_p_value`, `adf_statistic`, `alpha`,
`interpretation`, `kpss_p_value`, `kpss_statistic`, `quadrant`,
`recommendation`.
Further arguments, with defaults: `alpha` (0.05).
phillips_perron¶
def phillips_perron(
y: _ArrayLike,
regression: str = ...,
test_type: str = ...,
lags: int | None = ...,
) -> dict[str, Any]:
Phillips-Perron unit-root test (Z-tau/Z-alpha) with MacKinnon p-values.
Returned keys: `crit`, `lags`, `nobs`, `pvalue`, `stat`, `zalpha`,
`ztau`.
Further arguments, with defaults: `regression` ("c"), `test_type`
("tau"), `lags` (None).
dfgls¶
def dfgls(
y: _ArrayLike,
regression: str = ...,
lags: int | None = ...,
max_lags: int | None = ...,
method: str = ...,
) -> dict[str, Any]:
DF-GLS unit-root test (Elliott-Rothenberg-Stock 1996; null: unit root).
The ADF test run on a GLS-detrended series (quasi-differenced at the ERS
local alternative, cbar = -7.0 for "c", -13.5 for "ct") with no
deterministics in the test regression — near-optimal local power, the
recommended default over plain ADF. `regression`: "c" (constant, default)
or "ct" (constant + trend). `lags`: fixed lag count; None selects it by
`method` ("aic" default, "bic", "t-stat") on the OLS-detrended series
(Perron-Qu 2007) searching 0..=`max_lags` (default: Schwert's
ceil(12*(n/100)^(1/4)), capped at (n-1)/2 - 1). When `lags` is given,
`method`/`max_lags` are ignored (arch behavior). Returns `statistic`,
`p_value`, `used_lag`, `nobs` (= n - 1 - used_lag), `crit`
({"1%","5%","10%"}), `trend`. Statistic and selected lag match
arch.unitroot.DFGLS (< 1e-10); p-values/critical values are arch's DF-GLS
response surfaces (Sheppard's MacKinnon-style simulations, transcribed).
ng_perron¶
def ng_perron(
y: _ArrayLike,
trend: str = ...,
lags: int | str | None = ...,
max_lags: int | None = ...,
) -> dict[str, Any]:
Ng-Perron (2001) M unit-root tests (MZa, MZt, MSB, MPT; null: unit root).
GLS-detrends `y` through the same engine as `dfgls` (cbar = -7.0 for
"c", -13.5 for "ct"), selects the ADF lag by the paper's MAIC on the
detrended series (`lags=None` or `"maic"`, searching 0..=`max_lags`;
default Schwert's ceil(12*(n/100)^(1/4)) capped at (n-1)/2 - 1) or uses
a fixed integer `lags`, estimates the autoregressive spectral density at
frequency zero `s2_ar = sigma2_e / (1 - b(1))^2`, and forms the four M
statistics. All four reject the unit-root null when SMALL (below the
critical value); `mzt == mza * msb` exactly. No p-values: no published
response surface exists for the M tests, so compare each statistic
against its own critical values (Ng-Perron 2001 Table 1, asymptotic,
transcribed). Returns dict keys: `mza`, `mzt`, `msb`, `mpt`,
`used_lag`, `nobs` (= n - 1 - used_lag), `s2_ar`, `crit`
({"mza","mzt","msb","mpt"} each {"1%","5%","10%"}), `trend`. Prefer
this battery over `dfgls` under a suspected large negative MA root;
caveat (Perron-Qu 2007): on data far from the null MAIC drives the lag
to its maximum and power collapses — cap `max_lags` or fix `lags`
there.
phillips_ouliaris¶
def phillips_ouliaris(
y: _ArrayLike,
x: _ArrayLike,
trend: str = ...,
test_type: str = ...,
bandwidth: int | None = ...,
) -> dict[str, Any]:
Phillips-Ouliaris residual cointegration test (Zt/Za) with MacKinnon N-surfaces.
Returned keys: `crit`, `lags`, `n_vars`, `nobs`, `pvalue`, `stat`.
Further arguments, with defaults: `trend` ("c"), `test_type` ("Zt"),
`bandwidth` (None).
zivot_andrews¶
def zivot_andrews(
y: _ArrayLike,
regression: str = ...,
trim: float = ...,
max_lags: int | None = ...,
autolag: str | None = ...,
lags: int | None = ...,
) -> dict[str, Any]:
Zivot-Andrews unit-root test with one endogenous break.
Null: unit root with no break; alternative: stationary around one broken
deterministic component — `regression` "c" (intercept shift, default),
"t" (trend-slope shift), "ct" (both); the regression itself always has a
constant and a trend. The statistic is the minimum t on the lagged level
over candidate break dates inside the `trim` window (default 0.15, must
be in (0, 1/3] — 0 itself is unreachable); `break_index` is the last pre-break observation (the
shift begins at `break_index + 1`). Lag selection follows the
statsmodels/Baum single up-front convention on the "ct" base ADF:
`autolag` "aic" (default) / "bic" / "t-stat" capped at `max_lags`, or
`autolag=None` with `lags` fixed (both None: int(12*(n/100)**0.25)).
Pass either `lags` or `autolag`, not both. P-values and critical values
interpolate the statsmodels-simulated null table. Returns dict keys:
`stat`, `pvalue`, `crit` {"1%","5%","10%"}, `break_index`, `lags`,
`nobs`, `trim`, `regression`. Matches statsmodels `zivot_andrews`.
ndiffs¶
def ndiffs(
y: _ArrayLike, test: str = ..., alpha: float = ..., max_d: int = ...
) -> dict[str, Any]:
How many differences a series needs, with the per-order test evidence.
Returned keys: `alpha`, `d`, `interpretation`, `max_d`, `steps`, `stop`,
`test`.
Further arguments, with defaults: `alpha` (0.05), `max_d` (2).
nsdiffs¶
How many SEASONAL differences a series needs (Hyndman-Khandakar rule).
D += 1 while the STL seasonal strength is >= 0.64, capped at max_d
(the forecast::nsdiffs test="seas" rule; alpha is validated but unused
by this threshold rule, as in forecast). Returns `d`, `period`,
`threshold`, `alpha`, `max_d`, `stop`, per-order `steps`, and an
`interpretation`.
box_cox_lambda¶
def box_cox_lambda(
y: _ArrayLike,
method: str = ...,
bounds: tuple[float, float] | None = ...,
period: int | None = ...,
) -> dict[str, Any]:
Variance-stabilising Box-Cox lambda (MLE or Guerrero) with its objective.
`bounds` (None = (-2.0, 2.0)) are hard bounds on lambda; an optimum on a
bound is reported via `at_bound`.
Returned keys: `at_bound`, `interpretation`, `lambda`, `loglik_at_one`,
`loglik_at_zero`, `lower`, `lr_vs_one`, `lr_vs_zero`, `method`, `n`,
`objective`, `period`, `upper`.
Further arguments, with defaults: `method` ("mle"), `period` (None).
check_series¶
def check_series(
data: npt.ArrayLike,
seasonal_period: int | None = ...,
lags: int | None = ...,
alpha: float = ...,
max_breaks: int = ...,
trim: float = ...,
) -> dict[str, Any]:
One-call diagnostic battery with model recommendations (the Module 01 flagship).
Pure Python over the compiled tests, so plain lists are coerced. 1D input
runs descriptives, outliers, the ADF+KPSS quadrant, Ljung-Box/ACF/PACF,
ARCH-LM, Jarque-Bera, a sup-F/Bai-Perron mean-shift scan, GPH long memory,
and seasonality evidence; 2D (n, k) input runs per-series integration,
Johansen, and VAR lag selection with a stability check. Evidence is
reported in families with the multiple-testing arithmetic shown — never
silently corrected — and the report ends in an ordered `recommendations`
list routing to concrete tsecon calls. JSON-serializable throughout.
`lags` is shape-dependent: the Ljung-Box horizon for 1D input (default
min(10, n//5)), the VAR lag-search cap for 2D input (default 8). `alpha`
must lie in (0.01, 0.10] — the compiled KPSS p-value is clamped to that
range. `seasonal_period` must be an integer >= 2 with at least two full
cycles in sample.
Returned keys: `alpha`, `analysis_scale`, `arch_effects`, `breaks`,
`descriptives`, `kind`, `long_memory`, `multiple_testing`, `n`,
`normality`, `outliers`, `recommendations`, `seasonality`,
`serial_correlation`, `stationarity`, `tests_run`.
Further arguments, with defaults: `max_breaks` (5), `trim` (0.15).
summarize¶
Render any tsecon output as a readable results object (opt-in).
`print(tsecon.summarize(tsecon.adf(y)))` works for every function: a plain
dict becomes a generic `tsecon.results.Result` with an aligned `.summary()`,
while a bespoke `tsecon.results.*` object is returned unchanged. Additive —
the returned object is a `dict` subclass, so the plain-dict contract holds.
`wrap="generic"` forces the structural dump even on a bespoke object.
robust inference¶
long_run_variance¶
Kernel long-run variance of a series (demeaned internally).
Further arguments, with defaults: `kernel` ("bartlett"), `bandwidth`
(None).
ols¶
def ols(
y: _ArrayLike,
x: _ArrayLike,
se_type: str = ...,
maxlags: int | None = ...,
use_correction: bool = ...,
) -> dict[str, Any]:
OLS with nonrobust / HC0 / HC1 / HC2 / HC3 / HAC standard errors.
The leverage-corrected hc2/hc3 are what matter in small samples with
influential points; hc1's n/(n-k) factor barely moves. HC is
heteroskedasticity-robust only -- under serial correlation use "hac".
HAC matches statsmodels cov_type="HAC" when use_correction is matched;
the DEFAULTS differ deliberately (tsecon True, statsmodels False), so
pass use_correction=False to reproduce a default statsmodels call.
Returned keys: `bse`, `params`, `se_type`, `tvalues`.
Further arguments, with defaults: `se_type` ("hac"), `maxlags` (None).
bootstrap¶
bootstrap_indices¶
def bootstrap_indices(
n: int,
scheme: str = ...,
seed: int = ...,
block_length: int | None = ...,
p: float | None = ...,
) -> npt.NDArray[np.uint64]:
Bootstrap resampling indices (iid/moving/circular/stationary).
Further arguments, with defaults: `scheme` ("stationary"), `seed` (0),
`block_length` (None), `p` (None).
optimal_block_length¶
Politis-White (2004) automatic block length (stationary, circular).
Returned keys: `circular`, `stationary`.
philox_uniforms¶
Uniform draws from the Philox stream; bit-identical to NumPy.
state space¶
local_level_smooth¶
Exact-diffuse local-level Kalman filter + smoother (NaN = missing).
Keys: `loglik`, `filtered_state`/`filtered_state_var`,
`smoothed_state`/`smoothed_state_var` (each length n), and `d_diffuse`
(the number of initial observations spent in the exact-diffuse period).
ar_loglik¶
def ar_loglik(
y: _ArrayLike, coeffs: Sequence[float], sigma2: float, intercept: float = ...
) -> float:
Exact Gaussian log-likelihood of an AR(p) at fixed parameters.
Further arguments, with defaults: `intercept` (0.0).
NaN entries in `y` are treated as missing observations: the Kalman
filter skips their update and the log-likelihood sums the remaining
innovations (with `coeffs=[0]` it equals the log-likelihood of the
series with those entries deleted). Infinite entries are rejected.
ARIMA¶
arima_fit¶
def arima_fit(
y: _ArrayLike,
p: int = ...,
d: int = ...,
q: int = ...,
seasonal: tuple[int, int, int, int] | None = ...,
constant: bool = ...,
forecast_steps: int = ...,
conf_alpha: float | None = ...,
drift_uncertainty: bool = ...,
) -> dict[str, Any]:
Exact-MLE ARIMA(p,d,q) fit, with optional forecast + conf_alpha bands.
seasonal=(P, D, Q, s) fits the multiplicative SARIMA(p,d,q)(P,D,Q)_s —
the airline model is seasonal=(0, 1, 1, 12) on the logged series with
p=0, d=1, q=1, constant=False. Seasonal parameters are named
statsmodels-style (ar.S.L12, ma.S.L12); differencing (regular and
seasonal) is simple differencing, losing d + D*s observations.
Forecast standard errors treat parameters as known by default (the
statsmodels get_forecast convention). With d >= 1 and constant=True that
omits the estimated drift's own uncertainty, which grows like h^2 and
measurably under-covers: 90.2% at h=24, T=60 against a nominal 95%. Pass
drift_uncertainty=True to add it (94.5% on the same design).
Also returns bse / param_cov from the observed information, or None with
cov_ok=False when that matrix is too ill-conditioned to invert honestly.
converged reports the optimizer's convergence certificate (False = best
point found, not a certified optimum). boundary flags, per parameter,
AR/MA blocks whose fitted polynomial has a root within 0.1% of the unit
circle (auto_arima's admissibility epsilon): no classical SE exists
there — those bse entries are NaN with se_valid False, and
boundary_note (str | None) explains (an MA root at the unit circle is
the classic over-differencing signature). Interior bse still come from
the full-vector observed information, which a boundary degrades — treat
as approximate; reduced-Hessian interior SEs are a stated follow-up.
Returned keys: `aic`, `bic`, `boundary`, `boundary_note`, `bse`,
`conf_alpha`, `converged`, `cov_ok`, `drift_uncertainty`,
`forecast_lower`, `forecast_mean`, `forecast_se`, `forecast_upper`,
`loglik`, `param_cov`, `param_names`, `params`, `residuals`, `se_valid`.
Further arguments, with defaults: `forecast_steps` (0).
auto_arima¶
def auto_arima(
y: _ArrayLike,
seasonal_period: int = ...,
ic: str = ...,
stepwise: bool = ...,
max_p: int = ...,
max_q: int = ...,
max_P: int = ...,
max_Q: int = ...,
max_order: int = ...,
max_d: int = ...,
max_D: int = ...,
d: int | None = ...,
D: int | None = ...,
alpha: float = ...,
forecast_steps: int = ...,
conf_alpha: float | None = ...,
) -> dict[str, Any]:
Automatic ARIMA order selection (Hyndman-Khandakar 2008 stepwise).
D from the STL seasonal-strength rule (nsdiffs, when
seasonal_period >= 2), d from successive KPSS tests (ndiffs) on the
seasonally differenced series, then a stepwise search over
(p, q, P, Q, constant) minimizing `ic` ("aicc" default, "aic",
"bic") at those fixed differencing orders; stepwise=False fits the
exhaustive grid subject to max_order instead (like R, max_order
binds only the grid). Near-unit-root fits are recorded but never
selected; failed fits steer the search rather than aborting it.
Every candidate is fit by the exact-MLE engine behind arima_fit, so
the search is deterministic. No exogenous regressors in this slice.
Returns the arima_fit result dict for the selected model plus:
`order`, `seasonal_order`, `constant`, `converged`, `ic`,
`ic_value`, `aicc`, `stepwise`, `n_models`, `budget_exhausted`,
`trace` (every candidate tried, with its criterion and status),
`d_test` / `D_test` (the full ndiffs / nsdiffs evidence, None when
fixed or not applicable), and `interpretation`. Honest grading:
candidate fits are statsmodels-pinned; the selection loop itself is
graded by Monte-Carlo order recovery (rates in the model card), not
R/pmdarima parity.
Returned keys: `D_test`, `aic`, `aicc`, `bic`, `boundary`,
`boundary_note`, `bse`, `budget_exhausted`, `constant`, `converged`,
`cov_ok`, `d_test`, `ic`, `ic_value`, `interpretation`, `loglik`,
`n_models`, `order`, `param_cov`, `param_names`, `params`, `residuals`,
`se_valid`, `seasonal_order`, `stepwise`, `trace`.
Further arguments, with defaults: `max_p` (5), `max_q` (5), `max_P` (2),
`max_Q` (2), `max_d` (2), `max_D` (1), `alpha` (0.05), `forecast_steps`
(0), `conf_alpha` (None).
GARCH¶
garch_fit¶
def garch_fit(
y: _ArrayLike,
vol: str = ...,
mean: str = ...,
dist: str = ...,
p: int = ...,
o: int | None = ...,
q: int = ...,
forecast_horizon: int = ...,
) -> dict[str, Any]:
GARCH/GJR/EGARCH QMLE with MLE and Bollerslev-Wooldridge robust SEs.
Filter timing (matches `arch`): conditional_volatility[t] is the
one-step-ahead volatility FOR period t, formed from information through
t-1 (sigma2_t is built from eps_{t-1} and sigma2_{t-1}); the post-sample
continuation of that step is `variance_forecast`.
Named access: `params_named` is dict(zip(param_names, params)) — use
fit["params_named"]["omega"] on the raw dict (fit["omega"] is a
deliberate KeyError; the params/param_names parallel arrays stay the
positional source of truth).
`o` is the asymmetry (threshold) order and only `vol="gjr"`/`"egarch"`
have an asymmetry term, so `o > 0` with `vol="garch"` raises instead of
being silently discarded — the porting trap this guards: in the `arch`
package, `arch_model(y, p=1, o=1, q=1)` silently switches the model to
GJR-GARCH, while tsecon keeps that model choice explicit (pass
`vol="gjr"`). `o=None` (default) means no asymmetry term under
`vol="garch"` and one asymmetry lag under `vol="gjr"`/`"egarch"`.
Boundary fits (a coefficient at its sign constraint, persistence at 1)
carry per-parameter `se_valid`/`boundary` flags and a `boundary_note`:
boundary parameters have NaN standard errors (no classical asymptotics
exist there), interior parameters keep finite ones. `converged` reports
the optimizer's own verdict.
EGARCH has no closed-form multi-step variance forecast: with
`vol="egarch"` only `forecast_horizon` 0 or 1 is accepted, and
`forecast_horizon >= 2` raises a teaching ValueError (GARCH and GJR
forecast analytically at every horizon).
Returned keys: `aic`, `bic`, `boundary`, `boundary_note`,
`conditional_volatility`, `converged`, `loglik`, `param_names`,
`params`, `params_named`, `se_mle`, `se_robust`, `se_valid`,
`std_residuals`, `variance_forecast`.
Further arguments, with defaults: `mean` ("zero"), `dist` ("normal").
VAR¶
var_fit¶
Fit a VAR(p) by OLS; params, sigma_u, ICs, residuals, and stability.
Read `is_stable` for the stability verdict. `min_root`/`max_root` are the
smallest/largest moduli of the reciprocal characteristic roots — stable iff
`min_root > 1`, so `max_root` alone is not a verdict.
Keys: params, sigma_u, llf, aic, bic, hqic, resid ((T, k) nested list —
the OLS residuals over the effective sample, row t = observation lags+t),
fitted ((T, k) — the one-step fitted values, defined as
data[lags:] - resid, i.e. the OLS projection Z @ B; fitted + resid
reproduces data[lags:] exactly), nobs (T = len(data) - lags), df_resid
(T - m regressors per equation), max_root, min_root, is_stable.
Further arguments, with defaults: `trend` ("c").
var_irf¶
def var_irf(
data: _ArrayLike,
lags: int = ...,
horizon: int = ...,
orth: bool = ...,
trend: str = ...,
cumulative: bool = ...,
) -> list[list[list[float]]]:
Impulse responses [h][response][shock]; cumulative gives running sums.
Point path only. For frequentist confidence bands use `var_irf_bands`.
Further arguments, with defaults: `lags` (2), `horizon` (10), `orth`
(True), `trend` ("c").
var_irf_bands¶
def var_irf_bands(
data: _ArrayLike,
lags: int = ...,
horizon: int = ...,
orth: bool = ...,
method: str = ...,
alpha: float = ...,
cumulative: bool = ...,
n_boot: int = ...,
seed: int = ...,
trend: str = ...,
bias_correct: bool = ...,
band: str = ...,
band_scope: str = ...,
band_seed: int = ...,
band_n_sim: int = ...,
) -> dict[str, Any]:
Frequentist confidence bands on VAR impulse responses — the banded companion to var_irf.
Returns a dict with `point`/`se`/`lower`/`upper`, each `[h][response][shock]`
(same layout as `var_irf`), plus echoed `method`/`alpha`/`n_boot`/`band`.
`method`: "asymptotic" (Lütkepohl 1990 delta-method SEs, Wald bands
`point ± z_{1-alpha/2}·se`; `n_boot` is `None`) or "bootstrap" (residual
Efron/Kilian bootstrap, percentile bands, optional Kilian 1998
`bias_correct`). `orth` and `cumulative` behave exactly as in `var_irf`.
**Simultaneous bands.** `lower`/`upper` are POINTWISE whatever you pass:
each covers one `(horizon, response, shock)` cell and promises nothing
about the path as a whole. Set `band` to `"sup-t"`, `"sidak"` or
`"bonferroni"` to also get `sim_lower`/`sim_upper` — the same `point` and
the same `se` with a larger multiplier — plus `critical_value` (a k x k
grid), `pointwise_critical_value`, `band_scope`, `n_cells` (K) and
`n_cells_used`. `band="pointwise"` is the default and adds nothing.
Simultaneous **over what** is your choice and is always reported back:
`band_scope="horizon"` (default; `K = horizon+1`, one family per
response-shock pair — the object the coverage audit measured),
`"shock"` (`K = k(horizon+1)`) or `"all"` (`K = k²(horizon+1)`). Every cell
added to a family widens the band for every other cell in it.
**What it fixes and what it does not.** Audit design, nominal 90%, T=500,
h=0..12, 3000 replications, asymptotic branch: the pointwise band contained
the whole path in 70.4% ± 0.8 of samples, the sup-t band in 84.8% ± 0.7
(crate Monte Carlo, T=500, 3000 reps). The published harness measures the
same shape on its own BASE design: 71.7% ± 1.4 pointwise, 85.2% ± 1.1
sup-t. The sup-t rate **does not reach nominal**, and the residual is not
multiplicity — the pointwise band it is built from covers only about 91%
marginally at h=0 falling to 85.3% at h=12 against nominal 90%. sup-t fixes
multiplicity exactly and inherits everything else, so what is left needs a
better standard error, not a bigger multiplier.
**Shape, on the bootstrap branch.** `lower`/`upper` there are Efron
*percentile* bounds and pick up bootstrap skewness; the simultaneous band is
the symmetric `point ± c·se`. These are different shapes of interval, so
`sim_lower` is **not** guaranteed to sit below `lower` cell by cell. What it
is guaranteed to contain is the symmetric pointwise band
`point ± pointwise_critical_value·se` — the like-for-like comparator, in
which only the multiplier differs.
`band_seed`/`band_n_sim` drive the Gaussian simulation behind `"sup-t"` on
the asymptotic branch only, where the band is a pure function of
`band_seed`. On the bootstrap branch sup-t reads its quantile off the
bootstrap replications, so `seed` alone reproduces it (use `n_boot` ≥ 999);
Šidák and Bonferroni are closed forms in K and need neither. Method:
Montiel Olea and Plagborg-Møller.
Further arguments, with defaults: `lags` (2), `trend` ("c").
var_fevd¶
def var_fevd(
data: _ArrayLike, lags: int = ..., horizon: int = ..., trend: str = ...
) -> list[list[list[float]]]:
Forecast-error variance decomposition [h][variable][shock].
Horizon-first — `fevd[h][i][j]` is the share of variable i's (h+1)-step
forecast-error variance attributed to shock j — the same axis order as
`var_irf` and `structural_fevd`, with `horizon` outer entries; each
`fevd[h][i]` sums to 1 across shocks. (Before 0.6.0 the emitted list was
variable-major, contradicting this docstring; at k == horizon the two
layouts silently alias.) statsmodels stores the same numbers
variable-major: `np.transpose(res.fevd(horizon).decomp, (1, 0, 2))`
equals this output.
Further arguments, with defaults: `lags` (2), `trend` ("c").
var_forecast¶
def var_forecast(
data: _ArrayLike,
lags: int = ...,
steps: int = ...,
alpha: float = ...,
trend: str = ...,
band: str = ...,
band_scope: str = ...,
band_seed: int = ...,
band_n_sim: int = ...,
) -> dict[str, Any]:
Iterated VAR point forecasts with (1-alpha) intervals.
**Simultaneous bands.** `lower`/`upper` are MARGINAL whatever you pass: each
covers one `(horizon, series)` cell. Read as a statement about a whole fan
chart they are the worst offender in the library — the interval-coverage
audit, nominal 95% at T=100 over 12 horizons x 2 series, 6000 replications,
measured the marginal bands containing every cell at once in 41.2% ± 0.6 of
samples, and still only 48.1% at T=800. That is multiplicity, not a small
sample.
Set `band` to `"sup-t"`, `"sidak"` or `"bonferroni"` to also get `se`,
`sim_lower`/`sim_upper` (same `point`, same `se`, larger multiplier),
`critical_value` (one per series), `pointwise_critical_value`, `band_scope`,
`n_cells` (K) and `n_cells_used`. `band="pointwise"` is the default and adds
nothing. `band_scope="all"` (default) is `K = steps*k`, every horizon of
every series as one statement — the object the audit measured; `"horizon"`
is `K = steps`, one family per series.
**What it fixes and what it does not.** On that design the sup-t joint rate
was 90.5% ± 0.4 against a nominal 95%. It **does not reach nominal**, and
the residual is not multiplicity: these intervals are a plug-in treating the
coefficients as known, so their measured *marginal* rate is 93.3%, not 95%.
sup-t fixes multiplicity exactly and inherits that approximation unchanged.
`band_seed`/`band_n_sim` drive the Gaussian simulation behind `"sup-t"`, so
that band is a pure function of `band_seed`; the closed forms use neither.
Method: Montiel Olea and Plagborg-Møller.
Further arguments, with defaults: `lags` (2), `trend` ("c").
var_granger¶
def var_granger(
data: _ArrayLike,
caused: Sequence[int],
causing: Sequence[int],
lags: int = ...,
trend: str = ...,
) -> dict[str, Any]:
Granger-causality F test (matches statsmodels test_causality).
Returned keys: `df_den`, `df_num`, `p_value`, `statistic`.
Further arguments, with defaults: `lags` (2), `trend` ("c").
Bayesian VAR¶
bvar_fit¶
def bvar_fit(
data: _ArrayLike,
lags: int = ...,
lambda0: float = ...,
lambda1: float = ...,
lambda3: float = ...,
delta: float = ...,
scale_ar: int = ...,
) -> dict[str, Any]:
Minnesota-NIW conjugate BVAR posterior + log marginal likelihood. scale_ar sets the lag order of the AR residual-variance scale regressions (4 = default; 1 = the GLP 2015 convention).
Keys: posterior_mean_coefs (Bbar, k x K with k = 1 + p*K), sigma_posterior_mean,
log_marginal_likelihood, and the full NIW posterior: omega_bar (k x k),
s_bar (K x K), v_bar. Convention: vec(B)|Sigma,Y ~ N(vec(Bbar),
Sigma (x) Obar) with column-stacked vec (B.flatten(order="F"); Kronecker
order np.kron(sigma, omega_bar)), Sigma|Y ~ InvWishart(s_bar, v_bar).
Marginal coefficient posterior sd (v_bar > K + 1):
sd = np.sqrt(np.outer(np.diag(omega_bar), np.diag(s_bar)) / (v_bar - K - 1)),
a (k, K) array aligned with posterior_mean_coefs.
Further arguments, with defaults: `lags` (2), `lambda0` (100.0),
`lambda1` (0.2), `lambda3` (1.0), `delta` (0.0).
bvar_irf_draws¶
def bvar_irf_draws(
data: _ArrayLike,
lags: int = ...,
horizon: int = ...,
n_draws: int = ...,
seed: int = ...,
lambda0: float = ...,
lambda1: float = ...,
lambda3: float = ...,
delta: float = ...,
cumulative: bool = ...,
scale_ar: int = ...,
) -> list[list[list[list[float]]]]:
Posterior Cholesky-IRF draws [draw][h][variable][shock] for credible bands.
Further arguments, with defaults: `lags` (2), `horizon` (16), `n_draws`
(500), `seed` (0), `lambda0` (100.0), `lambda1` (0.2), `lambda3` (1.0),
`delta` (0.0), `cumulative` (False), `scale_ar` (4).
bvar_hierarchical¶
def bvar_hierarchical(
data: _ArrayLike,
lags: int = ...,
delta: float = ...,
lambda0: float = ...,
lambda3: float = ...,
lambda1_init: float = ...,
lambda1_lo: float = ...,
lambda1_hi: float = ...,
optimize: str = ...,
hyperprior: str = ...,
n_grid: int = ...,
max_iter: int = ...,
tol: float = ...,
scale_ar: int = ...,
) -> dict[str, Any]:
Empirical-Bayes Minnesota-BVAR: pick lambda1 by maximizing the marginal likelihood (Giannone-Lenza-Primiceri 2015). Default hyperprior="glp" (MAP-II under the GLP Gamma hyperprior) — pure ML-II (hyperprior="none") collapses lambda1 to the search-box floor on ~a fifth to a quarter of in-model datasets (audit round 6); a lambda1_opt at the box bottom is a red flag, not a selection. scale_ar=1 switches the prior's residual-scale regressions to GLP's own AR(1) convention (default 4).
Returned keys: `converged`, `grid_lambda1`, `grid_log_ml`,
`lambda1_fixed_log_ml`, `lambda1_opt`, `lambda3_opt`,
`log_marginal_likelihood`, `log_posterior`, `n_evals`,
`posterior_mean_coefs`, `sigma_posterior_mean`.
Further arguments, with defaults: `lags` (2), `delta` (0.0), `lambda0`
(100.0), `lambda3` (1.0), `lambda1_init` (0.2), `lambda1_lo` (0.0001),
`lambda1_hi` (10.0), `optimize` ("lambda1"), `n_grid` (25), `max_iter`
(200), `tol` (1e-08).
bvar_ssvs¶
def bvar_ssvs(
data: _ArrayLike,
lags: int = ...,
n_draws: int = ...,
burn: int = ...,
seed: int = ...,
c0: float = ...,
c1: float = ...,
prior_inclusion: float = ...,
ssvs_cov: bool = ...,
kappa0: float | None = ...,
kappa1: float | None = ...,
prior_inclusion_cov: float = ...,
gamma_a: float = ...,
gamma_b: float | None = ...,
horizon: int = ...,
thin: int = ...,
n_chains: int = ...,
) -> dict[str, Any]:
SSVS-BVAR (George-Sun-Ni 2008): spike-and-slab stochastic-search selection of VAR (and error-precision) restrictions by Gibbs; posterior inclusion probabilities, coef/Sigma means, and orthogonalized IRF draws. Default hyperpriors are unit-adaptive (None = scale by the per-equation OLS residual variance); explicit gamma_b/kappa0/kappa1 floats pin absolute prior scales.
Returned keys: `coef_mean`, `diagnostics`, `inclusion_prob`,
`irf_draws`, `sigma_mean`.
Further arguments, with defaults: `lags` (2), `n_draws` (10000), `burn`
(2000), `seed` (0), `c0` (0.1), `c1` (10.0), `prior_inclusion` (0.5),
`ssvs_cov` (False), `prior_inclusion_cov` (0.5), `gamma_a` (0.01),
`horizon` (16), `thin` (1), `n_chains` (1).
mcmc_diagnostics¶
Rank-normalized split R-hat and bulk/tail ESS (ArviZ-exact).
Returned keys: `ess_bulk`, `ess_tail`, `rhat`.
filters¶
hp_filter¶
Hodrick-Prescott filter (O(n)); one_sided=True for the real-time variant.
Keys: `trend`, `cycle` (= y - trend), `first_index` (0: full sample).
Further arguments, with defaults: `lamb` (1600.0).
bk_filter¶
def bk_filter(
y: _ArrayLike, low: float = ..., high: float = ..., k: int = ...
) -> dict[str, Any]:
Baxter-King band-pass filter (loses k observations at each end).
Keys: `cycle` (length n - 2k) and `first_index` (= k); no trend is
returned.
Further arguments, with defaults: `low` (6.0), `high` (32.0).
cf_filter¶
def cf_filter(
y: _ArrayLike, low: float = ..., high: float = ..., drift: bool = ...
) -> dict[str, Any]:
Christiano-Fitzgerald asymmetric band-pass filter (full sample).
Keys: `trend`, `cycle`, `first_index` (0: full sample).
Further arguments, with defaults: `low` (6.0), `high` (32.0), `drift`
(True).
hamilton_filter¶
def hamilton_filter(
y: _ArrayLike,
h: int = ...,
p: int = ...,
method: str = ...,
se: str | None = ...,
maxlags: int | None = ...,
use_correction: bool | None = ...,
) -> dict[str, Any]:
Hamilton (2018) regression filter — the modern HP alternative.
Frequency defaults (h, p): quarterly (8, 4), monthly (24, 12), annual
(2, 1). `method="random_walk"` is the short-sample variant `cycle =
y_t - y_{t-h}` (no regression). `se="hac"` adds Newey-West standard
errors on `beta` via the shared HAC engine with the h-overlap default
`maxlags = h` (the overlapping h-step residuals are MA(h-1) by
construction); `se="nonrobust"` the classical comparison point.
`maxlags` and `use_correction` parameterize the HAC covariance only,
so passing either explicitly under anything but `se="hac"` raises
(use_correction=None, the default, means True where HAC applies).
Returns `trend`, `cycle`, `first_index`, `beta` (regression method),
plus `bse`, `tvalues`, `se_type`, `maxlags`, `use_correction` when
`se` is requested. Defaults are bit-identical to earlier releases.
bn_filter¶
def bn_filter(
y: _ArrayLike,
p: int = ...,
delta: float | None = ...,
demean: str = ...,
d0: float | None = ...,
dt: float | None = ...,
) -> dict[str, Any]:
Kamber-Morley-Wong (2018) BN filter: Beveridge-Nelson output gap
with the signal-to-noise ratio pinned (delta=None: their
amplitude-to-noise selection on the grid (d0, dt); d0/dt default to
0.01/0.0005 when None and parameterize only that grid, so passing
either explicitly with a fixed delta= raises); demean "sm"
(sample mean, the baseline) or "nd" (no drift removal). Returns
trend, cycle, first_index (=1), delta, ar, cycle_se,
amplitude_to_noise, drift. Reference-run-validated against the
authors' R replication code.
Further arguments, with defaults: `p` (12).
bn_decomposition¶
def bn_decomposition(
y: _ArrayLike,
p: int = ...,
q: int = ...,
ar: _ArrayLike | None = ...,
ma: _ArrayLike | None = ...,
drift: float | None = ...,
) -> dict[str, Any]:
Classic Beveridge-Nelson (1981) decomposition from an
ARIMA(p, 1, q) — fit by the library's exact MLE (default: the
Morley-Nelson-Zivot p=2, q=2), or at fixed coefficients when any of
ar/ma/drift is passed (p/q are ignored on that
fixed-coefficient path — the decomposition is computed exactly at the
coefficients you supply). Returns trend, cycle,
innovations, first_index (=1), long_run_multiplier
(psi(1) = theta(1)/phi(1)), drift, ar, ma, mode, and for
the fit path sigma2, loglik, aic, bic, converged;
trend + cycle reconstructs y[1:] to within a final rounding.
stl¶
def stl(
y: _ArrayLike,
period: int,
seasonal: int = ...,
trend: int | None = ...,
low_pass: int | None = ...,
seasonal_deg: int = ...,
trend_deg: int = ...,
low_pass_deg: int = ...,
robust: bool = ...,
seasonal_jump: int = ...,
trend_jump: int = ...,
low_pass_jump: int = ...,
inner_iter: int | None = ...,
outer_iter: int | None = ...,
) -> dict[str, Any]:
STL seasonal-trend decomposition using LOESS (Cleveland et al. 1990).
Mirrors statsmodels.tsa.seasonal.STL parameter semantics and defaults
exactly (matched elementwise at 1e-8; observed ~1e-12); requires
n >= 2*period. Returns `seasonal`, `trend`, `resid` (y = seasonal +
trend + resid), `weights` (bisquare robustness weights; all 1 unless
the outer loop runs), `period`, and `config` (the resolved windows,
degrees, jumps, and inner/outer iteration counts).
Further arguments, with defaults: `low_pass` (None), `seasonal_deg` (1),
`trend_deg` (1), `low_pass_deg` (1), `robust` (False), `seasonal_jump`
(1), `trend_jump` (1), `low_pass_jump` (1), `inner_iter` (None),
`outer_iter` (None).
mstl¶
def mstl(
y: _ArrayLike,
periods: Sequence[int],
windows: Sequence[int] | None = ...,
iterate: int = ...,
trend: int | None = ...,
low_pass: int | None = ...,
seasonal_deg: int = ...,
trend_deg: int = ...,
low_pass_deg: int = ...,
robust: bool = ...,
seasonal_jump: int = ...,
trend_jump: int = ...,
low_pass_jump: int = ...,
inner_iter: int | None = ...,
outer_iter: int | None = ...,
) -> dict[str, Any]:
MSTL — STL iterated over multiple seasonal periods
(Bandara-Hyndman-Bergmeir 2021), e.g. periods=[24, 168] for hourly
data with daily and weekly cycles.
Matches statsmodels.tsa.seasonal.MSTL elementwise at 1e-8: periods
sorted ascending, any period >= n/2 dropped (reported in
`dropped_periods`), per-period seasonal windows from `windows` (None:
the 7 + 4*k rule -> 11, 15, 19, ...), `iterate` refinement rounds
(default 2; 1 for a single period), remaining STL keywords forwarded
to every pass. statsmodels' Box-Cox `lmbda` option is not implemented
(pre-transform y instead); duplicate periods and iterate=0 are
refused. Returns `seasonal` (dict of per-period arrays keyed
"seasonal_<period>"), `trend`, `resid`, `weights` (from the final
pass), the resolved `periods`/`windows`, `iterate`,
`dropped_periods`, and per-period `seasonal_strength` (None for a
constant series).
Further arguments, with defaults: `low_pass` (None), `seasonal_deg` (1),
`trend_deg` (1), `low_pass_deg` (1), `robust` (False), `seasonal_jump`
(1), `trend_jump` (1), `low_pass_jump` (1), `inner_iter` (None),
`outer_iter` (None).
seasonal_strength¶
Wang-Smith-Hyndman seasonal/trend strength from a default STL fit.
strength = max(0, 1 - var(resid)/var(component + resid)), sample
variances; near 1 means the component dominates. Returns
`seasonal_strength`, `trend_strength`, `period`.
forecasting / evaluation¶
dm_test¶
Diebold-Mariano test with the Harvey-Leybourne-Newbold correction.
Returned keys: `dm_stat`, `hln_stat`, `mean_loss_diff`, `p_value`.
Further arguments, with defaults: `h` (1), `loss` ("squared").
accuracy¶
def accuracy(
actual: _ArrayLike,
forecast: _ArrayLike,
insample: _ArrayLike | None = ...,
period: int = ...,
) -> dict[str, float]:
Forecast accuracy measures (ME/RMSE/MAE/MAPE/sMAPE/MASE/RMSSE).
Returned keys: `mae`, `mape`, `mase`, `me`, `rmse`, `rmsse`, `smape`.
Further arguments, with defaults: `insample` (None), `period` (1).
theta_forecast¶
The Theta method (Assimakopoulos-Nikolopoulos 2000).
Matches statsmodels `ThetaModel(deseasonalize=True, use_test=False)`;
statsmodels' default additionally pre-tests seasonality and skips
deseasonalization when the test fails, so the two defaults diverge on
weakly-seasonal data declared with `period > 1`.
local projections¶
lp¶
def lp(
y: _ArrayLike,
shock: _ArrayLike,
horizons: int = ...,
n_lag_controls: int = ...,
se: str | None = ...,
maxlags: int | None = ...,
cumulative: bool | str | None = ...,
band: str | None = ...,
band_alpha: float = ...,
band_seed: int = ...,
band_n_sim: int = ...,
) -> dict[str, Any]:
Local projection IRFs; se is None (auto), "lag_augmented" or "hac".
`se=None` (the default) resolves to "lag_augmented" — except under
`cumulative="both"`, where it resolves to "hac": the cumulated impulse
`sum_(j=0..h) shock_(t+j)` shares FUTURE shocks across base times up to h
apart, which past-lag augmentation cannot project out, so lag-augmented
HC1 standard errors are inconsistent there (audit: 0.507 coverage at a
nominal 95%, h=12, flat in T) and `se="lag_augmented"` with
`cumulative="both"` raises. The method actually used is returned as
`se_method`.
`cumulative`: False/"none" (level), True/"outcome" (cumulated outcome on
the contemporaneous impulse — a cumulative IRF, NOT a multiplier), or
"both" (cumulated outcome on cumulated impulse). For an identified
multiplier use `lp_multiplier`.
**Bands.** `band=None` (default) returns the point path and its standard
errors only, exactly as before. `"pointwise"`, `"sup-t"`, `"sidak"` or
`"bonferroni"` add `lower`/`upper`, `critical_value`,
`pointwise_critical_value`, `band_scope`, `n_cells` (K), `n_cells_used`
and `cov_se_max_rel_diff` (largest relative gap between the band
covariance's sqrt(diag) and the reported `se`; ~machine epsilon on the
lag-augmented sup-t path, up to a few percent on the HAC sup-t path,
None where no covariance is built).
The family is **the horizons of this one response**, `K = horizons + 1`
(`band_scope` reports `"horizon"`). A pointwise band covers one horizon at a
time; the other three cover every horizon at once at `1 - band_alpha`.
LP is the clean case for this. Audit, nominal 90% over 13 horizons, 400
replications: pointwise contained the whole path in 36.5% of samples at
T=240 and sup-t in 81.8%; at T=720, where the per-horizon marginals sit on
nominal, sup-t lands on nominal too (89.5%) while pointwise reached only
42.7%. Tripling the sample moved pointwise 36.5% → 42.7% — it is not
converging, because the problem is multiplicity, not consistency.
`"sup-t"` builds the cross-horizon covariance and simulates `band_n_sim`
Gaussian draws from `band_seed`, so the band is a **pure function** of that
seed; the closed forms use neither. Measured at K=13, alpha=0.10: pointwise
1.6449, sup-t 2.20–2.65 depending on persistence, Šidák 2.6490, Bonferroni
2.6653. Method: Montiel Olea and Plagborg-Møller.
Returned keys: `horizons`, `irf`, `se`, `se_method`.
Further arguments, with defaults: `n_lag_controls` (4), `maxlags`
(None).
lp_iv¶
def lp_iv(
y: _ArrayLike,
impulse: _ArrayLike,
instrument: _ArrayLike,
horizons: int = ...,
n_lag_controls: int = ...,
cumulative: bool | str | None = ...,
band: str | None = ...,
band_alpha: float = ...,
) -> dict[str, Any]:
LP-IV: instrumented local projections with a first-stage F diagnostic.
`cumulative` takes False/"none", True/"outcome" or "both". True/"outcome"
cumulates only the OUTCOME, giving cumulated response per unit of
*contemporaneous* impulse — that grows without bound in the horizon and is
not a multiplier. Use `lp_multiplier` for the Ramey-Zubairy integral
multiplier.
**Bands.** `band=None` (default) returns no band. `"pointwise"`, `"sidak"`
and `"bonferroni"` add `lower`/`upper` over the horizons of this response
(`K = horizons + 1`, `band_scope="horizon"`) with `critical_value`,
`pointwise_critical_value`, `n_cells`, `n_cells_used` and
`cov_se_max_rel_diff` (always None here: no covariance is built).
`band="sup-t"` is **refused** here with an error saying why: sup-t needs the
covariance ACROSS horizons and tsecon estimates none for LP-IV, so `lp_iv`,
`lp_multiplier` and `lp_state` get the **closed-form** simultaneous routes
only. Šidák and Bonferroni need nothing but K, are valid under arbitrary
dependence, and are simply wider than a sup-t band would be — never describe
a band from this function as sup-t. For sup-t use `lp` or `smooth_lp`.
Returned keys: `first_stage_f`, `horizons`, `irf`, `se`.
Further arguments, with defaults: `n_lag_controls` (4), `band_alpha`
(0.1).
lp_multiplier¶
def lp_multiplier(
y: _ArrayLike,
impulse: _ArrayLike,
instrument: _ArrayLike,
horizons: int = ...,
n_lag_controls: int = ...,
maxlags: int | None = ...,
band: str | None = ...,
band_alpha: float = ...,
) -> dict[str, Any]:
Ramey-Zubairy (2018) integral multiplier by one-step LP-IV.
Regresses the cumulated outcome on the cumulated impulse, instrumented by
the contemporaneous instrument, controlling for lags of both series. Both
sides accumulate over the same window, so the coefficient is a multiplier
rather than a cumulative impulse response. `se` is the kernel-HAC standard
error of that single 2SLS coefficient — inference on the multiplier
itself, not a delta-method ratio and not a leg's SE relabelled.
**Bands.** `band=None` (default) returns no band. `"pointwise"`, `"sidak"`
and `"bonferroni"` add `lower`/`upper` around `multiplier` over the horizons
of this path (`K = horizons + 1`, `band_scope="horizon"`) with
`critical_value`, `pointwise_critical_value`, `n_cells`, `n_cells_used` and
`cov_se_max_rel_diff` (always None here: no covariance is built).
`band="sup-t"` is **refused**: no cross-horizon covariance is estimated for
the multiplier path, so this function (like `lp_iv` and `lp_state`) gets the
closed-form routes only. Do not call such a band sup-t.
Returned keys: `cumulative_impulse`, `cumulative_outcome`,
`first_stage_f`, `horizons`, `multiplier`, `nobs_per_h`, `se`.
Further arguments, with defaults: `n_lag_controls` (4), `maxlags`
(None), `band_alpha` (0.1).
penalized regression¶
ridge¶
Ridge regression (closed form); scikit-learn Ridge objective.
elastic_net¶
def elastic_net(
x: _ArrayLike,
y: _ArrayLike,
alpha: float,
l1_ratio: float = ...,
tol: float = ...,
max_iter: int = ...,
) -> dict[str, Any]:
Elastic-net via coordinate descent; scikit-learn objective.
Keys: `coef`, `n_iter`, `max_change` (largest absolute coefficient update
in the final sweep, in coefficient units) and `max_rel_change` (that
update scaled as max_j |Δb_j|·‖x_j‖/‖y‖ — the scale-free quantity the
stopping rule compares with `tol`).
Further arguments, with defaults: `l1_ratio` (0.5), `max_iter` (100000).
lasso¶
def lasso(
x: _ArrayLike,
y: _ArrayLike,
alpha: float,
tol: float = ...,
max_iter: int = ...,
) -> dict[str, Any]:
Lasso (elastic net with l1_ratio = 1.0).
Keys: `coef`, `n_iter`, `max_change` and `max_rel_change` (as in
`elastic_net`: the scale-free update the stopping rule compares with
`tol`).
Further arguments, with defaults: `max_iter` (100000).
structural identification¶
sign_restricted_svar¶
def sign_restricted_svar(
data: _ArrayLike,
restrictions: Sequence[tuple[int, int, int, str]],
lags: int = ...,
horizon: int = ...,
n_draws: int = ...,
max_tries: int = ...,
seed: int = ...,
lambda1: float = ...,
) -> dict[str, Any]:
Sign-restricted Bayesian SVAR: identified-set bands + acceptance diagnostics.
`restrictions` are (variable, shock, horizon, sign) tuples with sign in
{"+", "-"}. Returns per-(horizon, variable, shock) `quantiles` at
`probs=[0.05,0.16,0.50,0.84,0.95]`, the identified-set envelope
(`set_min`/`set_max`), and `diagnostics`.
Further arguments, with defaults: `lags` (2), `n_draws` (500),
`max_tries` (400), `seed` (0), `lambda1` (0.2).
Returned keys: `diagnostics`, `probs`, `quantiles`, `set_max`,
`set_min`.
zero_sign_svar¶
def zero_sign_svar(
data: _ArrayLike,
sign_restrictions: Sequence[tuple[int, int, int, str]],
zero_restrictions: Sequence[tuple[int, int, int]],
lags: int = ...,
horizon: int = ...,
n_draws: int = ...,
max_tries: int = ...,
seed: int = ...,
lambda1: float = ...,
weighted: bool = ...,
) -> dict[str, Any]:
Zero + sign restricted Bayesian SVAR: exact zeros by construction + sign rejection.
`sign_restrictions` are (variable, shock, horizon, sign) tuples with sign in
{"+", "-"} (may be empty); `zero_restrictions` are (variable, shock, horizon)
tuples imposing an exact zero on `Theta_h[variable, shock]` (horizon 0 =
impact). At least one list must be non-empty. Returns per-(horizon, variable,
shock) `quantiles` at `probs=[0.05,0.16,0.50,0.84,0.95]` (ARW-2018 importance-
weighted when `weighted=True`), the weight-invariant identified-set envelope
(`set_min`/`set_max`), per-accepted-draw `weights` (normalized to sum to 1) and
their effective sample size `ess`, the acceptance `diagnostics`, and
`arw_weighted` (whether the ARW weights were applied to `quantiles`). With
strict-upper-triangle impact zeros and no sign restrictions the rotation at
every draw is pinned to Q=I, so each posterior draw's structural IRF equals
that draw's recursive Cholesky IRF (a per-draw identity checked to ~1e-10 in
the crate golden); the posterior of the bands therefore coincides with the
recursive-Cholesky posterior, and the reported `set_min`/`set_max` span
reflects posterior (not identified-set) uncertainty since the rotation is
fixed. The ARW weight is exactly 1 for impact-only zero patterns.
Further arguments, with defaults: `lags` (2), `n_draws` (500),
`max_tries` (400), `seed` (0), `lambda1` (0.2).
Returned keys: `arw_weighted`, `diagnostics`, `ess`, `probs`,
`quantiles`, `set_max`, `set_min`, `weights`.
structural_fevd¶
def structural_fevd(
data: _ArrayLike,
lags: int = ...,
horizon: int = ...,
trend: str = ...,
impact: _ArrayLike | None = ...,
sigma: str = ...,
) -> dict[str, Any]:
Structural FEVD for an arbitrary structural impact matrix A0 (the gap var_fevd, recursive-Cholesky only, leaves).
`impact` is an optional (n, n) structural impact A0 (columns = one-SD
structural shocks, A0 A0' = Sigma; from any identification scheme). If None,
A0 is the lower Cholesky of the innovation covariance and the result equals
`var_fevd` exactly. `sigma` ("dfadj"|"mle") sets the default Cholesky's df
scaling; the FEVD shares are invariant to it (it only rescales the reported
`impact`). Returns `fevd` [horizon+1][variable][shock] (each row sums to 1)
and `impact` [n][n] (the A0 used).
Further arguments, with defaults: `lags` (2), `trend` ("c").
historical_decomposition¶
def historical_decomposition(
data: _ArrayLike,
restrictions: Sequence[tuple[int, int, int, str]] | None = ...,
lags: int = ...,
horizon: int | None = ...,
identification: str = ...,
n_draws: int = ...,
max_tries: int = ...,
seed: int = ...,
lambda1: float = ...,
narrative_restrictions: list[dict] | None = ...,
n_weight_draws: int = ...,
) -> dict[str, Any]:
Historical decomposition: per-(time, variable, shock) structural-shock contributions.
Splits each variable into a deterministic/initial-condition `baseline` plus the
cumulated contribution `hd[time][variable][shock]` of each structural shock,
obeying the exact adding-up identity y = baseline + sum_j hd (validated to ~1e-10
against a NumPy reference). `times` are 0-based effective-sample indices
(= data_row - lags).
identification="cholesky" (default): a point decomposition at the OLS VAR with
Q=I; returns `times`, `baseline` [T_eff][n], `hd` [T_eff][n][n] indexed
[time][variable][shock], and the structural `shocks` [T_eff][n].
identification="sign": the importance-weighted SET decomposition over sign- (and
optionally narrative-) restricted rotations; returns `times`, `baseline`
(posterior-mean), `probs`, `hd_quantiles` [T_eff][n][n][len(probs)] (weighted
type-7), the weight-free identified-set envelope `hd_set_min`/`hd_set_max`,
per-draw `weights`, and `diagnostics`.
`narrative_restrictions` (sign mode) is a list of dicts with 0-based effective
indices:
{"type":"shock_sign","shock":int,"period":int,"sign":"+"|"-"}
{"type":"contribution","variable":int,"shock":int,"start":int,"end":int,
"rule":"most"|"least","strong":bool}
{"type":"contribution_sign","variable":int,"shock":int,"start":int,"end":int,
"sign":"+"|"-"}
Further arguments, with defaults: `horizon` (None), `n_draws` (500),
`max_tries` (400), `seed` (0), `lambda1` (0.2), `n_weight_draws` (200).
`restrictions` (None = no restrictions; used under identification="sign")
is the `sign_restricted_svar` list of `(variable, shock, horizon, sign)`
tuples with `sign` in {"+", "-"}.
narrative_svar¶
def narrative_svar(
data: _ArrayLike,
sign_restrictions: Sequence[tuple[int, int, int, str]] | None = ...,
narrative_restrictions: list[dict] | None = ...,
lags: int = ...,
horizon: int = ...,
n_draws: int = ...,
max_tries: int = ...,
seed: int = ...,
lambda1: float = ...,
n_weight_draws: int = ...,
) -> dict[str, Any]:
Narrative sign-restricted Bayesian SVAR (Antolín-Díaz & Rubio-Ramírez 2018).
Augments traditional sign restrictions with restrictions on named historical
episodes — shock signs and "most/least important contributor" statements (see
`historical_decomposition` for the `narrative_restrictions` dict schema) —
imposed by importance-reweighting the accepted rotations with weight = 1/P̂(N|S).
Returns per-(horizon, variable, shock) `quantiles` (weighted type-7) at
`probs=[0.05,0.16,0.50,0.84,0.95]`, the weight-free identified-set envelope
`set_min`/`set_max`, per-draw `weights` (mean 1), and `diagnostics` (with `ess`,
`narrative_acceptance_rate`, `min_ptilde`). With no narrative restrictions every
weight is 1 and it reproduces `sign_restricted_svar` bit-for-bit.
Further arguments, with defaults: `lags` (2), `n_draws` (500),
`max_tries` (400), `seed` (0), `lambda1` (0.2), `n_weight_draws` (200).
`sign_restrictions` (None = no restrictions) is the `sign_restricted_svar`
list of `(variable, shock, horizon, sign)` tuples with `sign` in
{"+", "-"}.
Returned keys: `diagnostics`, `probs`, `quantiles`, `set_max`,
`set_min`, `weights`.
fry_pagan_svar¶
def fry_pagan_svar(
data: _ArrayLike,
restrictions: Sequence[tuple[int, int, int, str]],
lags: int = ...,
horizon: int = ...,
n_draws: int = ...,
max_tries: int = ...,
seed: int = ...,
lambda1: float = ...,
target: str = ...,
) -> dict[str, Any]:
Fry-Pagan (2011) median-target SVAR: the single coherent draw closest to the median band.
Sign restrictions set-identify a *set* of structural models; the pointwise
median band mixes responses from mutually inconsistent draws and is not
itself a model. This returns instead the single accepted, sign-normalized
draw whose structural IRFs jointly minimize the Fry-Pagan criterion -- the
sum, over the target cells, of squared deviations from the pointwise median,
each standardized by that cell's across-draw dispersion. `restrictions` are
(variable, shock, horizon, sign) tuples with sign in {"+", "-"}; `target` is
"restricted" (response cells of the sign-restricted shocks; default) or
"all". Returns the coherent `median_target_irf` [horizon+1][n][n], the
incoherent pointwise `median_irf` (for comparison), the selected `mt_index`
(0-based into the accepted set), its `mt_statistic`, `n_accepted`, and the
acceptance `diagnostics`. Reproducible at a fixed `seed` (substream
contract). The selected draw is a descriptive summary -- one interior point
of the identified set, dependent on the informative Haar prior -- not a
prior-free point estimate.
Further arguments, with defaults: `lags` (2), `n_draws` (500),
`max_tries` (400), `lambda1` (0.2).
robust_svar_bounds¶
def robust_svar_bounds(
data: _ArrayLike,
restrictions: Sequence[tuple[int, int, int, str]],
lags: int = ...,
horizon: int = ...,
n_draws: int = ...,
seed: int = ...,
lambda1: float = ...,
alpha: float = ...,
) -> dict[str, Any]:
Giacomini-Kitagawa prior-robust identified-set bounds for a sign-restricted SVAR.
`restrictions` are (variable, shock, horizon, sign) tuples with sign in
{"+", "-"}. For each restricted shock, the per-draw identified set of the
structural IRF is computed exactly over the admissible rotation set and
summarized over the reduced-form posterior, removing the informative-Haar-
prior artifact that pointwise `sign_restricted_svar` bands carry. Returns
per (horizon, variable, shock): `set_lower_mean`/`set_upper_mean` (posterior-
mean identified-set edges), `robust_ci_lower`/`robust_ci_upper` (the level-
`alpha` robust credible region), and `lower_quantiles`/`upper_quantiles` at
`probs=[0.05,0.16,0.50,0.84,0.95]`. Unrestricted shocks are NaN;
`restricted_shocks` lists the valid shock indices; `diagnostics` reports
`empty_set_rate` (the share of draws whose restrictions were mutually
infeasible). Exact for a single restricted shock (Gafarov-Meier-Montiel-Olea
2018 closed form); with multiple jointly-restricted shocks each bound is that
shock's marginal identified set — a conservative outer approximation of the
joint set, since the cross-shock orthogonality coupling is not imposed.
Further arguments, with defaults: `lags` (2), `n_draws` (500), `seed`
(0), `lambda1` (0.2).
Returned keys: `alpha`, `diagnostics`, `lower_quantiles`, `probs`,
`restricted_shocks`, `robust_ci_lower`, `robust_ci_upper`,
`set_lower_mean`, `set_upper_mean`, `upper_quantiles`.
long_run_svar¶
def long_run_svar(
data: _ArrayLike,
lags: int = ...,
horizon: int = ...,
trend: str = ...,
restrictions: Sequence[tuple[int, int]] | None = ...,
normalize: str = ...,
) -> dict[str, Any]:
Blanchard-Quah long-run SVAR: closed-form structural IRFs under frequency-zero restrictions.
`restrictions` is a list of (variable, shock) long-run zero pairs (None =>
classic recursive BQ); `normalize` is "long_run" (positive LR diagonal;
default) or "impact" (positive B diagonal). Returns `impact` (B),
`long_run` (LR = C(1) B), `long_run_multiplier` (C(1)), `irf`
[horizon+1][i][j], `cumulative_irf`, and `fevd`. Point estimate, no RNG.
Further arguments, with defaults: `lags` (2), `trend` ("c").
max_share_svar¶
def max_share_svar(
data: _ArrayLike,
lags: int = ...,
target: int = ...,
h0: int = ...,
h1: int = ...,
horizon: int = ...,
trend: str = ...,
exclude_impact: bool = ...,
weighting: str = ...,
sign: str = ...,
) -> dict[str, Any]:
Max-share / maximum-FEV structural shock (Uhlig 2004; Francis et al 2014; Barsky-Sims 2011 news).
Identifies the single UNIT-VARIANCE structural shock maximizing the `target`
variable's forecast-error variance accumulated over the window `[h0, h1]`.
`weighting="window"` selects the Uhlig/Francis objective (incremental
windowed FEV; `share_window` is an exact accumulated-FEV fraction),
`"cumulative"` the Barsky-Sims objective (window-mean cumulative FEV share).
`exclude_impact=True` imposes zero impact on the target (Barsky-Sims news
shock). `sign` pins the identified sign ("cumsum"|"impact"|"none").
Returns `irf` [horizon+1][k], `impact` [k], `q` [k], `share_window` (float),
`fev_share` [horizon+1], and `eigenvalues` (ascending; length k, or k-1 when
`exclude_impact`).
Further arguments, with defaults: `lags` (2), `trend` ("c").
proxy_svar_bands¶
def proxy_svar_bands(
data: _ArrayLike,
proxy: _ArrayLike,
lags: int = ...,
horizon: int = ...,
norm_var: int = ...,
unit: float = ...,
trend: str = ...,
alpha: float = ...,
n_boot: int = ...,
seed: int = ...,
bands: str = ...,
block_length: int | None = ...,
robust_f: bool = ...,
) -> dict[str, Any]:
Confidence bands for a proxy-SVAR impulse response.
bands="moving_block" (default) is the Jentsch-Lunsford moving-block
bootstrap: (u_t, m_t) resampled jointly, the VAR reconstructed and
re-estimated per draw, the unit-effect normalization re-imposed per draw.
bands="wild" reproduces Mertens-Ravn / Gertler-Karadi but is NOT
asymptotically valid here -- a common Rademacher draw leaves the
identifying moment bit-identical across draws, so it carries no bootstrap
variability. Check asymptotically_valid / validity_note.
Returns lower/upper (Hall, recommended) and lower_efron/upper_efron.
The h=0 entry for norm_var is degenerate at `unit` by construction.
Bands are pointwise, not joint. Failed draws are counted by reason in
`failures`, never dropped; a nonzero n_failed means the instrument may be
too weak for a Wald band -- see proxy_ar_sets.
Returned keys: `alpha`, `asymptotically_valid`, `block_length`,
`failure_warning`, `failures`, `first_stage_f_draws`,
`gamma_norm_draws`, `lower`, `lower_efron`, `method`, `n_boot`,
`n_failed`, `n_proxy`, `n_used`, `point`, `point_first_stage_f`,
`point_gamma_norm`, `point_reliability`, `reliability_draws`,
`rho_draws`, `se`, `upper`, `upper_efron`, `validity_note`.
Further arguments, with defaults: `lags` (2), `horizon` (12), `trend`
("c"), `alpha` (0.1), `n_boot` (2000), `seed` (0), `block_length`
(None), `robust_f` (True).
proxy_ar_sets¶
def proxy_ar_sets(
data: _ArrayLike,
proxy: _ArrayLike,
lags: int = ...,
horizon: int = ...,
norm_var: int = ...,
unit: float = ...,
trend: str = ...,
alpha: float = ...,
variance: str = ...,
hac_lags: int | None = ...,
reduced_form_uncertainty: bool = ...,
rf_method: str = ...,
rf_draws: int | None = ...,
rf_seed: int | None = ...,
) -> dict[str, Any]:
Weak-instrument-robust (Anderson-Rubin) confidence SETS for a proxy SVAR.
Under weak identification no bounded set can be honest (Dufour 1997), so a
cell may be a bounded interval, the COMPLEMENT of an interval (kind
"exterior", two rays), a single ray ("ray_below"/"ray_above"), the whole
line, empty, or a point. That shape is the answer.
Do not read an "exterior" set as an interval -- `lower`/`upper` are the
set's own bounds (+/-inf there) and `excluded_lower`/`excluded_upper` are
the rejected middle. `excludes_zero` on an unbounded set does NOT establish
a sign: both signs can be members.
Reduced-form uncertainty is propagated by default. Omitting it is
catastrophic on an estimated VAR -- measured coverage 0.952 at h=0 falling
to 0.119 by h=8 against nominal 0.95, versus 0.952 to 0.913 with it. When
reduced_form_uncertainty=False the returned `level` is None, because a set
conditional on the reduced form has no honest 1-alpha label.
rf_method="second_order" (with rf_draws/rf_seed) replaces the first-order
delta propagation with seeded exact simulation of the coefficient
uncertainty through the nonlinear MA map -- the measured long-horizon
repair (h=12 coverage 0.889 -> 0.964 on the card's VAR(2) at T=300, 0.830
-> 0.932 on a routine VAR(1) at T=250; median width ~1.15x at h=8, ~1.45x
at h=12; weak-instrument boundedness bit-identical). Default "delta" is
unchanged. `rf_seed=None` (the default) means seed 0 — reproducible, not
fresh entropy — and `rf_draws=None` means 256 draws.
rf_method="second_order_bc" centres the same seeded simulation at Pope
(1990) bias-corrected coefficients (Kilian stationarity shrinkage) --
measured at-or-above nominal at EVERY horizon on both DGPs (h=12: 0.982
card VAR(2) / 0.966 routine VAR(1)) at a further width price (median
~1.8x the delta width at h=12). A conservative floor, not a calibration:
it overshoots where "second_order" already reaches nominal. Boundedness
is again bit-identical.
Proxy missingness follows the family convention: NaN rows are treated as
dates where the instrument is unavailable and dropped from the moments
(`n_proxy` reports the kept count); an infinite proxy value is refused as
corruption (0.7.0 -- previously dropped as if missing).
Returned keys: `ar_bound_stat`, `ar_bounded_all`, `cells`,
`critical_value`, `first_stage_f`, `impact`, `level`, `n_proxy`,
`reduced_form_uncertainty`.
Further arguments, with defaults: `lags` (2), `norm_var` (0), `unit`
(1.0), `trend` ("c"), `variance` ("hc0"), `hac_lags` (None).
proxy_svar¶
def proxy_svar(
data: _ArrayLike,
proxy: _ArrayLike,
lags: int = ...,
horizon: int = ...,
norm_var: int = ...,
unit: float = ...,
trend: str = ...,
robust_f: bool = ...,
) -> dict[str, Any]:
Proxy SVAR (external-instrument SVAR-IV): one shock from one instrument.
The residual-instrument covariance identifies the target shock's impact
column up to scale; the unit-effect normalization sets its impact on
`norm_var` to `unit` (sign pinned). `proxy` aligns to `data` rows (NaN
outside the instrument window is dropped, `n_proxy` counts the kept
rows; an infinite proxy value is refused as corruption — 0.7.0).
Returns `irf` (horizon+1, n),
`impact`, `relative_impact`, `cov_um`, `first_stage_f` (HC1-robust when
`robust_f`), `reliability` = Corr(m, u_norm)^2, `n_proxy`, the estimated
`shock` (length T - lags, the residual sample), and `first_stage`: the proxy_first_stage diagnostics
dict (the MOP effective F with its tau-based critical values --
mop_cv_tau10 = 23.11 is the conventional bar, not the folklore 10).
Point estimate only; see proxy_svar_bands for moving-block bands
(strong instrument) and proxy_ar_sets for weak-IV-robust sets (use when
first_stage["weak_mop_tau10"] is True).
Further arguments, with defaults: `trend` ("c").
proxy_first_stage¶
def proxy_first_stage(
data: _ArrayLike,
proxy: _ArrayLike,
lags: int = ...,
norm_var: int = ...,
trend: str = ...,
variance: str = ...,
hac_lags: int | None = ...,
) -> dict[str, Any]:
First-stage strength diagnostics: the Montiel Olea-Pflueger effective F.
With one instrument the MOP effective F equals the robust F (the squared
robust t of the first-stage slope; Windmeijer 2025), reported under
variance="hc1" (default), "hac" (Bartlett/Newey-West, hac_lags defaulting
to the Newey-West rule -- for serially correlated proxies), or
"classical" (for comparison with published homoskedastic tables).
Returns `beta`, `se`, `effective_f`, `f_classical`, `f_hc1`,
`reliability`, `n_proxy`, `hac_lags`, the MOP critical values at 5% test
level (`mop_cv_tau5/10/20/30` = 37.42 / 23.11 / 15.06 / 12.05 -- the
null "worst-case relative bias > tau"), `tau_bound` (the smallest tau the
observed effective F rejects; +inf when even zero relevance cannot be
rejected), and the verdicts `weak_mop_tau10` (the honest bar) and
`weak_folklore` (F < 10, kept only because the literature reports it).
When weak_mop_tau10 is True do not trust Wald-type bands
(proxy_svar_bands); use proxy_ar_sets.
Proxy missingness follows the family convention: NaN rows are treated as
dates where the instrument is unavailable and dropped from the
first-stage sample (`n_proxy` reports the kept count); an infinite proxy
value is refused as corruption (0.7.0 -- previously dropped as if
missing).
Returned keys: `beta`, `effective_f`, `f_classical`, `f_hc1`,
`hac_lags`, `mop_cv_tau10`, `mop_cv_tau20`, `mop_cv_tau30`,
`mop_cv_tau5`, `n_proxy`, `reliability`, `se`, `tau_bound`,
`weak_folklore`, `weak_mop_tau10`.
Further arguments, with defaults: `lags` (2), `norm_var` (0), `trend`
("c").
nongaussian_svar¶
def nongaussian_svar(
data: _ArrayLike,
lags: int = ...,
horizon: int = ...,
trend: str = ...,
contrast: str = ...,
max_iter: int = ...,
tol: float = ...,
order_by: str = ...,
) -> dict[str, Any]:
Non-Gaussian / independent-component SVAR identification (Lanne-Meitz-Saikkonen 2017; Gourieroux-Monfort-Renne 2017; FastICA).
Point-identifies the structural impact matrix B in u_t = B eps_t from the
reduced-form residuals ALONE -- no sign, zero, long-run, or proxy
restriction -- by exploiting the statistical INDEPENDENCE and NON-GAUSSIANITY
of the structural shocks (at most one Gaussian). Whitens by Sigma_u^{-1/2},
finds the orthogonal rotation maximizing non-Gaussianity via a deterministic
symmetric FastICA fixed point (log-cosh contrast, identity init -- bit-
reproducible), then B = Sigma_u^{1/2} Q. Columns are ordered by `order_by`
("kurtosis" = descending |excess kurtosis|, or "colnorm") and signed max-abs-
positive; both are CONVENTIONS, not economics. This is STATISTICAL
identification: it FAILS if the shocks are Gaussian, and a `shock_kurtosis`
near zero flags a weakly identified (near-Gaussian) column. Returns `impact`
(B, [var][shock]), `irf` ([horizon+1][var][shock], Theta_h = Psi_h B),
`rotation` (Q, [whitened][shock]), `shock_kurtosis` [k] (identified order),
`converged` (bool), `n_iter` (int), and `order` [k] (raw FastICA index per
identified position).
Further arguments, with defaults: `lags` (2), `trend` ("c"), `max_iter`
(200), `tol` (1e-08).
hetero_svar¶
def hetero_svar(
data: _ArrayLike,
regime_labels: npt.NDArray[np.integer] | Sequence[int],
lags: int = ...,
horizon: int = ...,
trend: str = ...,
base_regime: int | None = ...,
sign_normalization: str = ...,
) -> dict[str, Any]:
SVAR identification through heteroskedasticity (Rigobon 2003; Lanne-Lutkepohl 2008), two known variance regimes.
`data` is (T, n); `regime_labels` is an array-like of length T with EXACTLY
two distinct integer values (labels align to observations; the first `lags`
are dropped to match residuals). `base_regime` is the label normalized to
Lambda=I (default: the smaller label); the other regime's shock-variance
ratios are reported. `sign_normalization`: "max" (largest-|entry| per B
column made positive; default) or "diag" (B[j,j] >= 0).
Returns a dict with `B` (n x n impact matrix = Theta_0, columns in
ascending variance-ratio order), `variance_ratios` (the n generalized
eigenvalues, ascending), `structural_irf` ([h][i][j] = Theta_h = Psi_h B),
`min_ratio_gap` and `ratio_dist_from_unity` (identification margins),
`identified` (bool heuristic), `covariance_equality` (Bartlett-corrected
Box's M: statistic/dof/pvalue/distinct_regimes), `sigma_regime1`,
`sigma_regime2`, `regime1_label`, `regime2_label`, `regime_sizes`,
`n_vars`, `horizon`, `lags`, `sign_convention`.
Point-identified IF AND ONLY IF the variance ratios are pairwise distinct
(min_ratio_gap > 0); the shocks come out ordered by variance ratio and
carry no economic labels. Standard errors on B/Theta_h are not provided in
this closed-form build. The >2-regime and Markov-switching/GARCH variants
are deferred.
Further arguments, with defaults: `trend` ("c").
panel¶
panel_fe¶
def panel_fe(
outcome: _ArrayLike,
regressors: _ArrayLike,
se_type: str = ...,
bandwidth: float | None = ...,
mask: _ArrayLike | None = ...,
) -> dict[str, Any]:
Fixed-effects panel OLS; outcome is N x T, regressors is k x N x T.
`se_type`: "nonrobust", "cluster" (by entity), or "driscoll_kraay".
`bandwidth` is the Driscoll-Kraay lag truncation and acts ONLY under
`se_type="driscoll_kraay"` (4.0 when omitted there); passing it
explicitly with any other `se_type` raises instead of being silently
absorbed — those estimators use no kernel. `mask` (default None = a
balanced panel) is an N x T array of 0/1 (False/True) flags, 1 where
the entity is observed in that period, for an UNBALANCED panel: cells
outside the mask are ignored and may hold NaN; without a mask a NaN
anywhere is refused. Validated against linearmodels PanelOLS on the
Arellano-Bond EmplUK panel (fixtures/panel_unbalanced.json); a mask
that is 1 everywhere is bit-identical to no mask.
Returned keys: `bse`, `params`, `se_type`, `tvalues`.
panel_lp¶
def panel_lp(
outcome: _ArrayLike,
shock: _ArrayLike,
horizon: int = ...,
n_lag_controls: int = ...,
se_type: str = ...,
bandwidth: float | None = ...,
cumulative: bool = ...,
jackknife: bool = ...,
bias_correction: str = ...,
band: str | None = ...,
band_alpha: float = ...,
mask: _ArrayLike | None = ...,
) -> dict[str, Any]:
Panel local projection of a common shock with fixed effects.
`mask` (default None = a balanced panel) is an N x T array of 0/1
flags, 1 where the entity is observed in that period, for an UNBALANCED
panel: the horizon-h regression keeps the rows whose target (or
cumulated window) and lagged-outcome controls are observed, so `nobs`
shrinks with the gaps as well as the horizon (validated per horizon
against linearmodels PanelOLS, fixtures/panel_unbalanced.json). The
half-panel jackknives (`jackknife=True`, `bias_correction="dj"`/"spj")
raise on an unbalanced panel.
`bandwidth` is the Driscoll-Kraay lag truncation and acts ONLY under
`se_type="driscoll_kraay"` (the default se_type; 4.0 when omitted);
passing it explicitly with `se_type="cluster"`/`"nonrobust"` raises
instead of being silently absorbed.
`outcome` is N x T; `shock` is length T. Fixed effects + lagged outcomes
+ short T carry Nickell bias (horizon-amplified); two half-panel
corrections are offered. `jackknife=True` (equivalently
`bias_correction="dj"`) is the Dhaene-Jochmans half-panel jackknife:
corrected point estimates, full-sample plug-in standard errors
(measured cost: the estimator's variance inflates at short T while se
is unchanged — coverage 0.88 -> 0.80 at T=60, equivalence by T ~ 240).
`bias_correction="spj"` is the Mei-Sheng-Shi (2026, J. Int. Economics)
split-panel jackknife for panel LPs: leads/lags stay full-panel, the
regression rows split at the median usable period, and the standard
errors are recomputed for the corrected estimator (adjusted-score
cluster or Driscoll-Kraay sandwich, matching their pLP reference
implementation; `se_type="nonrobust"` is refused under "spj").
Combining `jackknife=True` with `bias_correction="spj"` raises.
Returns a dict with `irf`, `se`, `nobs` (each length horizon+1) and the
stamped `se_type`, `cumulative`, `jackknife`, `bias_correction`.
**Bands.** `band=None` (default) returns no band. `"pointwise"`, `"sidak"`
and `"bonferroni"` add `lower`/`upper` over the horizons of this response
(`K = horizon + 1`, `band_scope="horizon"`) at level `band_alpha`, with
`critical_value`, `pointwise_critical_value`, `n_cells`, `n_cells_used`
and `cov_se_max_rel_diff` (always None here: no covariance is built;
`band_n_sim`/`band_seed` come back 0 — no simulation ran). A pointwise
band covers one horizon at a time; the closed-form simultaneous routes
cover every horizon at once at `1 - band_alpha` (Montiel Olea and
Plagborg-Møller's simultaneous-bands framework; joint coverage measured
in `test_simultaneous_bands.py` — see the panel model card).
`band="sup-t"` is **refused** with an error saying why: sup-t needs the
covariance ACROSS horizons and tsecon estimates none for the panel LP
(a cross-horizon panel covariance is a documented follow-up), so
`panel_lp` gets the closed-form routes only, like `lp_iv`,
`lp_multiplier` and `lp_state`. Never describe such a band as sup-t.
Further arguments, with defaults: `n_lag_controls` (2).
lp_did¶
def lp_did(
outcome: _ArrayLike,
treatment: _ArrayLike,
pre_window: int = ...,
post_window: int = ...,
absorbing: bool = ...,
nonabsorbing_lag: int = ...,
reweight: bool = ...,
pooled: bool = ...,
never_treated_only: bool = ...,
mask: _ArrayLike | None = ...,
) -> dict[str, Any]:
LP-DiD event-study difference-in-differences (Dube-Girardi-Jordà-Taylor).
`mask` (default None) is accepted for symmetry with the other panel
callables, but LP-DiD needs a BALANCED panel: an unbalanced mask raises
with the reason (contiguous outcome paths; the reference fixest run was
made on balanced panels) — trim to a common window first.
`outcome` and `treatment` are N x T (treatment binary 0/1). Per horizon,
regresses `y[i, t+h] - y[i, t-1]` on the treatment switch with period
effects, using only clean controls (not-yet-treated; stabilized units
under `absorbing=False` with `nonabsorbing_lag`; never-treated when
`never_treated_only=True`) — avoiding the negative-weight comparisons of
TWFE event studies. `reweight=True` gives the equally-weighted ATT;
`pooled=True` adds pooled post/pre estimates. Entity-clustered SEs in
the authors' fixest/reghdfe convention.
Returns a dict with `horizons` (-pre_window..post_window; -1 is the
omitted baseline, stored as zeros), `coef`, `se`, `nobs`, `n_switchers`
(clean samples shrink with |h| — read them), pooled keys
(`pooled_post_att`, `pooled_post_se`, `pooled_post_nobs`,
`pooled_post_n_switchers`, and `pooled_pre_*` when `pre_window >= 2`)
only when `pooled=True`, and the stamped `absorbing`,
`nonabsorbing_lag`, `reweight`, `pooled`, `never_treated_only`,
`se_type`.
forecast comparison¶
cw_test¶
def cw_test(
e_small: _ArrayLike,
e_large: _ArrayLike,
yhat_small: _ArrayLike,
yhat_large: _ArrayLike,
lrv_lags: int = ...,
) -> dict[str, float]:
Clark-West test for nested-model equal predictive accuracy.
Returned keys: `cw_stat`, `mean_adj_diff`, `p_value`.
Further arguments, with defaults: `lrv_lags` (0).
gw_test¶
Giacomini-White unconditional test of equal predictive ability.
Returned keys: `df`, `gw_stat`, `p_value`.
Further arguments, with defaults: `lrv_lags` (0).
var_backtest¶
def var_backtest(
returns_or_hits: _ArrayLike,
var_forecasts: _ArrayLike | None = ...,
alpha: float = ...,
dq_lags: int = ...,
input: str = ...,
) -> dict[str, Any]:
VaR backtest battery: Kupiec unconditional coverage, Christoffersen independence/conditional coverage, and the Engle-Manganelli DQ test.
Sign convention: returns and VaR forecasts on the same (return) scale,
`var_forecasts[t]` the alpha-quantile of the conditional return
distribution (negative for small alpha); a violation is return < VaR.
`alpha` is the VaR coverage level (0.05 for a 95% VaR), not a test
size. With `var_forecasts` the first argument is a return series;
without, a pre-computed 0/1 violation sequence (`input="hits"`
combines pre-computed hits WITH VaR forecasts so the DQ regression
keeps its VaR regressor). Returns the three statistics with p-values,
the violation counts/transition cells, and a teaching `verdict`.
Returned keys: `alpha`, `dq_df`, `dq_includes_var`, `dq_lags`,
`dq_stat`, `dq_var_dropped`, `expected_violations`, `hit_rate`, `lr_cc`,
`lr_ind`, `lr_uc`, `n`, `n00`, `n01`, `n10`, `n11`, `n_violations`,
`p_cc`, `p_dq`, `p_ind`, `p_uc`, `pi01`, `pi11`, `verdict`.
Further arguments, with defaults: `dq_lags` (4).
spectral analysis¶
periodogram¶
def periodogram(
x: _ArrayLike, fs: float = ..., window: str = ..., detrend: str = ...
) -> dict[str, _F64]:
Periodogram PSD (freqs, psd); matches scipy.signal.periodogram.
Default `detrend="constant"` (mean removal) is SciPy's own default, so
default call matches default call; `"none"` / `"linear"` as in SciPy.
Further arguments, with defaults: `fs` (1.0), `window` ("boxcar").
Returned keys: `freqs`, `psd`.
welch¶
def welch(
x: _ArrayLike,
nperseg: int = ...,
fs: float = ...,
noverlap: int | None = ...,
window: str = ...,
detrend: str = ...,
) -> dict[str, _F64]:
Welch averaged-periodogram PSD; matches scipy.signal.welch.
Default `detrend="constant"` (per-segment mean removal) is SciPy's own
default, so default call matches default call.
Returned keys: `freqs`, `psd`.
Further arguments, with defaults: `nperseg` (256), `fs` (1.0),
`noverlap` (None), `window` ("hann").
coherence¶
def coherence(
x: _ArrayLike,
y: _ArrayLike,
nperseg: int = ...,
fs: float = ...,
noverlap: int | None = ...,
window: str = ...,
detrend: str = ...,
) -> dict[str, _F64]:
Magnitude-squared coherence in [0,1]; matches scipy.signal.coherence.
Default `detrend="constant"` (per-segment mean removal) is SciPy's own
default, so default call matches default call.
Returned keys: `coherence`, `freqs`.
Further arguments, with defaults: `nperseg` (256), `fs` (1.0),
`noverlap` (None), `window` ("hann").
cointegration¶
johansen¶
Johansen cointegration test (data is T x k); trace + max-eig + rank + evec.
Matches statsmodels ``coint_johansen(det_order=0)`` — the *unrestricted
constant* convention. Warning: ``vecm``'s default is ``deterministic="n"``
(no deterministic terms), a different case; fit the VECM this test ranks
with ``vecm(..., deterministic="co")``.
Returned keys: `eig`, `evec`, `max_eig_crit_90_95_99`, `max_eig_stat`,
`rank_max_eig_5pct`, `rank_trace_5pct`, `trace_crit_90_95_99`,
`trace_stat`.
Further arguments, with defaults: `k_ar_diff` (1).
engle_granger¶
def engle_granger(
data: _ArrayLike,
trend: str = ...,
autolag: str | None = ...,
maxlag: int | None = ...,
) -> dict[str, Any]:
Engle-Granger two-step cointegration test: stat + MacKinnon p-value/crit (statsmodels coint).
`data` is T x k, column 0 the regressand. Keys: `stat`, `pvalue`, `crit`,
`coint_coefs` (step-1 coefficients, deterministics first), `resid`
(length `nobs`), `used_lag`/`adf_nobs` (the residual ADF's lag and
sample), `n_vars`, `nobs`.
Further arguments, with defaults: `trend` ("c"), `maxlag` (None).
vecm¶
def vecm(
data: _ArrayLike,
k_ar_diff: int = ...,
coint_rank: int = ...,
deterministic: str = ...,
seasons: int = ...,
first_season: int | None = ...,
) -> dict[str, Any]:
VECM ML estimation: alpha, beta, det_coef_coint, gamma, det_coef, sigma_u, llf.
``deterministic`` names the statsmodels VECM case (all nine accepted):
``"n"`` (default) — no deterministic terms; ``"co"``/``"ci"`` — constant
outside/inside the cointegration relation; ``"lo"``/``"li"`` — linear
trend outside/inside; combinations ``"colo"``/``"coli"``/``"cilo"``/
``"cili"``. Restricted (inside) terms widen the cointegrating matrix —
their coefficients are returned as the rows of ``det_coef_coint``
(constant first, then trend; statsmodels ``VECMResults.det_coef_coint``);
unrestricted terms land in ``det_coef`` (statsmodels column order:
constant, seasons-1 centered seasonal dummies, trend). ``seasons``/
``first_season``: statsmodels-style centered seasonal dummies —
``first_season`` (0-based season of the first row, default 0 when None)
is taken modulo ``seasons`` (statsmodels-compatible), and passing it
explicitly with ``seasons=0`` raises (no cycle to phase). Warning:
``johansen`` assumes the unrestricted constant (det_order=0), NOT this
function's ``"n"`` default — pass ``deterministic="co"`` when the rank
came from ``johansen`` (det_order -1/0/1 ↔ ``"n"``/``"co"``/``"colo"``).
Further arguments, with defaults: `k_ar_diff` (1), `coint_rank` (1).
Returned keys: `alpha`, `beta`, `det_coef`, `det_coef_coint`, `gamma`,
`llf`, `sigma_u`.
threshold_vecm¶
def threshold_vecm(
data: _ArrayLike,
k_ar_diff: int = ...,
trim: float = ...,
n_grid_gamma: int = ...,
n_grid_beta: int | None = ...,
beta_span: float | None = ...,
beta: _ArrayLike | None = ...,
) -> dict[str, Any]:
Hansen-Seo (2002) two-regime threshold VECM (threshold
cointegration): the error-correction term w_{t-1} = beta' y_{t-1}
drives the regime split; estimation is the concentrated Gaussian MLE —
grid search over (beta, gamma) with per-cell two-regime OLS minimizing
ln det of the pooled residual covariance. trim is Hansen-Seo's pi0
(default 0.05, their suggestion). beta=None estimates the
cointegrating vector (BIVARIATE only, grid centered on the linear
Johansen estimate); for k > 2 pass beta= explicitly. With beta=
supplied the beta grid search never runs (beta_grid comes back
empty), so passing n_grid_beta/beta_span explicitly alongside
it raises (they default to 50/10.0 when None).
Keys: beta, threshold, params_low/params_high (k x n_regressors, rows
= equations, columns [const, ect, lagged diffs]) with EICKER-WHITE
bse_low/bse_high, n_low/n_high/nobs/frac_low, sigma, log_det_sigma,
llf, llf_linear, beta_linear, beta_grid, ect, min_regime, neqs,
n_regressors, k_ar_diff.
Further arguments, with defaults: `n_grid_gamma` (300).
hansen_seo_test¶
def hansen_seo_test(
data: _ArrayLike,
k_ar_diff: int = ...,
trim: float = ...,
n_grid: int = ...,
n_boot: int = ...,
seed: int = ...,
beta: _ArrayLike | None = ...,
) -> dict[str, Any]:
Hansen-Seo (2002) sup-LM test of linear vs two-regime THRESHOLD cointegration, p-valued by their fixed-regressor bootstrap. beta is fixed at the null (linear Johansen ML) estimate unless supplied. The threshold is unidentified under the null (Davies problem), so NO chi-squared p-value exists; p_value = (1 + #{LM* >= LM})/(n_boot + 1) — seeded, parallel, bit-identical at any thread count. Presumes the series ARE cointegrated (test that first: johansen/engle_granger). Unlike threshold_vecm — whose beta=None grid SEARCH is bivariate-only — the null beta here is the linear Johansen ML estimate, defined for any k, so k > 2 with an estimated (null) beta is accepted.
Keys: stat, p_value, threshold, beta, n_boot, nobs, thresholds,
lm_path, boot_stats, min_regime, neqs, n_regressors, k_ar_diff.
Further arguments, with defaults: `trim` (0.05), `n_grid` (300), `seed`
(0).
ou_fit¶
Ornstein-Uhlenbeck mean-reversion fit for a spread (exact-discretization MLE).
``dX = kappa (mu - X) dt + sigma dW`` observed at step ``dt`` is exactly
the AR(1) ``X_{t+1} = c + phi X_t + eps`` with ``phi = exp(-kappa dt)``;
the MLE is the closed-form AR(1) OLS (statsmodels ``AutoReg(x, lags=1)``)
mapped back — no optimizer. Returns ``kappa``/``mu``/``sigma`` with
delta-method ``kappa_se``/``mu_se``/``sigma_se``; ``half_life``
(= ln 2 / kappa) with ``half_life_ci`` at ``level`` — the level-scale
kappa interval mapped through ln 2 / kappa, with an ``inf`` upper
endpoint when the interval crosses zero (Monte-Carlo measured against
the log-scale alternative and shipped because it covers closer to
nominal in every cell; the cointegration model card has the table and
quantifies the well-known upward finite-sample kappa bias of roughly
4 / (sample time span)); ``stationary_sd`` (= sigma / sqrt(2 kappa));
``mean_reverting``; the AR(1) leg ``phi``/``phi_se``/``c``/``c_se``/
``eta2``/``loglik``; and the echoed call inputs ``n_obs``/``dt``/
``level`` (the confidence level ``half_life_ci`` was built at, echoed
back so a stored result dict stays self-describing). A fit with
``phi >= 1`` is returned honestly: ``mean_reverting=False``,
``half_life=inf``, ``half_life_ci=None``, ``stationary_sd=None``.
spread_zscore¶
def spread_zscore(
x: _ArrayLike,
kappa: float | None = ...,
mu: float | None = ...,
sigma: float | None = ...,
dt: float | None = ...,
) -> dict[str, Any]:
Z-score of a spread against the stationary OU law N(mu, sigma^2/(2 kappa)).
``zscore = (x - mu) / stationary_sd``, ``stationary_sd = sigma /
sqrt(2 kappa)``. Pass all three of ``kappa``/``mu``/``sigma`` (a frozen
``ou_fit``) or none (fitted from ``x`` at step ``dt``); partial
specification is refused, and ``dt`` (default 1.0 when None)
parameterizes only the internal fit, so passing it explicitly with a
frozen triple raises. Returns ``zscore``, the ``kappa``/``mu``/
``sigma`` used, ``stationary_sd``, ``fitted``. Refuses ``kappa <= 0``:
a non-mean-reverting process has no stationary distribution to score
against.
regime switching¶
markov_switching_ar¶
def markov_switching_ar(
y: _ArrayLike,
k_regimes: int = ...,
order: int = ...,
switching_variance: bool = ...,
max_iter: int = ...,
tol: float = ...,
) -> dict[str, Any]:
Markov-switching AR fitted by EM (Hamilton 1989); regimes + durations.
ar is the estimated common AR block (phi_1, .., phi_p), a length-order
array shared across regimes (the binding fits Hamilton's common-AR
specification on deviations y_t - mu_{S_t}). smoothed_prob /
filtered_prob are the full (n, k_regimes) probability matrices,
n = len(y) - order; smoothed_prob_last_regime keeps the 0.2.0 scalar
path (= smoothed_prob[:, -1]). The EM run reports converged and
iterations (EM steps actually run — converged=False with
iterations == max_iter means the cap bound).
Returned keys: `ar`, `converged`, `expected_durations`, `filtered_prob`,
`iterations`, `loglik`, `means`, `regimes`, `smoothed_prob`,
`smoothed_prob_last_regime`, `transition`, `variances`.
Further arguments, with defaults: `switching_variance` (True), `tol`
(1e-06).
setar¶
def setar(
y: _ArrayLike,
p: int,
delay: int = ...,
trim: float = ...,
delays: Sequence[int] | None = ...,
ic: str = ...,
constant: bool = ...,
) -> dict[str, Any]:
Two-regime SETAR(p) (Tong-Lim 1980) by concentrated least squares
(Hansen 1997): grid over the trimmed order statistics of y_{t-delay},
per-candidate OLS in each regime, pooled-SSR-minimizing threshold (and
delay, when delays is a list — all candidates then share the common
sample t >= max(p, max(delays)) so SSRs are comparable; delays
overrides delay).
Keys: threshold, delay, params_low/params_high (constant first) with
classical nonrobust bse_low/bse_high, n_low/n_high/nobs, pooled ssr and
sigma2 = SSR/(nobs - 2k), sigma2_low/sigma2_high, aic/bic (n ln(SSR/n) +
penalty * m, m = 2k + 1 counting the threshold), ic/ic_used (`ic`
selects which criterion is *reported* — with p fixed the SSR ranking and
the IC ranking coincide), min_regime, k, and the candidate grid
thresholds with its ssr_path. Validated against an independent NumPy
transcription of the published algorithm (fixtures/setar.json).
Further arguments, with defaults: `trim` (0.15).
setar_test¶
def setar_test(
y: _ArrayLike,
p: int,
delay: int = ...,
trim: float = ...,
n_boot: int = ...,
seed: int = ...,
) -> dict[str, Any]:
Hansen (1996) sup-F linearity test against a two-regime SETAR(p): stat = nobs (ssr_linear - ssr_setar)/ssr_setar over the trimmed threshold grid. The threshold is unidentified under the null (Davies problem), so NO chi-squared p-value exists; p_value = (1 + #{F >= F}) / (n_boot + 1) from the fixed-regressor wild bootstrap (y = resid * eta, eta iid N(0,1), same fixed regressors, same grid) — seeded, parallel, bit-identical at any thread count.
Keys: stat, p_value, threshold, delay, n_boot, nobs, ssr_linear,
ssr_setar, thresholds, f_path, boot_stats.
Further arguments, with defaults: `trim` (0.15), `seed` (0).
star¶
def star(
y: _ArrayLike,
p: int,
model: str = ...,
delay: int = ...,
trim: float = ...,
delays: Sequence[int] | None = ...,
constant: bool = ...,
n_gamma: int = ...,
n_c: int = ...,
) -> dict[str, Any]:
Smooth-transition AR (Terasvirta 1994): y_t = phi1'x_t + G(gamma, c; y_{t-delay}) phi2'x_t + e_t, model "lstar" (G = 1/(1+exp(-gamma(s-c)))) or "estar" (G = 1 - exp(-gamma(s-c)^2)). gamma is RAW (tsDyn convention); gamma_standardized is Terasvirta's gammasd(s) (lstar) / gammavar(s) (estar). Concentrated NLS: (gamma, c) grid (standardized gamma log-spaced [0.5, 100]; c on trimmed order statistics) + Nelder-Mead refinement. Run star_test first — STAR on linear data happily "finds" a transition.
Keys: model, gamma, gamma_standardized, c, delay, s_sd, converged,
gamma_at_boundary (standardized gamma at the searchable range's edge:
read gamma as a bound, not a point estimate), params_linear,
params_nonlinear (high-regime coefficients are the sum), bse_linear,
bse_nonlinear, se_gamma, se_c, se_valid (Gauss-Newton; NaN + False on
a degenerate J'J), ssr, sigma2, loglik, aic, bic, nobs, k, transition,
grid_gamma, grid_c, ssr_grid (n_gamma x n_c), best_cell, fevals.
Further arguments, with defaults: `trim` (0.15), `delays` (None),
`constant` (True).
star_eval¶
def star_eval(
y: _ArrayLike,
p: int,
gamma: float,
c: float,
model: str = ...,
delay: int = ...,
constant: bool = ...,
) -> dict[str, Any]:
Concentrated STAR fit at FIXED (gamma, c) (raw gamma, as in star): OLS of y_t on [x_t, G_t x_t] with Gauss-Newton SEs — for scoring a published parameterization (SSR/loglik comparison is robust to optimizer differences).
Keys: params_linear, params_nonlinear, bse_linear, bse_nonlinear,
se_gamma, se_c, se_valid, ssr, sigma2, loglik, aic, bic, nobs, k,
transition.
Further arguments, with defaults: `model` ("lstar"), `delay` (1),
`constant` (True).
star_test¶
def star_test(
y: _ArrayLike,
p: int,
delay: int = ...,
delays: Sequence[int] | None = ...,
) -> dict[str, Any]:
Terasvirta STAR modeling-cycle battery: LM3 linearity test against
STAR (Luukkonen-Saikkonen-Terasvirta 1988; chi-squared form lm3_stat
with 3q df and small-sample F form lm3_f_stat — no bootstrap needed,
the auxiliary regression is linear so the null is standard, unlike
setar_test) plus the H03/H02/H01 sequence: suggested = "estar" iff
the H02 p-value is strictly smallest. delays evaluates each
candidate delay; best indexes the smallest F-form LM3 p-value and the
top-level scalars are that battery's.
Keys: delay, nobs, q, k0, lm3_stat, lm3_p_value, lm3_f_stat,
lm3_f_p_value, h3_f_stat, h3_p_value, h2_f_stat, h2_p_value,
h1_f_stat, h1_p_value, ssr0, ssr1, ssr2, ssr3, suggested, best,
tests.
threshold_var¶
def threshold_var(
data: _ArrayLike,
p: int,
threshold_index: int = ...,
delay: int = ...,
trim: float = ...,
delays: Sequence[int] | None = ...,
constant: bool = ...,
) -> dict[str, Any]:
Two-regime threshold VAR (the multivariate SETAR) by concentrated
least squares / Gaussian MLE: regime split by z_t =
y[threshold_index]_{t-delay} <= threshold; per-candidate two-regime
OLS minimizing ln det of the pooled residual covariance over the
trimmed order-statistic grid (delays as a list searches the delay
jointly on the common sample and overrides delay). Two regimes
only; regime-dependent generalized IRFs are deliberately not provided
(see the model card).
Keys: threshold, delay, threshold_index, params_low/params_high (k x
n_regressors, rows = equations, columns [const?, y_{t-1}..,
y_{t-p}..]) with classical bse_low/bse_high, n_low/n_high/nobs,
sigma/sigma_low/sigma_high, log_det_sigma, llf, aic/bic, thresholds,
logdet_path, min_regime, neqs, n_regressors.
Further arguments, with defaults: `trim` (0.1), `constant` (True).
threshold_var_test¶
def threshold_var_test(
data: _ArrayLike,
p: int,
threshold_index: int = ...,
delay: int = ...,
trim: float = ...,
n_grid: int = ...,
n_boot: int = ...,
seed: int = ...,
constant: bool = ...,
) -> dict[str, Any]:
Robust sup-Wald (score-form) linearity test of a linear VAR(p) against the two-regime threshold VAR — the multivariate analogue of the Hansen-Seo sup-LM — with a Hansen (1996) fixed-regressor wild bootstrap p-value. The threshold is unidentified under the null (Davies problem), so NO chi-squared p-value exists; p_value = (1 + #{W* >= W})/(n_boot + 1) — seeded, parallel, bit-identical at any thread count.
Keys: stat, p_value, threshold, delay, threshold_index, n_boot, nobs,
thresholds, wald_path, boot_stats, min_regime, neqs, n_regressors.
Further arguments, with defaults: `trim` (0.1), `n_grid` (300), `seed`
(0), `constant` (True).
MIDAS¶
midas_weights¶
MIDAS weights (sum to 1); scheme "exp_almon" or "beta".
umidas¶
def umidas(
y: _ArrayLike, hf_lags: _ArrayLike, se_type: str = ..., maxlags: int | None = ...
) -> dict[str, Any]:
U-MIDAS: unrestricted mixed-frequency regression (hf_lags is nobs x K).
Returned keys: `bse`, `params`, `rsquared`.
Further arguments, with defaults: `se_type` ("hac"), `maxlags` (None).
multivariate GARCH¶
ccc_garch¶
def ccc_garch(
returns: _ArrayLike,
forecast_horizon: int = ...,
vol: str = ...,
mean: str = ...,
univariate_dist: str = ...,
p: int = ...,
o: int | None = ...,
q: int = ...,
) -> dict[str, Any]:
CCC-GARCH (Bollerslev 1990); returns is T x k, H_t = D_t R D_t with
a constant correlation R. The univariate stage defaults to the
historical zero-mean Normal GARCH(1,1) (default calls bit-identical)
and takes garch_fit's knobs: vol/mean/p/o/q plus univariate_dist
("normal" | "t", the per-series innovation density — named apart from
dcc_garch's dist=, which is the second-stage correlation likelihood).
o follows garch_fit's sentinel: None (default) means no asymmetry
term under vol="garch" and one lag under "gjr"/"egarch"; explicit
o > 0 under vol="garch" raises instead of being silently dropped.
Keys: correlation (k x k), loglik, sigma2 ((T, k) per-series
conditional variance paths), covariance ((T, k, k) in-sample
H_t = D_t R D_t), and with forecast_horizon > 0 covariance_forecast
((horizon, k, k)) and variance_forecast ((horizon, k)).
Timing (identical to dcc_garch): sigma2[t] — hence covariance[t] —
conditions on information through t-1 (the standard GARCH filter);
covariance[-1] is the last IN-SAMPLE matrix, and the one-step-ahead
H_{T+1} is covariance_forecast[0]. Because R is constant the forecast
is analytic and exact at every horizon (no DCC-style h >= 2
approximation), with the variance forecasts identical to each series'
own garch_fit forecast path. The univariate stage inherits garch_fit's
EGARCH limit: vol="egarch" with forecast_horizon >= 2 raises.
dcc_garch¶
def dcc_garch(
returns: _ArrayLike,
variant: str = ...,
dist: str = ...,
forecast_horizon: int = ...,
vol: str = ...,
mean: str = ...,
univariate_dist: str = ...,
p: int = ...,
o: int | None = ...,
q: int = ...,
) -> dict[str, Any]:
DCC-GARCH (Engle 2002); returns is T x k. variant: "dcc" | "cdcc" (Aielli 2013 consistent targeting) | "adcc" (Cappiello-Engle-Sheppard 2006 asymmetric); dist: "normal" | "t" (second-stage likelihood; "t" adds nu). The univariate first stage defaults to the historical zero-mean Normal GARCH(1,1) and takes garch_fit's knobs: vol/mean/p/o/q plus univariate_dist ("normal" | "t") — the per-series innovation density, distinct from dist=, which configures the second-stage correlation likelihood; o follows garch_fit's sentinel (None default; explicit o > 0 under vol="garch" raises). Keys: a, b, g, qbar, loglik, converged, variant, dist, correlation ((T, k, k) nested list -- the in-sample conditional correlation path), correlation_last, sigma2 ((T, k) per-series conditional variance paths), covariance ((T, k, k) in-sample H_t = D_t R_t D_t), univariate (a list of k dicts -- the full per-series stage-1 GARCH results with exactly garch_fit's keys, bit-identical to garch_fit on that column under the same spec), std_residuals ((T, k) -- the stacked stage-1 standardized residuals z[t][i] = eps_{i,t} / sqrt(sigma2[t][i]) that drive the correlation recursion), nu (dist="t" only), nbar (variant="adcc" only -- the (k, k) asymmetric targeting matrix Nbar = mean of n_t n_t', n_t = min(z_t, 0)), and with forecast_horizon > 0 correlation_forecast / covariance_forecast ((horizon, k, k)) and variance_forecast ((horizon, k)).
Timing: correlation[t] = R_t conditions on information through t-1
(filter convention; Q_0 = Qbar), and sigma2[t] / covariance[t] follow
the same convention. correlation_last = correlation[-1] is
the last IN-SAMPLE matrix, not a forecast; the one-step-ahead R_{T+1}
also uses the final residual z_T and is correlation_forecast[0].
h >= 2 forecasts use the Engle-Sheppard (2001) recursion on E[Q]
normalized each step (an approximation), converging to corr(qbar).
The univariate stage inherits garch_fit's EGARCH limit: vol="egarch"
with forecast_horizon >= 2 raises. The default call is bit-identical to
earlier releases.
dcc_test¶
def dcc_test(
returns: _ArrayLike,
lags: int = ...,
vol: str = ...,
mean: str = ...,
univariate_dist: str = ...,
p: int = ...,
o: int | None = ...,
q: int = ...,
) -> dict[str, Any]:
Engle-Sheppard (2001) test of constant conditional correlation (CCC vs DCC); returns is T x k. A univariate GARCH per series (the same first stage ccc_garch/dcc_garch use — zero-mean Normal GARCH(1,1) by default, configurable via vol/mean/univariate_dist/p/o/q; o follows garch_fit's sentinel — None default, explicit o > 0 under vol="garch" raises), joint standardization by the symmetric inverse square root of the constant correlation, pooled AR(lags) on the stacked off-diagonal outer products. Keys: stat, df (= lags + 1), p_value (small rejects constant correlation), lags, nobs, n_stacked.
realized volatility / HAR¶
realized_measures¶
Realized variance, bipower variation, and jump component (BNS 2004).
Returned keys: `bipower`, `jump`, `rv`.
har_rv¶
def har_rv(
rv: _ArrayLike,
start: int = ...,
variant: str = ...,
hac_maxlags: int = ...,
use_correction: bool = ...,
) -> dict[str, Any]:
HAR-RV (Corsi 2009): RV_t on [const, daily, weekly, monthly], HAC SEs.
The aggregates follow Corsi's definition and INCLUDE the daily lag:
weekly = mean(RV_{t-1}..RV_{t-5}), monthly = mean(RV_{t-1}..RV_{t-22}).
(Changed in 0.5: through 0.4.0 the windows mistakenly excluded RV_{t-1};
coefficients on the same data shift.)
variant is "level", "log", or "sqrt". use_correction now defaults True
(False through 0.2.0): bse/tvalues carry the finite-sample sqrt(n/(n-k))
factor by default. statsmodels cov_type="HAC" defaults the correction
off -- pass use_correction=False to match it (and the old numbers).
Returned keys: `bse`, `nobs`, `params`, `rsquared`, `tvalues`.
Further arguments, with defaults: `start` (22), `hac_maxlags` (5).
connectedness¶
connectedness¶
def connectedness(
data: _ArrayLike, lags: int = ..., horizon: int = ..., trend: str = ...
) -> dict[str, Any]:
Diebold-Yilmaz connectedness (percent) from a VAR's GFEVD.
total, to_others, from_others, net, gfevd, pairwise_net (data is T x k).
Further arguments, with defaults: `lags` (2), `horizon` (10), `trend`
("c").
Returned keys: `from_others`, `gfevd`, `net`, `pairwise_net`,
`to_others`, `total`.
factor model¶
factor_model¶
PCA factor model (T x N) + Bai-Ng (2002) factor selection.
factors, loadings, eigenvalues, icp1/icp2/pcp1/pcp2 and the
Ahn-Horenstein eigenvalue-ratio factor count `er` with the ratios
`er_ratios` (length kmax) it was read from.
Further arguments, with defaults: `n_factors` (2).
Returned keys: `eigenvalues`, `er`, `er_ratios`, `factors`, `icp1`,
`icp2`, `loadings`, `pcp1`, `pcp2`.
term structure¶
nelson_siegel¶
def nelson_siegel(
maturities: _ArrayLike,
yields: _ArrayLike,
decay: float = ...,
optimal_lambda: bool = ...,
) -> dict[str, Any]:
Nelson-Siegel yield-curve fit (Diebold-Li 2006).
level/slope/curvature factors, lambda, residuals, rsquared.
optimal_lambda=True estimates the decay by NLS.
Returned keys: `curvature`, `factors`, `lambda`, `level`, `residuals`,
`rsquared`, `slope`.
svensson¶
def svensson(
maturities: _ArrayLike, yields: _ArrayLike, lambda1: float, lambda2: float
) -> dict[str, Any]:
Svensson (1994) four-factor yield-curve fit; nests Nelson-Siegel.
Returned keys: `factors`, `lambda1`, `lambda2`, `residuals`, `rsquared`.
GMM / IV-GMM¶
iv_gmm¶
def iv_gmm(
x: _ArrayLike,
z: _ArrayLike,
y: _ArrayLike,
method: str = ...,
weight: str = ...,
bandwidth: float | None = ...,
tol: float = ...,
max_iter: int = ...,
) -> dict[str, Any]:
Linear IV-GMM (Hansen 1982) with robust or HAC weighting.
POSITIONAL ORDER IS (x, z, y): regressors, instruments, outcome. x and z
are both 2-D float matrices, so swapping them coerces cleanly and returns
plausible-looking garbage -- prefer keywords: iv_gmm(x=X, z=Z, y=y).
bandwidth defaults to None, which selects the Newey-West rule of thumb.
It previously defaulted to 0.0 -- a Bartlett kernel truncated at zero
lags IS White, so weight="hac" used to be a silent no-op returning
results bit-identical to weight="robust". An explicit bandwidth=0.0 now
raises. The truncation actually used comes back as hac_bandwidth.
Neither setting restores nominal coverage under persistent moments: the
audit measured 0.868 against nominal 0.95 at bandwidth=10.
Also returns first_stage, a list of per-regressor weak-instrument F
diagnostics keyed by "regressor". Entries are omitted where the
statistic is undefined, so it may be shorter than the regressor count,
and a missing entry is not a failed fit. With two or more endogenous
regressors these are NOT a weak-identification test -- all can clear 10
while the system is under-identified. F > 10 is not a safety threshold
even with one: coverage was 0.915 at a median F of 10.5.
method is "2sls", "2step", or "iterated"; weight is "robust" or "hac".
Z must include the exogenous regressor columns. Returns params, bse, cov,
residuals, nobs, nmoments, nparams, steps, hac_bandwidth, first_stage,
and (over-identified) the Hansen j_stat/j_dof/j_pval.
Further arguments, with defaults: `tol` (1e-08), `max_iter` (100).
Returned keys: `bse`, `cov`, `first_stage`, `hac_bandwidth`, `j_dof`,
`j_pval`, `j_stat`, `nmoments`, `nobs`, `nparams`, `params`,
`residuals`, `steps`.
leakage-safe time-series CV¶
cv_splits¶
def cv_splits(
n: int,
scheme: str = ...,
train: int = ...,
horizon: int = ...,
step: int = ...,
k: int = ...,
purge: int = ...,
embargo: int = ...,
) -> list[dict[str, list[int]]]:
Leakage-safe CV split indices for sequential data.
scheme is "expanding", "rolling", or "purged_kfold". Returns a list of
{"train": [...], "test": [...]} index dicts. `train` defaults to 0, which
the two walk-forward schemes refuse as an empty first window — pass it
explicitly there ("purged_kfold" ignores it). purge drops the last purge
indices from the end of every training window (all schemes; set it >=
horizon - 1 for h-step-ahead labels). embargo excludes training rows
after the test block, which only exist under "purged_kfold"; nonzero
embargo raises on "expanding"/"rolling". Under "purged_kfold" the
embargo is measured from the end of the purged window (Lopez de Prado
2018, ch. 7), so the right-hand gap is purge + embargo indices.
Further arguments, with defaults: `k` (5).
penalized ML (paths)¶
adaptive_lasso¶
def adaptive_lasso(
x: _ArrayLike,
y: _ArrayLike,
alpha: float,
l1_ratio: float = ...,
gamma: float = ...,
tol: float = ...,
max_iter: int = ...,
) -> dict[str, Any]:
Adaptive LASSO (Zou 2006): oracle-property weighted-L1 penalty.
Keys: `coef`, `n_iter`, `max_change` and `max_rel_change` (as in
`elastic_net`: the scale-free update the stopping rule compares with
`tol`).
Further arguments, with defaults: `l1_ratio` (1.0), `gamma` (1.0),
`max_iter` (100000).
lasso_path¶
def lasso_path(
x: _ArrayLike,
y: _ArrayLike,
l1_ratio: float = ...,
n_lambdas: int = ...,
eps: float = ...,
tol: float = ...,
max_iter: int = ...,
) -> dict[str, Any]:
Elastic-net regularization path with AIC/BIC selection.
lambdas, coefs, rss, df, aic, bic, aic_best, bic_best.
Further arguments, with defaults: `l1_ratio` (1.0), `n_lambdas` (100),
`eps` (0.001), `tol` (1e-07), `max_iter` (100000).
Returned keys: `aic`, `aic_best`, `bic`, `bic_best`, `coefs`, `df`,
`lambdas`, `rss`.
forecast backtest¶
backtest¶
def backtest(
y: _ArrayLike,
window: str = ...,
train: int = ...,
horizon: int = ...,
refit_every: int = ...,
forecaster: str | Callable[[_F64, int], _ArrayLike] | None = ...,
period: int | None = ...,
insample_period: int = ...,
) -> dict[str, Any]:
Rolling/expanding pseudo-out-of-sample backtest.
window is "expanding" or "rolling"; forecaster is one of naive, drift,
mean, seasonal_naive, theta (None means naive) — or any Python callable.
period (default 1 when None) is the seasonal cycle length of
seasonal_naive/theta only; passing it explicitly with any other
forecaster (callables included) raises — the MASE/RMSSE scale period
is the separate insample_period. forecaster may be any Python callable
f(train, horizon) -> array-like of exactly `horizon` finite point
forecasts, where train is a read-only float64 ndarray holding only the
training window for that origin (the engine's leakage discipline; with
refit_every > 1 the callable is asked for up to refit_every - 1 + horizon
steps at each refit origin). Exceptions raised inside the callable
re-raise naming the failing origin and window, with the original chained
as __cause__; wrong-length / non-finite returns raise teaching errors.
Returns origins, per-horizon forecasts and
targets, and a per-horizon accuracy table.
Returned keys: `accuracy`, `forecasts`, `horizon`, `n_origins`,
`origins`, `targets`.
conformal forecast intervals¶
conformal_forecast¶
def conformal_forecast(
y: _ArrayLike,
horizon: int = ...,
method: str = ...,
base: str | Callable[[_F64, int], _ArrayLike] | None = ...,
alpha: float = ...,
calib: int | None = ...,
mode: str = ...,
period: int = ...,
gamma: float | None = ...,
n_eval: int | None = ...,
lags: int | None = ...,
n_boot: int | None = ...,
seed: int | None = ...,
optimize_beta: bool | None = ...,
order: tuple[int, int, int] | None = ...,
) -> dict[str, Any]:
Distribution-free conformal forecast intervals around a point forecaster.
method is "split" (finite-sample-corrected residual-quantile calibration
on held-out origins; mode "symmetric" or "asymmetric"), "enbpi" (Xu-Xie
2021 bootstrap-ensemble batch prediction intervals; base must be "ar"),
or "aci" (Gibbs-Candes 2021 adaptive conformal inference,
alpha_{t+1} = alpha_t + gamma (alpha - err_t), gamma default 0.005 from
the paper). base wraps "theta" (None means "theta"), "naive", "drift",
"mean", "seasonal_naive", "ar", or "arima" (order=(p, d, q)) — or, for
split/aci, any Python callable base(train, horizon) -> array-like of
horizon point forecasts with the backtest contract (train is a read-only
float64 ndarray of the training window only). calib defaults to
n // 4 residuals per horizon; n_eval (aci) to n // 5. `seed=None` (the
default) is NOT fresh entropy — under "enbpi" it means seed 0, so two
default calls are bit-identical (`n_boot=None` means 25). Inert kwargs
raise (0.7.0; defaults bit-identical): order needs base="arima"; lags
needs base="ar" (or EnbPI's ensemble); gamma needs method="aci";
n_boot/seed/optimize_beta need method="enbpi"; n_eval is ACI-only in
this function (conformal_backtest uses it for every method); calib
never reaches EnbPI. Returns mean,
lower, upper, level, plus per-method calibration diagnostics (split:
q_lower/q_upper/scores/finite_sample_level; enbpi: beta/residuals;
aci: alpha_final/alpha_trajectory/err/realized_coverage).
Returned keys: `alpha`, `base`, `finite_sample_level`, `horizon`,
`level`, `lower`, `mean`, `method`, `mode`, `n_calib`, `q_lower`,
`q_upper`, `scores`, `upper`.
Further arguments, with defaults: `period` (1).
conformal_backtest¶
def conformal_backtest(
y: _ArrayLike,
horizon: int = ...,
method: str = ...,
base: str | Callable[[_F64, int], _ArrayLike] | None = ...,
alpha: float = ...,
calib: int | None = ...,
mode: str = ...,
period: int = ...,
gamma: float | None = ...,
n_eval: int | None = ...,
lags: int | None = ...,
n_boot: int | None = ...,
batch: int | None = ...,
seed: int | None = ...,
optimize_beta: bool | None = ...,
order: tuple[int, int, int] | None = ...,
) -> dict[str, Any]:
Online out-of-sample evaluation of conformal intervals ("split",
"aci", or "enbpi") over the last n_eval origins: per-origin intervals
formed from information available then, miss indicators, and realized
coverage per horizon. base as in conformal_forecast, including a Python
callable base(train, horizon) for split/aci. ACI adds its alpha_t
trajectory; EnbPI is the published one-step online algorithm with the
residual window sliding by batch. seed=None (the default) means seed 0
under "enbpi", not fresh entropy (n_boot=None means 25). The same inert-kwarg refusals as
conformal_forecast apply (order/lags/gamma/n_boot/seed/optimize_beta,
EnbPI's calib), plus batch, which is EnbPI-only; n_eval is live for
every method here. Defaults stay bit-identical.
Returned keys: `alpha`, `base`, `err`, `horizon`, `level`, `lower`,
`mean`, `method`, `n_eval`, `origins`, `realized_coverage`, `upper`.
Further arguments, with defaults: `alpha` (0.1), `mode` ("symmetric"),
`period` (1).
nonlinear GMM (callback)¶
gmm_nonlinear¶
def gmm_nonlinear(
moments_fn: Callable[[_F64], _ArrayLike],
initial: _ArrayLike,
weight: _ArrayLike | None = ...,
) -> dict[str, Any]:
Nonlinear GMM (Hansen 1982) via Nelder-Mead over a Python moment function.
moments_fn maps a parameter vector (a 1-D float64 array) to an n-by-m matrix
of per-observation moment contributions (rows = observations, cols = moments),
returned as a NumPy array or list of lists -- the return must be 2-D even
for a single moment condition (reshape with g.reshape(-1, 1)); a 1-D return
raises a TypeError naming moments_fn. weight is the flattened m*m
weighting matrix (row-major) or None for the identity. Returns params,
objective, gbar, converged, iterations, fevals, nmoments, nparams.
Returned keys: `converged`, `fevals`, `gbar`, `iterations`, `nmoments`,
`nparams`, `objective`, `params`.
weighted MIDAS¶
weighted_midas¶
def weighted_midas(
y: _ArrayLike,
hf_lags: _ArrayLike,
scheme: str = ...,
weight_start: tuple[float, float] | None = ...,
) -> dict[str, Any]:
Weighted MIDAS by NLS (Ghysels et al. 2007); exp_almon/beta weights, hf_lags is nobs x K.
Returned keys: `converged`, `fitted`, `intercept`, `iterations`,
`residuals`, `rsquared`, `scheme`, `slope`, `ssr`, `weight_params`,
`weights`.
Further arguments, with defaults: `scheme` ("exp_almon"), `weight_start`
(None).
state-dependent LP¶
lp_state¶
def lp_state(
y: _ArrayLike,
shock: _ArrayLike,
state_indicator: _ArrayLike,
horizons: int = ...,
n_lag_controls: int = ...,
se: str | None = ...,
maxlags: int | None = ...,
cumulative: bool | str | None = ...,
band: str | None = ...,
band_alpha: float = ...,
) -> dict[str, Any]:
State-dependent (interacted) local projections (Ramey-Zubairy 2018); per-regime IRFs and SEs.
`cumulative` takes False/"none", True/"outcome" or "both", as in `lp`.
`se=None` (the default) resolves to "lag_augmented" — except under
`cumulative="both"`, where it resolves to "hac" for the same reason as in
`lp` (the cumulated impulse shares future shocks across nearby base
times, so lag-augmented HC1 is inconsistent there; audit: 0.640 coverage
at a nominal 95%, h=12) — and `se="lag_augmented"` with
`cumulative="both"` raises. The method actually used is returned as
`se_method`.
**Bands.** `band=None` (default) returns no band. `"pointwise"`, `"sidak"`
and `"bonferroni"` add one band PER REGIME —
`lower_state1`/`upper_state1` and `lower_state0`/`upper_state0`, with
`critical_value_state1`/`critical_value_state0`,
`n_cells_used_state1`/`n_cells_used_state0` and
`cov_se_max_rel_diff_state1`/`cov_se_max_rel_diff_state0` (always None
here: no covariance is built) — over the horizons of that
regime's own response (`K = horizons + 1`, `band_scope="horizon"`). The two
regimes are banded separately; nothing here is simultaneous *across*
regimes.
`band="sup-t"` is **refused**: no cross-horizon covariance is estimated for
the interacted regressions, so `lp_state` (like `lp_iv` and `lp_multiplier`)
gets the closed-form simultaneous routes only. Report such a band as Šidák
or Bonferroni, never as sup-t.
Returned keys: `horizons`, `irf_state0`, `irf_state1`, `se_method`,
`se_state0`, `se_state1`.
Further arguments, with defaults: `n_lag_controls` (4), `maxlags`
(None), `band_alpha` (0.1).
mean-group panel VAR¶
mean_group_var¶
def mean_group_var(
entities: Sequence[_ArrayLike],
lags: int = ...,
trend: str = ...,
horizon: int = ...,
response: int = ...,
impulse: int = ...,
) -> dict[str, Any]:
Pesaran-Smith (1995) mean-group panel VAR over per-entity T_i x k matrices.
Returned keys: `coefs`, `coefs_se`, `intercept`, `intercept_se`,
`irf_path`, `irf_path_se`, `lags`, `n_entities`, `neqs`, `orth_irfs`,
`orth_irfs_se`.
Further arguments, with defaults: `lags` (1), `trend` ("c"), `horizon`
(10), `response` (0), `impulse` (0).
dynamic Nelson-Siegel¶
dynamic_ns¶
Dynamic Nelson-Siegel factors + one-step forecast (Diebold-Li 2006).
panel is T x n_maturities. Returns maturities, lambda, factors (T x 3),
rsquared, level/slope/curvature series, and a forecast dict.
Further arguments, with defaults: `decay` (0.0609).
Returned keys: `curvature`, `factors`, `forecast`, `lambda`, `level`,
`maturities`, `rsquared`, `slope`.
FAVAR¶
favar¶
def favar(
panel: _ArrayLike,
policy: _ArrayLike,
n_factors: int = ...,
lags: int = ...,
trend: str = ...,
slow_indices: list[int] | None = ...,
horizon: int = ...,
orth: bool = ...,
) -> dict[str, Any]:
Two-step factor-augmented VAR (Bernanke-Boivin-Eliasz 2005).
factors (T x r), params, sigma_u, n_factors, n_endog, policy_index, and
the recursive policy-shock IRFs irf_panel (N x horizon+1) / irf_policy.
Further arguments, with defaults: `lags` (2), `trend` ("c"),
`slow_indices` (None), `orth` (True).
Returned keys: `factors`, `irf_panel`, `irf_policy`, `n_endog`,
`n_factors`, `params`, `policy_index`, `sigma_u`.
realized-volatility extras¶
realized_quarticity¶
Realized quarticity RQ = (n/3) sum r^4 (BNS 2002).
tripower_quarticity¶
Jump-robust tripower quarticity of integrated quarticity (BNS 2004).
bns_jump_test¶
BNS ratio jump test in the Huang & Tauchen (2005) form; dict with 'ratio'.
The HT finite-sample M/(M-1) and M/(M-2) scalings on BV/TQ are applied
inside the statistic; the exported measures stay unadjusted BNS 2004.
Returned keys: `ratio`.
realized_range¶
def realized_range(
high: _ArrayLike,
low: _ArrayLike,
method: str = ...,
open: _ArrayLike | None = ...,
close: _ArrayLike | None = ...,
) -> float:
Range variance from OHLC bars; method is "parkinson" or "garman_klass".
Further arguments, with defaults: `open` (None), `close` (None).
score-driven models (GAS/DCS)¶
gas_volatility¶
GAS(1,1) score-driven volatility (Creal-Koopman-Lucas 2013).
density is "gaussian" or "student_t". Returns omega/a/b (+ nu),
variance, std_resid, loglik, aic, bic, next_variance, `converged` and
`iterations` (the optimizer's certificate and step count), and
(horizon>0) a forecast.
Returned keys: `a`, `aic`, `b`, `bic`, `converged`, `forecast`,
`iterations`, `loglik`, `next_variance`, `omega`, `std_resid`,
`variance`.
dcs_local_level¶
DCS robust local level mu_{t+1} = mu_t + kappa*u_t (Harvey-Luati 2014).
MLE of (kappa, scale[, nu]). density is "t" (default; bounded redescending
score — robust to additive outliers), "laplace" (sign filter, tracks a
local median), or "gaussian" (exactly the steady-state Kalman local level;
kappa = steady-state gain). Returns kappa/scale (+ nu) with
observed-information *_se, the one-step-predicted level path, resid,
next_level, loglik, aic, bic, honest converged, iterations, n_obs,
density.
Returned keys: `aic`, `bic`, `converged`, `density`, `iterations`,
`kappa`, `kappa_se`, `level`, `loglik`, `n_obs`, `next_level`, `nu`,
`nu_se`, `resid`, `scale`, `scale_se`.
heterogeneous panel (MG)¶
panel_mean_group¶
def panel_mean_group(
ys: Sequence[_ArrayLike], xs: Sequence[_ArrayLike], method: str = ...
) -> dict[str, Any]:
Mean-group (Pesaran-Smith 1995) / CCE-MG (Pesaran 2006) panel estimator.
method is "mg" or "cce". ys/xs are per-unit response vectors and T_i x k
regressor matrices. Returns coef, se, tstat, coef_per_unit, n_units, k.
Returned keys: `coef`, `coef_per_unit`, `k`, `n_units`, `se`, `tstat`.
panel_pmg¶
def panel_pmg(
ys: Sequence[_ArrayLike],
xs: Sequence[_ArrayLike],
tol: float = ...,
max_iter: int = ...,
) -> dict[str, Any]:
Pooled Mean Group ARDL(1,1) panel estimator (Pesaran-Shin-Smith 1999).
Pools the long-run coefficient across units by ML; error-correction speed
and short-run dynamics stay unit-specific. Returns theta, theta_se,
phi_bar, phi, sigma2, loglik, iterations, n_units, k.
tol (default 3e-13) is the RELATIVE stopping tolerance of the
back-substitution: stop when |dtheta|_inf <= tol * (1 + |theta|_inf),
within max_iter (default 1000) iterations per pass. Relative because
the float noise floor of the pooled solve scales with |theta|; 3e-13
is the historical absolute 1e-12 rule's measured effective relative
stringency on the O(1)-theta panels it was validated on, so the
golden fixture stops at the identical iterate (bit-identical). The
iteration runs from the deterministic theta = 0 start and, when that
start diverges (routine on textbook I(1) panels — measured 14/20
seeds of a stable I(1) battery pre-fix, 0/20 post-fix), it is rerun
once from the Pesaran-Shin-Smith unrestricted-ARDL start.
Non-convergence from both starts raises (the last iterate is not a
verified fixed point); raise max_iter or loosen tol for genuinely
slow-mixing panels.
Returned keys: `iterations`, `k`, `loglik`, `n_units`, `phi`, `phi_bar`,
`sigma2`, `theta`, `theta_se`.
panel_unit_root¶
def panel_unit_root(
data: _ArrayLike | Sequence[_ArrayLike],
test: str = ...,
lags: str | int | None = ...,
regression: str = ...,
max_lags: int | None = ...,
lrv_kernel: str = ...,
lrv_bandwidth: float | None = ...,
) -> dict[str, Any]:
First-generation panel unit-root tests (LLC, IPS, Fisher/Maddala-Wu-Choi).
data is a balanced N x T array (rows = units) or a list of 1-D per-unit
series (unbalanced OK for "ips"/"fisher"; "llc" needs a common T). test is
"ips" (default), "llc", or "fisher"; regression is "c"/"ct"/"n" ("n" is
invalid for "ips"); lags is None (per-unit auto AIC), an int (fixed common
lag), or "aic"/"bic"/"t-stat". Returns statistic, p_value,
per_unit_tstat/pvalue/lags/nobs, n_units, regression, plus test-specific
extras: ips -> t_bar; llc -> delta_hat, t_delta, s_n, t_bar_periods;
fisher -> maddala_wu, choi_z, choi_z_pvalue.
Returned keys: `n_units`, `p_value`, `per_unit_lags`, `per_unit_nobs`,
`per_unit_pvalue`, `per_unit_tstat`, `regression`, `statistic`, `t_bar`,
`test`.
Further arguments, with defaults: `max_lags` (None), `lrv_kernel`
("bartlett"), `lrv_bandwidth` (None).
DFM nowcasting¶
dfm_nowcast¶
def dfm_nowcast(
data: _ArrayLike,
n_factors: int = ...,
factor_order: int = ...,
method: str = ...,
) -> dict[str, Any]:
Dynamic-factor-model nowcast; data is T x N with an optional NaN edge.
method is "two_step" (Doz-Giannone-Reichlin 2011) or "mle" (exact
one-step Gaussian MLE, single factor). Returns nowcast, edge_factor,
loglik, fit_loglik, smoothed_factors ((T_b, r): one row per BALANCED-panel
observation, i.e. T minus the ragged-edge rows), n_factors, factor_order,
and the fitted model itself (both methods, same surface): loadings
((N, r)), factor_ar ((r, r*p) stacked [A_1 | ... | A_p]), factor_cov
((r, r)), idiosyncratic (length N), center / scale (length N training
moments; scale is all ones for "mle"). Mapping factors to series
levels is exact: nowcast == center + scale * (loadings @ edge_factor),
and center + scale * (F @ L.T) is the common-component fit of the
balanced panel. Also,
Returned keys: `center`, `edge_factor`, `factor_ar`, `factor_cov`,
`factor_order`, `fit_loglik`, `idiosyncratic`, `loadings`, `loglik`,
`n_factors`, `nowcast`, `scale`, `smoothed_factors`.
dfm_news¶
def dfm_news(
old_vintage: _ArrayLike,
new_vintage: _ArrayLike,
target_series: int = ...,
target_period: int | None = ...,
n_factors: int = ...,
factor_order: int = ...,
) -> dict[str, Any]:
News/update decomposition of a DFM nowcast revision (Banbura-Modugno 2014).
Splits the target-series nowcast revision between two data vintages into
per-datapoint contributions (weight*news). Returns old_nowcast,
new_nowcast, total_revision, and contributions (a list of dicts).
Returned keys: `contributions`, `new_nowcast`, `old_nowcast`,
`target_period`, `target_series`, `total_revision`.
Further arguments, with defaults: `target_series` (0), `target_period`
(None), `n_factors` (1), `factor_order` (2).
predictive regressions / IVX¶
predictive_regression¶
def predictive_regression(
r: _ArrayLike, x: _ArrayLike, cz: float | None = ..., alpha: float = ...
) -> dict[str, Any]:
Predictive regression with a persistent regressor.
Returns ols, stambaugh (bias-corrected), and ivx (Kostakis-Magdalinos-
Stamatogiannis 2015, Wald test valid uniformly over persistence).
Returned keys: `ivx`, `nobs`, `ols`, `stambaugh`.
Further arguments, with defaults: `cz` (None = -1.0), `alpha` (0.95).
ivx_test¶
def ivx_test(
r: _ArrayLike,
xs: _ArrayLike,
cz: float | None = ...,
alpha: float = ...,
joint: str = ...,
) -> dict[str, Any]:
Joint IVX predictability test for several persistent predictors (xs is T x k).
Returns beta_ivx, the joint wald/pvalue, rz, nregressors, nobs. The
default is joint="bonferroni" (changed in 0.5; through 0.4.0 the default
was "chi2"): per-predictor scalar IVX tests combined at level/k, whose
measured size is at or below nominal for every measured k, with power on
par with a size-corrected chi-square test for sparse alternatives. It
adds wald_scalar/pvalue_scalar/joint keys, and its `wald` is the LARGEST
scalar statistic (chi-square(1) scale) with `pvalue` already
Bonferroni-adjusted. The flip is measured, not stylistic: the
joint="chi2" chi-square(k) Wald's size degrades in k near a unit root
(0.28 at k=8, n=250, nominal 0.05) and n does not repair it (alpha=0.5
restores convergence but still ~0.13 at k=8, n=250); chi2 stays
available for small k or rho safely below 1 — see the
predictive-regressions model card.
`cz` (None = -1.0) and `alpha` (0.95) tune the IVX instrument's persistence
`rho_z = 1 + cz / n^alpha` (Kostakis-Magdalinos-Stamatogiannis 2015),
exactly as in `predictive_regression`; `alpha` here is not a
significance level (no level is passed; `pvalue` is returned).
Returned keys: `beta_ivx`, `joint`, `nobs`, `nregressors`, `pvalue`,
`pvalue_scalar`, `rz`, `wald`, `wald_scalar`.
recession probability¶
recession_probit¶
def recession_probit(
y: _ArrayLike, x: _ArrayLike, link: str = ..., dynamic: bool = ...
) -> dict[str, Any]:
Probit/logit of a binary recession indicator (Kauppi-Saikkonen dynamic option).
link is "probit" or "logit". Returns params, bse, zstats, probabilities,
loglik, pseudo_r2, converged (and rho for dynamic=True).
Returned keys: `bse`, `converged`, `loglik`, `params`, `probabilities`,
`pseudo_r2`, `zstats`.
survey expectations¶
cg_regression¶
def cg_regression(
errors: _ArrayLike,
revisions: _ArrayLike,
maxlags: int | None = ...,
use_correction: bool = ...,
) -> dict[str, Any]:
Coibion-Gorodnichenko (2015) information-rigidity regression (OLS-HAC).
Returns `intercept` and `slope` with HAC `se_intercept`/`se_slope`, the
slope's `t_slope`/`p_slope`, `r_squared`, `implied_rigidity`, and the
echoed `maxlags`/`nobs`.
use_correction defaults True (the n/(n-k) HAC scaling); statsmodels
cov_type="HAC" defaults it off -- match it when comparing.
forecast_efficiency¶
def forecast_efficiency(
errors: _ArrayLike,
regressors: _ArrayLike,
maxlags: int | None = ...,
use_correction: bool = ...,
) -> dict[str, Any]:
Mincer-Zarnowitz forecast-efficiency Wald test (OLS-HAC); regressors is T x k.
use_correction defaults True (the n/(n-k) HAC scaling); statsmodels
cov_type="HAC" defaults it off.
Returned keys: `bse`, `params`, `pvalues`, `r_squared`, `tvalues`,
`wald`, `wald_df`, `wald_pvalue`.
Further arguments, with defaults: `maxlags` (None).
forecast_disagreement¶
Forecast-disagreement measures (per-period std/quartiles/iqr) from a forecaster panel.
Returned keys: `counts`, `iqr`, `p25`, `p50`, `p75`, `std`.
Further arguments, with defaults: `ddof` (1).
long memory¶
frac_diff¶
Fractional differencing (1-L)^d via the binomial expansion.
frac_integrate¶
Fractional integration (1-L)^-d, the inverse of frac_diff.
long_memory_d¶
Estimate the memory parameter d; method is "gph" or "local_whittle".
Returns d, se, se_asymptotic and m for both methods, plus se_regression for
method="gph". BUILD INTERVALS FROM `se`: it is the standard error at the
bandwidth actually used. `se_asymptotic` is the textbook large-m closed form
(pi/sqrt(24m) for GPH, 1/(2*sqrt(m)) for local Whittle), kept for reference
-- at the default bandwidth it is materially too NARROW, measured about 25%
at n=512.
Returned keys: `d`, `m`, `se`, `se_asymptotic`, `se_regression`.
specification tests¶
heteroskedasticity_test¶
Heteroskedasticity test (test="white" or "breusch_pagan"); x is T x k with a constant.
Returned keys: `df`, `f_pvalue`, `fstat`, `pvalue`, `statistic`.
reset_test¶
Ramsey RESET functional-form F-test; x is T x k.
Returned keys: `df_den`, `df_num`, `fstat`, `pvalue`.
Further arguments, with defaults: `max_power` (3).
chow_test¶
Chow structural-break F-test at a known 0-indexed split; x is T x k.
Returned keys: `df_den`, `df_num`, `fstat`, `pvalue`, `ssr1`, `ssr2`,
`ssr_pooled`.
cusum_test¶
CUSUM parameter-stability test (Brown-Durbin-Evans); returns the path and 5% bounds.
Returned keys: `bound_lower`, `bound_upper`, `path`, `sigma`.
arbitrage-free NS¶
afns_adjustment¶
Arbitrage-free Nelson-Siegel yield adjustment (Christensen-Diebold-Rudebusch 2011); sigma has 3 elements.
Further arguments, with defaults: `decay` (0.0609).
ACM term premium¶
acm_term_premium¶
def acm_term_premium(
yields: _ArrayLike,
maturities: Sequence[int],
n_factors: int = ...,
periods_per_year: float = ...,
) -> dict[str, Any]:
ACM regression-based term premium (Adrian-Crump-Moench 2013).
The three-step estimator: PCA factors from the yield panel, a factor
VAR(1), excess-return regressions on lagged factors and contemporaneous
innovations, the convexity-adjusted lambda0/lambda1 price-of-risk OLS,
then affine log-price recursions with and without the prices of risk.
Decomposes fitted yields into risk-neutral (expected-short-rate) yields
and the term premium.
UNITS: `yields` is T x M of ANNUALIZED continuously-compounded zero-coupon
log yields in DECIMAL (divide percent by 100 — the convexity terms are
quadratic, so percent input misprices them, it does not just rescale).
`maturities` are integer PERIODS (months for monthly data), ascending,
containing 1; excess returns need n - 1 in the grid for each return
maturity n (contiguous grid or pairs; interpolate the curve first if
needed). Returns factors, factor_loadings, mu/phi/sigma, rx_maturities,
a/beta/c, sigma2, lambda0/lambda1, delta0/delta1, A/B, A_rn/B_rn,
fitted / risk_neutral / term_premium (T x M, fitted = risk_neutral +
term_premium), var/rx/short_rate/yield R-squareds, and the echoed
inputs maturities / n_factors / periods_per_year. The premium's
LEVEL is estimation-sample sensitive; compare only across models fit on
the same sample.
Returned keys: `A`, `A_rn`, `B`, `B_rn`, `a`, `beta`, `c`, `delta0`,
`delta1`, `factor_loadings`, `factors`, `fitted`, `lambda0`, `lambda1`,
`maturities`, `mu`, `n_factors`, `periods_per_year`, `phi`,
`risk_neutral`, `rx_maturities`, `rx_rsquared`, `short_rate_rsquared`,
`sigma`, `sigma2`, `term_premium`, `var_rsquared`, `yield_rsquared`.
DSGE-lite¶
dsge_solve¶
def dsge_solve(
a: _ArrayLike, b: _ArrayLike, c: _ArrayLike, n_predetermined: int
) -> dict[str, Any]:
Blanchard-Kahn solution of a linear RE model A E[y_{t+1}] = B y_t + C z.
Returns the decision rule g, the law of motion p/q, eigenvalue_moduli, and verdict.
Returned keys: `eigenvalue_moduli`, `g`, `p`, `q`, `verdict`.
quantile & growth-at-risk¶
quantile_regression¶
def quantile_regression(
y: _ArrayLike,
x: _ArrayLike,
taus: Sequence[float] | None = ...,
se: str = ...,
) -> dict[str, Any]:
Linear quantile regression (statsmodels QuantReg, all defaults).
IRLS check-loss coefficients with Powell kernel-sandwich standard errors
(Epanechnikov kernel, Hall-Sheather bandwidth; `se="robust"` is the only
flavor). Include the constant column in `x`. Returns per-tau `params`,
`bse`, `tvalues`, `iterations`, `bandwidth`, `sparsity`, plus a single
`converged` bool over all taus.
Returned keys: `bandwidth`, `bse`, `converged`, `iterations`, `params`,
`sparsity`, `taus`, `tvalues`.
quantile_lp¶
def quantile_lp(
y: _ArrayLike,
shock: _ArrayLike,
taus: Sequence[float] | None = ...,
horizons: int = ...,
n_lag_controls: int = ...,
) -> dict[str, Any]:
Quantile local projections: irf[tau][h] with Powell-sandwich se[tau][h].
Per horizon, `y_{t+h}` on `[shock_t, const, p lags of y and shock]` at
each tau (tsecon-lp design conventions); matches statsmodels QuantReg on
the identical design. `converged[tau][h]` is the per-fit IRLS flag: a
False entry hit the 1000-iteration cap before the 1e-6 coefficient
tolerance, so that point of the IRF is the last iterate, not a verified
check-loss minimum — do not quote it without refitting.
Returned keys: `converged`, `horizons`, `irf`, `se`, `taus`.
Further arguments, with defaults: `taus` (None), `horizons` (12),
`n_lag_controls` (4).
growth_at_risk¶
def growth_at_risk(
y: _ArrayLike,
conditions: _ArrayLike,
horizon: int = ...,
taus: Sequence[float] | None = ...,
rearrange: bool = ...,
) -> dict[str, Any]:
Growth-at-risk (Adrian-Boyarchenko-Giannone 2019).
Conditional quantiles of the h-ahead outcome on `[const, conditions,
y_t]`, evaluated at every t — `current` is the latest risk read. `taus`
must be strictly increasing and `horizon >= 1`. `rearrange` applies the
Chernozhukov-Fernandez-Val-Galichon monotone sort across tau; `crossing`
reports whether the raw fitted quantile paths crossed either way. `bse`
carries the Newey-West overlap correction at `hac_lags = horizon - 1`
lags; `bse_powell` is the uncorrected Powell sandwich (the statsmodels
`QuantReg` number), identical to `bse` at `horizon = 1`. `converged` is
the per-tau IRLS flag (aligned with `params`): a False entry hit the
1000-iteration cap before the 1e-6 tolerance, so that tau's
coefficients — and the fitted quantiles and `current` risk read built
from them — are the last iterate, not a verified check-loss minimum.
Returned keys: `bse`, `bse_powell`, `converged`, `crossing`, `current`,
`fitted`, `fitted_raw`, `hac_lags`, `horizon`, `params`, `taus`.
functional shocks (FVAR / FLP)¶
functional_pca¶
Functional PCA of a T x M curve panel (Inoue-Rossi 2021).
Returns mean_curve, eigenfunctions (K x M), scores (T x K), eigenvalues,
explained, total_variance. Sign: each eigenfunction's largest-|.| entry
is positive.
Further arguments, with defaults: `n_factors` (3).
Returned keys: `eigenfunctions`, `eigenvalues`, `explained`,
`mean_curve`, `scores`, `total_variance`.
flp¶
def flp(
y: _ArrayLike,
scores: _ArrayLike,
horizons: int = ...,
n_lag_controls: int = ...,
hac_maxlags: int | None = ...,
) -> dict[str, Any]:
Functional local projection: y_{t+h} on ALL K scores jointly + const + lags of y, Newey-West HAC (maxlags = h + n_lag_controls default).
Returns horizons, n_factors, betas ((H+1) x K), covs (joint (H+1) x K x K),
se, nobs. Per-element se conditions on the scores: inconsistent for
functional_pca-estimated scores (generated regressors) — flp_scenario's
w'beta contrasts are immune; see the functional-shocks model card.
Further arguments, with defaults: `hac_maxlags` (None).
Returned keys: `betas`, `covs`, `horizons`, `n_factors`, `nobs`, `se`.
flp_scenario¶
def flp_scenario(
y: _ArrayLike,
curves: _ArrayLike,
delta: _ArrayLike,
n_factors: int = ...,
horizons: int = ...,
n_lag_controls: int = ...,
hac_maxlags: int | None = ...,
) -> dict[str, Any]:
IRF of y to a whole-curve scenario delta (length M): FPCA, joint FLP, then response w'beta_h with se sqrt(w' Cov_h w).
Returns horizons, weights, response, se, betas, explained.
Further arguments, with defaults: `n_factors` (3), `n_lag_controls` (2),
`hac_maxlags` (None).
Returned keys: `betas`, `explained`, `horizons`, `response`, `se`,
`weights`.
fvar_scenario¶
def fvar_scenario(
y: _ArrayLike,
curves: _ArrayLike,
delta: _ArrayLike,
n_factors: int = ...,
lags: int = ...,
horizon: int = ...,
) -> dict[str, Any]:
FVAR scenario: VAR([scores, y], scores FIRST) with Cholesky identification; score innovation set to w = phi'delta, outcome's own structural shock zero (impact response of y is a modeling assumption).
Returns horizons, weights, response_outcome, responses ((H+1) x (K+1),
scores first then outcome), implied_outcome_innovation.
Further arguments, with defaults: `n_factors` (3), `lags` (2), `horizon`
(10).
Returned keys: `horizons`, `implied_outcome_innovation`,
`response_outcome`, `responses`, `weights`.
structural breaks¶
bai_perron¶
def bai_perron(
y: _ArrayLike, x: _ArrayLike, max_breaks: int = ..., trim: float = ...
) -> dict[str, Any]:
Bai-Perron multiple breaks: DP global partitions, sequential supF(l+1|l) selection at 5%, per-regime OLS, and Bai (1997) break-date confidence intervals; x is T x q with all coefficients switching (include your constant).
Returned keys: `break_dates`, `break_dates_by_m`, `bse`, `ci_lower_90`,
`ci_lower_95`, `ci_scale`, `ci_upper_90`, `ci_upper_95`, `h`,
`n_breaks`, `params`, `regime_ends`, `regime_ssr`, `regime_starts`,
`ssr_path`, `sup_f_crit`, `sup_f_seq`.
Further arguments, with defaults: `max_breaks` (5), `trim` (0.15).
sup_f_test¶
Andrews sup-F (Quandt) unknown-break test with Hansen (1997) approximate p-value; returns stat, p_value, break_date, and the full f_path over the trimmed dates.
Returned keys: `break_date`, `dates`, `f_path`, `h`, `p_value`, `stat`.
Further arguments, with defaults: `trim` (0.15).
smooth local projections¶
smooth_lp¶
def smooth_lp(
y: _ArrayLike,
shock: _ArrayLike,
horizons: int = ...,
n_lag_controls: int = ...,
lam: float | str | None = ...,
degree: int = ...,
n_basis: int | None = ...,
penalty_order: int = ...,
lambda_grid: Sequence[float] | None = ...,
n_folds: int = ...,
hac_maxlags: int | None = ...,
band: str | None = ...,
band_alpha: float = ...,
band_seed: int = ...,
band_n_sim: int = ...,
) -> dict[str, Any]:
Smooth local projections (Barnichon-Brownlees 2019): the IRF as a penalized B-spline in the horizon, estimated jointly across horizons.
`lam`: a float fixes the smoothing parameter (0.0 reproduces the
per-horizon `lp(se="hac")` point estimates with the default basis);
"cv"/None cross-validates it by leave-h-block-out CV over `lambda_grid`.
`lambda_grid=None` uses the default **scale-relative** grid — a 17-point
log ladder spanning eight decades, anchored to the mean diagonal of the
spline block of the stacked X'X, so the selected smoothing (and the
unit-normalized IRF) is invariant to rescaling `y` and/or `shock`; an
explicit `lambda_grid` is absolute (in the units of your data) and used
verbatim, and `cv_grid` always reports the grid actually searched.
`penalty_order=2` shrinks the IRF toward
a straight line as `lam` grows. `se` conditions on `lam` and does not
account for shrinkage bias; `irf_raw`/`se_raw` are the unsmoothed
per-horizon HAC LP for comparison. Keys: horizons, irf, se, lambda_used,
cv_grid, cv_scores, theta, irf_raw, se_raw.
**Bands.** `band=None` (default) returns no band. `"pointwise"`, `"sup-t"`,
`"sidak"` or `"bonferroni"` add `lower`/`upper` over the horizons of this
response (`K = horizons + 1`, `band_scope="horizon"`) with
`critical_value`, `pointwise_critical_value`, `n_cells`, `n_cells_used`
and `cov_se_max_rel_diff` (~machine epsilon here: the band covariance IS
the delta-method matrix behind `se`; None where no covariance is built).
A pointwise band covers one horizon at a time; the other three cover every
horizon at once at `1 - band_alpha`.
Smooth LP is the one estimator here that already had the full cross-horizon
covariance — the path is `irf_h = B_h' theta` for a single jointly-estimated
coefficient vector — so `"sup-t"` needs no extra estimation and no
compromise. It simulates `band_n_sim` Gaussian draws from `band_seed`, so
the band is a **pure function** of that seed. The usual smooth-LP caveat
still applies and is not a band problem: `se` conditions on `lam` and
ignores the penalty's shrinkage bias, so any band here is centred on a
shrunk estimator. Method: Montiel Olea and Plagborg-Møller.
Further arguments, with defaults: `n_lag_controls` (4), `degree` (3),
`n_basis` (None), `n_folds` (5), `hac_maxlags` (None).
extreme value theory¶
gpd_fit¶
def gpd_fit(
y: _ArrayLike,
threshold: float | None = ...,
quantile: float = ...,
p_tail: Sequence[float] | None = ...,
) -> dict[str, Any]:
Peaks-over-threshold GPD tail fit with McNeil-Frey (2000) VaR/ES.
Fits a generalized Pareto distribution by MLE to the strict exceedances
of `y` over `threshold` (default: the empirical `quantile` of `y`,
numpy-linear convention; both the threshold and its quantile are
reported). `xi` is the tail index — scipy's `genpareto` `c` is the same
quantity (matches `scipy.stats.genpareto.fit(z, floc=0)`, polished, at
1e-6). Standard errors are observed-information; when `xi <= -0.5`
(Smith 1985 irregularity) they are reported but `se_valid` is False.
`var`/`es` are the McNeil-Frey POT tail quantiles at each `p_tail` entry
(default [0.99, 0.995, 0.999]; each must reach beyond the threshold:
`1 - p < n_exceed / n`) in the units of `y` — fit losses (`-returns` or
`abs(returns)`) to read them as risk numbers; `es` is NaN where
`xi >= 1`. At least 10 exceedances are required. Keys: threshold,
threshold_quantile, n, n_exceed, exceed_rate, xi, beta, se_xi, se_beta,
se_valid, loglik, converged, p_tail, var, es.
gev_fit¶
def gev_fit(
y: _ArrayLike,
block_size: int | None = ...,
return_periods: Sequence[float] | None = ...,
) -> dict[str, Any]:
GEV block-maxima fit with return levels.
With `block_size=None`, `y` IS the pre-computed block maxima; otherwise
`y` is cut into consecutive non-overlapping blocks of that length (a
trailing partial block is dropped) and each block contributes its
maximum. Fits GEV(`xi`, `mu`, `sigma`) by MLE — `xi` is the tail index;
scipy's `genextreme` shape is `c = -xi` (matches
`scipy.stats.genextreme.fit(maxima)`, polished, at 1e-6). Standard
errors are observed-information with the same `se_valid` certification
as `gpd_fit` (`xi <= -0.5` reported, not certified). `return_levels`
are the `1 - 1/T` GEV quantiles at each `return_periods` entry (default
[10, 50, 100] blocks; each `T > 1`). At least 10 maxima are required.
Keys: xi, mu, sigma, se_xi, se_mu, se_sigma, se_valid, loglik,
converged, n_maxima, block_size, return_periods, return_levels.
static copulas¶
pseudo_obs¶
Pseudo-observations: the average-rank probability-scale transform.
`u[i, j] = rank of x[i, j] within column j / (n + 1)`, ties assigned
their average rank — exactly scipy `rankdata(method="average")/(n+1)`
(golden-pinned, ties included). The `n + 1` denominator keeps every
value strictly inside (0, 1), which the copula quantile transforms
require. Ranks see only order, so any strictly INCREASING transform of
a margin (logs, standardization, exp) leaves the output — and any
copula fitted to it — bit-identical (property-tested). A strictly
decreasing transform instead reverses that margin's ranks (`u -> 1 - u`
when there are no ties), flipping the sign of the fitted dependence —
the standard copula invariance is increasing-only. This is the one-line
companion to `copula_fit`: `copula_fit(pseudo_obs(x))`. Accepts any
number of columns (the transform is columnwise); `copula_fit` itself
is bivariate in this slice.
copula_fit¶
Fits a bivariate copula to (n, 2) probability-scale pseudo-observations.
`u` must lie strictly inside (0, 1): rank/PIT-transform the raw margins
first — `pseudo_obs(x)` does it in one line, and the whole workflow is
then invariant to strictly increasing transforms of each margin (the
point of the copula decomposition; property-tested — a decreasing
transform flips the sign of the dependence instead). At least 20 pairs
required.
`family`: "gaussian" (param `rho`), "t" (`rho`, `nu`), "clayton"
(`theta` > 0, lower-tail), "gumbel" (`theta` >= 1, upper-tail), "frank"
(`theta`, either sign). Clayton/Gumbel model positive dependence only
in this slice (rotations deferred) and raise a teaching error when the
empirical Kendall tau is <= 0. `method`: "mle" (maximum likelihood,
observed-information SEs — matches a polished scipy optimum of the
statsmodels log-density at 1e-6) or "tau" (Kendall-tau inversion, the
statsmodels `fit_corr_param` route — for "t", tau pins `rho` and `nu`
is profiled by MLE; SEs are NaN with `se_valid` False, honestly, since
the moment-based SE is deferred).
Returns the named dependence parameter(s) (`rho` / `rho` + `nu` /
`theta`, also stacked in `params` with `param_names`), their SEs
(`se_rho` / `se_nu` / `se_theta`, stacked in `se`, certified by
`se_valid`), `loglik`, `aic`, `bic`, the empirical Kendall `tau` and
the fit-implied `tau_implied`, and the closed-form tail-dependence
coefficients `tail_lower`/`tail_upper` (Gaussian/Frank 0 — the classic
reason a Gaussian fit understates joint crashes; t symmetric
Demarta-McNeil; Clayton lower 2^(-1/theta); Gumbel upper
2 - 2^(1/theta)). Keys: family, method, n, params, param_names, rho,
nu, theta, se, se_rho, se_nu, se_theta, se_valid, loglik, aic, bic,
tau, tau_implied, tail_lower, tail_upper, converged (rho/nu/theta and
their se_* appear per family).
copula_select¶
def copula_select(
u: _ArrayLike,
families: Sequence[str] | None = ...,
method: str = ...,
) -> dict[str, Any]:
Fits several copula families to the same (n, 2) pseudo-observations and ranks them by AIC/BIC, with a teaching verdict.
`families`: list of names (default all five: gaussian, t, clayton,
gumbel, frank); `method` as in `copula_fit`. Families whose domain
excludes the data (Clayton/Gumbel under Kendall tau <= 0) are
*skipped with a reason* rather than failing the call, so the default
menu works on any data. Each entry of `fits` is a full `copula_fit`
dict; `ranking_aic`/`ranking_bic` list family names best-first;
`best_aic`/`best_bic` name the winners; `verdict` states who wins, by
how much, whether AIC and BIC agree (they differ exactly when the
extra parameter is not earning its keep by BIC), what the winner
implies for tail dependence, and what was skipped and why. Keys:
fits, skipped, best_aic, best_bic, ranking_aic, ranking_bic, verdict.
kernel methods¶
kernel_ridge¶
def kernel_ridge(
x: _ArrayLike,
y: _ArrayLike,
alpha: float = ...,
kernel: str = ...,
gamma: float | None = ...,
degree: float = ...,
coef0: float = ...,
x_test: _ArrayLike | None = ...,
rff_features: int | None = ...,
seed: int = ...,
) -> dict[str, Any]:
Kernel ridge regression: exact dual solve, or the Rahimi-Recht random-Fourier-feature approximation of the rbf kernel.
Minimizes `sum_i (y_i - f(x_i))^2 + alpha * ||f||_H^2` over the RKHS of
the kernel — scikit-learn's `KernelRidge` objective (no `1/n`, no
intercept: center `y` if the kernel does not model a level). The exact
solution is `(K + alpha I) a = y` by Cholesky. `x` is `(n, k)` (or 1-D
for one regressor). Kernels in scikit-learn's exact parameterization:
`kernel="rbf"` `exp(-gamma ||x-y||^2)`, `"laplacian"`
`exp(-gamma ||x-y||_1)`, `"polynomial"` `(gamma <x,y> + coef0)^degree`,
`"linear"` `<x,y>`. `gamma=None` resolves to `1 / n_features`
(scikit-learn's default); the linear kernel has no gamma and refuses
one. `degree`/`coef0` act on the polynomial kernel only and are refused
at non-default values elsewhere (nothing is silently ignored).
`rff_features=D` switches to the random-Fourier-feature primal
approximation (Rahimi & Recht 2007): `z(x) = sqrt(2/D) cos(Wx + b)`
with `W ~ N(0, 2 gamma I)`, `b ~ U[0, 2 pi)` drawn from a Philox stream
keyed by `seed` (same seed, bit-identical features), then ridge on `z`
— `O(n D^2)` instead of `O(n^3)`, converging to the exact fit as `D`
grows. rbf only; `seed` is refused in exact mode. `x_test` (`(m, k)`)
adds `predicted`. `alpha=0` is the interpolating fit and is refused
when `K` is not positive definite (scikit-learn silently falls back to
least squares there; tsecon raises and names `alpha`).
Keys: `dual_coef` (exact mode; the `a` of `f(x) = sum_i a_i k(x, x_i)`,
scikit-learn's `dual_coef_`) or `coef` (RFF mode; the `D` primal
weights), `fitted`, `predicted` (only when `x_test` is given),
`kernel`, `gamma` (resolved; None for linear), `n_rff_features` (None
in exact mode).
Validated against scikit-learn 1.9.0 `KernelRidge` — `dual_coef_`,
`predict(X)` and `predict(X_test)` for all four kernels at 1e-8
(independent package). The RFF approximation is a Monte-Carlo object
and is property-tested (seeded determinism; error against the exact
fit falling with `D`), not golden-pinned.
kernel_regression¶
def kernel_regression(
x: _ArrayLike,
y: _ArrayLike,
bandwidth: float | Sequence[float] | _ArrayLike | None = ...,
kind: str = ...,
kernel: str = ...,
bandwidth_method: str = ...,
block: int | None = ...,
x_test: _ArrayLike | None = ...,
) -> dict[str, Any]:
Nadaraya-Watson or local-linear kernel regression of y on x
((n, k), k <= 3, or 1-D for one regressor) with a product Gaussian
kernel, at a fixed or cross-validated bandwidth.
Conventions are statsmodels `KernelReg(reg_type="lc" | "ll",
var_type="c"*k)` exactly: `kind="nadaraya_watson"` is the local
constant `sum_i K_h(x_i - x) y_i / sum_i K_h(x_i - x)`;
`kind="local_linear"` (default — no boundary bias) is the intercept of
the kernel-weighted least squares of `y` on `[1, x_i - x]`, solved
through the pseudoinverse as statsmodels does. `kernel`: `"gaussian"`
only (the one statsmodels validates against; compact-support kernels
are deferred). The bandwidth is the kernel's standard deviation per
column, in the column's units.
`bandwidth_method="fixed"` (default) uses `bandwidth` (a positive
scalar broadcast to every column, or one value per column) as given.
`"loo_cv"` minimizes the leave-one-out least-squares criterion
`n^-1 sum_i (y_i - g_{-i}(x_i))^2` (statsmodels `cv_loo`).
`"block_cv"` minimizes the leave-block-out criterion (Chu & Marron
1991): predicting `y_i` drops the `2*block + 1` observations with
`|j - i| <= block` (default `block = ceil(n^(1/3))`), so serially
correlated neighbours never vote on their own errors — leave-one-out
undersmooths badly under autocorrelated errors, and this is the method
to use for time-series regressors. Selection is a 21-point log grid
on a common multiple of the Scott reference `1.06 sd(x_j) n^(-1/(4+k))`
over `[0.05, 20]`, golden-section refinement, then per-column
coordinate refinement for `k >= 2`; deterministic, and not statsmodels'
Nelder-Mead path (the criterion value at any bandwidth matches
statsmodels at 1e-10; the search reaches a criterion no worse than
fmin's). Under the CV methods `bandwidth` must be omitted and under
`"fixed"`/`"loo_cv"` `block` must be omitted — a conflicting argument
raises rather than being ignored.
Keys: `fitted` (at the training rows), `predicted` (only when `x_test`
is given; NaN where every training weight underflows), `bandwidth`
(resolved, one per column), `bandwidth_method`, `block` (resolved
half-width under `"block_cv"`, else None), `cv_criterion` (the
leave-one-out criterion under `"fixed"`/`"loo_cv"`, the leave-block-out
criterion under `"block_cv"`, at the reported bandwidth), `effective_df`
(`tr(S)` of the linear smoother: from `k+1` (local linear) or `1`
(Nadaraya-Watson) at huge bandwidths up to `n` at tiny ones), `kind`,
`kernel`, `bandwidth_at_boundary` (True when a selected bandwidth sits
on a wall of the search range — the criterion was still falling, so
the reported value is the search's limit, not an interior optimum;
typically a target with no detectable signal), and
`n_criterion_evaluations` (0 under `"fixed"`).
Validated against statsmodels 0.15.0 `KernelReg.fit()` at fixed
bandwidths (`k = 1, 2`, both estimators) at 1e-8 and `cv_loo` at 1e-10
(independent package); the leave-block-out criterion and
`effective_df` are documented-formula transcriptions (no package
computes them) pinned at 1e-10.
structured penalties and post-selection¶
group_lasso¶
def group_lasso(
x: _ArrayLike,
y: _ArrayLike,
groups: npt.NDArray[np.integer] | Sequence[int],
alpha: float,
l1_ratio: float = ...,
group_weights: str | _ArrayLike | None = ...,
tol: float = ...,
max_iter: int = ...,
) -> dict[str, Any]:
Group LASSO (Yuan-Lin 2006) / sparse-group LASSO (Simon et al. 2013) by block coordinate descent with exact per-block Lipschitz constants.
Objective: (1/(2n))||y - Xb||^2 + alpha*[(1 - l1_ratio)*sum_g w_g
||b_g||_2 + l1_ratio*||b||_1] — the crate's `lasso` scaling, so
`l1_ratio=1` IS `lasso(x, y, alpha)`. No intercept, no standardization
inside: center `y`, standardize `x`. `groups` is one integer label per
column (any integers, contiguous or not; integer arrays are passed
through untouched). `group_weights`: "sqrt_size" (default, w_g =
sqrt(|g|)), "none" (w_g = 1), or one positive weight per distinct label
in ascending label order. `tol` is the dimensionless coefficient-change
rule shared with `lasso` and also bounds the KKT residual relative to
max_j |x_j'y|/n; `max_iter` caps the block sweeps.
Returns `coef`, `n_iter`, `converged` (True only when the sweep rule AND
the KKT certificate are met; when False the last iterate is returned —
read `kkt_violation`), `active_groups` (labels with a nonzero block),
`active_set` (nonzero column indices), `objective`, `kkt_violation`
(largest subgradient KKT residual at `coef` — a readable optimality
certificate for this convex problem), `max_rel_change`, and `alpha_max`
(smallest alpha with the all-zero solution).
Validation: independent KKT certificate <= 1e-8 on every fixture case
(~2e-13 achieved) plus cross-package agreement with skglm at 1e-8
(~1.5e-12 achieved); reductions to `lasso` at 1e-8.
post_lasso¶
def post_lasso(
x: _ArrayLike,
y: _ArrayLike,
alpha: float,
l1_ratio: float = ...,
tol: float = ...,
max_iter: int = ...,
) -> dict[str, Any]:
Post-LASSO OLS refit (Belloni-Chernozhukov 2013): LASSO / elastic net
with elastic_net's objective, then OLS on the selected columns to
remove the shrinkage bias. No intercept, no standardization inside.
Returns `support` (selected indices), `coef_lasso` (first stage),
`coef_ols` (the refit, zeros off-support; minimum-norm least squares on
the support), `n_selected`, `rss`.
NO standard errors, deliberately: OLS standard errors after a
data-driven selection are invalid (the selection event depends on the
same sample). For inference on a target coefficient use `pds_lasso`.
Validation: refit pinned to scikit-learn LinearRegression
(fit_intercept=False) on the scikit-learn Lasso/ElasticNet support at
1e-10 (~8e-15 achieved).
Further arguments, with defaults: `l1_ratio` (1.0), `tol` (1e-08),
`max_iter` (100000).
pds_lasso¶
def pds_lasso(
y: _ArrayLike,
d: _ArrayLike,
x: _ArrayLike,
alpha: float | str | None = ...,
hac_lags: int | None = ...,
tol: float = ...,
max_iter: int = ...,
) -> dict[str, Any]:
Post-double-selection LASSO (Belloni-Chernozhukov-Hansen 2014) for
the coefficient on a treatment d with high-dimensional controls x,
with Newey-West (Bartlett) HAC inference from the shared HAC engine.
LASSO `y` on `x` and `d` on `x`, take the union of supports, OLS `y` on
[d, x_union]; the treatment is never penalized. Center `y` and `d`,
standardize `x` (no intercept, nothing standardized inside). `alpha`:
a float applied to both LASSOs, or "bic" (default) — the per-equation
BIC pick along `lasso_path`'s default grid. `hac_lags`: None (default)
= the Newey-West rule floor(4 (n/100)^(2/9)); a positive integer = that
Bartlett lag truncation; 0 = classical spherical-errors standard errors.
HAC covariance carries n/(n-k) (statsmodels HAC with
use_correction=True); `p_value`/`conf_int` use the standard normal in
both modes (statsmodels use_t=False).
Returns `coef`, `se`, `t_stat`, `p_value`, `conf_int` (95% (lo, hi)),
`support_y`, `support_d`, `union_support`, `n_controls_selected`,
`alpha_y`, `alpha_d`, `hac_lags_resolved`.
Validation: Monte-Carlo grade for coverage (R hdm / Stata pdslasso not
runnable here) — the seeded design in structured_properties.rs measures
the PDS interval's coverage against the single-selection interval's
undercoverage, numbers on the model card; exact leg against statsmodels
HAC / nonrobust OLS on the selected union at 1e-8 (~1e-14 achieved).
Further arguments, with defaults: `tol` (1e-08), `max_iter` (100000).
Trees and forests¶
regression_tree¶
def regression_tree(
x: _ArrayLike,
y: _ArrayLike,
max_depth: int | None = ...,
min_samples_leaf: int = ...,
min_samples_split: int = ...,
x_test: _ArrayLike | None = ...,
) -> dict[str, Any]:
CART regression tree (Breiman et al. 1984) with scikit-learn's best-split conventions.
Reproduces scikit-learn 1.9.0 `DecisionTreeRegressor(criterion=
"squared_error", splitter="best", max_features=None)` — an
independent-package golden (fixtures/trees.json): test predictions
at 1e-12, `n_leaves`/`depth` exact, `feature_importances_` at 1e-10,
the sorted (feature, threshold) multiset at 1e-12. Exact matching is
possible because the fixture proves every stored case tie-free (the
same tree under five sklearn `random_state` values): sklearn breaks an
exact tie between two features by its private RNG's visit order, this
tree by the lowest feature index, and only two-row nodes make such
ties likely. sklearn works in float32, this tree in float64.
Conventions: squared-error criterion; threshold = midpoint of the two
adjacent sorted distinct values (values within 1e-7 count as one); a
split leaves at least `min_samples_leaf` rows on both sides; a node
with fewer than `min_samples_split` rows, at `max_depth` (None =
unbounded), or pure is a leaf; leaves predict the training mean.
Returns `fitted` (n, leaf mean per training row), `predicted` (rows of
`x_test`, or None), `n_nodes`, `n_leaves`, `depth` (root = 0),
`feature_importance` (impurity-based, normalized to one; zeros if the
tree never split), and `splits` (list of [feature, threshold] pairs
sorted by (feature, threshold)). Keys: fitted, predicted, n_nodes,
n_leaves, depth, feature_importance, splits.
Raises ValueError for NaN/inf (naming the array), an `x_test` column
mismatch, `min_samples_leaf < 1`, `min_samples_split < 2`, and
`insufficient data: {got} observations, at least {needed} required`
when n < max(min_samples_split, 2 * min_samples_leaf).
random_forest¶
def random_forest(
x: _ArrayLike,
y: _ArrayLike,
n_trees: int = ...,
max_features: str | int | None = ...,
max_depth: int | None = ...,
min_samples_leaf: int = ...,
bootstrap: str = ...,
block_length: int | None = ...,
seed: int = ...,
x_test: _ArrayLike | None = ...,
quantiles: Sequence[float] | None = ...,
importance: str = ...,
importance_groups: Sequence[int] | None = ...,
permutation_block: int | None = ...,
n_permutations: int | None = ...,
) -> dict[str, Any]:
Random forest for regression (Breiman 2001) with time-series-aware resampling, out-of-bag error, quantile regression forests (Meinshausen 2006), and grouped block-permutation importance.
Each tree is the CART tree of `regression_tree` grown on a row resample
(drawn rows act as multiplicity weights; rows never drawn are the
tree's out-of-bag rows), visiting `max_features` random columns per
node; the forest averages the trees. Validation grade (honest): the
deterministic tree is golden-pinned to scikit-learn 1.9.0 and
`random_forest(bootstrap="none", max_features="all", n_trees=1,
min_samples_leaf=1)` reproduces `regression_tree` bit-for-bit, which is
how the forest inherits that golden; the full forest's randomness is
tsecon's own Philox stream (one SeedSequence substream per tree, so it
is bit-identical at any thread count; same `seed` same forest,
different `seed` different forest) and is validated by seeded
Monte-Carlo property tests whose measured numbers the model card
quotes (Friedman #1 out-of-sample R^2, autocorrelation preserved by
block resampling, out-of-bag optimism, quantile-band coverage,
importance recovery).
`n_trees` (default 500); `max_features` (None = "third") in {"sqrt",
"third" (max(1, p // 3)), "all", or an int in 1..=p}; `max_depth` (None =
unbounded); `min_samples_leaf` (default 5); `bootstrap` in {"iid"
(default, Efron), "block" (Künsch moving block), "stationary"
(Politis-Romano, geometric blocks of mean `block_length`), "none"
(every tree sees every row; no out-of-bag rows)} — `block_length` is
required for "block"/"stationary" and refused for "iid"/"none"; `seed`
(default 0); `x_test` (m, p) rows to predict; `quantiles` (strictly
inside (0, 1), strictly increasing; requires `x_test`) turns on the
quantile regression forest; `importance` in {"none" (default),
"impurity", "block_permutation"}; `importance_groups` (one integer
label per column — give all lags of one variable one label so they
are permuted and credited as one unit; needs `importance` != "none";
a label vector, not data — pass a list of ints); `permutation_block`
(rows per permuted block; None = ceil(n ** (1/3)); 1 = single-row) and
`n_permutations` (None = 10) act only under
importance="block_permutation" and are refused elsewhere.
Returns `fitted` (n, in-sample forest prediction), `predicted` (m or
None), `oob_prediction` (n; NaN where a row was never out-of-bag; None
under bootstrap="none"), `oob_mse` (over rows with an out-of-bag
prediction; None under bootstrap="none"), `importance` (per unit;
impurity sums to one, block_permutation is the mean out-of-bag MSE
increase in units of y^2, may be negative; None under "none"),
`importance_groups_resolved` (the unit label each `importance` entry
refers to), `quantile_predictions` ((m, len(quantiles)), never
crossing; None without `quantiles`), `n_trees`, `max_features_resolved`.
Keys: fitted, predicted, oob_prediction, oob_mse, importance,
importance_groups_resolved, quantile_predictions, n_trees,
max_features_resolved.
Gotchas, measured and quoted on the model card. (1) OUT-OF-BAG ERROR IS
OPTIMISTIC ON TIME SERIES: an out-of-bag row's temporal neighbours are
in-bag in the trees that score it and, with persistent predictors and
autocorrelated errors, carry its error — the property suite measures
OOB/POOS MSE ratios of about 0.70 under AR(0.9) errors vs about 0.84 under iid
errors on the same persistent design; report pseudo-out-of-sample
metrics and prefer bootstrap="block"/"stationary". (2) A PERSISTENT
IRRELEVANT PREDICTOR'S IMPORTANCE IS INFLATED when the relevant
predictors are persistent too (the forest uses it as a time proxy);
grouping the lags of a variable keeps permuted rows dynamically
possible and stops dilution across collinear lags, but block
permutation does NOT remove that inflation — for a row-wise forest
scored row-wise, single-row and block permutation give the same mean
importance; compare against a control instead. (3) Impurity
importance favours columns with many distinct values.
Raises ValueError for NaN/inf (naming the array), `insufficient data:
{got} observations, at least {needed} required` when n < 2 *
min_samples_leaf, unknown string options (listing the accepted
values), malformed `quantiles` (naming the fix), `importance_groups`
of the wrong length (naming both lengths), a block length outside
1..=n, and every inert-kwarg combination above.
Trend filtering and boosting¶
l1_trend_filter¶
def l1_trend_filter(
y: _ArrayLike,
lam: float,
order: int = ...,
penalty: str = ...,
tol: float | None = ...,
max_iter: int | None = ...,
) -> dict[str, Any]:
L1 trend filtering (Kim, Koh & Boyd 2009) — a piecewise-linear trend
with data-chosen knots — or, with penalty="l2", the Hodrick-Prescott
filter on the same objective.
Minimizes over the trend `x`, with `D` the `order`-th difference
operator: `penalty="l1"`: `(1/2)||y - x||^2 + lam * ||D x||_1` (most
`order`-th differences exactly zero — `order=2` a piecewise-linear
trend whose kinks are the `knots`, `order=1` piecewise-constant, the
fused LASSO on the level); `penalty="l2"`: `(1/2)||y - x||^2 +
(lam/2) * ||D x||^2`, which for `order=2` is exactly `hp_filter(y,
lam)` (same minimizer, same `lam`; 1600 quarterly), solved in closed
form. Scan an L1 `lam` downward from `lam_max`, the value at which the
trend collapses to the least-squares polynomial of degree `order - 1`.
Solver: Kim-Koh-Boyd primal-dual interior point on the banded dual
(O(n) per step, no n×n matrices) plus an exact active-set polish.
`tol` is the relative duality gap at which it stops (`duality_gap <=
tol * objective`), `max_iter` the Newton-step budget; both act only
under `penalty="l1"`, and passing either explicitly under `"l2"` — a
closed-form solve with nothing to iterate — raises rather than being
ignored (`None`, the default, means 1e-8 / 10000 where they apply).
Returns `trend`, `cycle` (`y - trend`), `knots` (indices into the
`order`-th differences where `|(D trend)_i|` exceeds `max(1e-6 *
max|D y|, 1e-12 * max|y|)`; under `"l2"` nearly every index),
`n_knots`, `duality_gap` (the certificate — an upper bound on
`objective - optimum`), `objective`, `converged` (`duality_gap <= tol *
objective`; always True on the closed-form paths), `n_iter` (0 on
closed-form paths), and `lam_max`. Keys: trend, cycle, knots, n_knots,
duality_gap, objective, converged, n_iter, lam_max.
Validation: an independent KKT / duality-gap certificate re-derived
in the tests for every fixture case (relative gap <= 1e-8 asserted),
cvxpy + Clarabel third-party trends at 1e-8, the `lam -> 0` and
`lam >= lam_max` limits, and the `hp_filter` identity at 1e-10.
boosting¶
def boosting(
x: _ArrayLike,
y: _ArrayLike,
learning_rate: float = ...,
n_steps: int = ...,
stop: str = ...,
x_test: _ArrayLike | None = ...,
) -> dict[str, Any]:
Componentwise L2 boosting with single-column least-squares base
learners (Buhlmann & Yu 2003; Buhlmann 2006 — the R mboost glmboost
engine): a slow-learning variable selector read as sequential ARDL
building.
From `F_0 = 0` (no intercept — pass a centered `y` and centered,
typically standardized, columns), each step regresses the current
residual on every column separately, picks the column with the
smallest residual sum of squares (ties to the smallest index), and
adds `learning_rate` times that fit. Seedless and deterministic.
`learning_rate` in (0, 1] (0.1 conventional; 1.0 unshrunk greedy);
`n_steps >= 1`. The boosting operator `B_m = B_{m-1} + nu H_j (I -
B_{m-1})` is tracked exactly in a rank-m factored form — no n×n matrix
— and its trace is the degrees of freedom in Buhlmann's (2006)
corrected AIC, `log(RSS_m/n) + (1 + df_m/n)/(1 - (df_m+2)/n)` (`+inf`
where `df_m + 2 >= n`). `stop="aic"` reports the AIC-minimizing step,
`stop="none"` the last; the paths are returned either way.
Returns `coef` (length p, at the reported step), `coef_path` (n_steps
× p; row m is the model after m + 1 iterations), `selected` (column
chosen at each step), `rss_path`, `df_path`, `aic_path`, `best_step`
(0-based index into the path arrays), `fitted` (`x @ coef`), and
`predicted` (`x_test @ coef`, or None). Keys: coef, coef_path,
selected, rss_path, df_path, aic_path, best_step, fitted, predicted.
Validation (graded honestly): a transcription of the published
algorithm into dense NumPy — the operator formed explicitly, so the
trace is exact by construction — pins `coef_path`, `selected`,
`df_path`, `aic_path` at 1e-12; R mboost is not runnable in the build
environment, so this is not a third-party run. Properties: RSS
nonincreasing, the small-step limit is OLS on the selected support,
AIC recovers a sparse truth's support.
Neural (MLP, echo state network)¶
mlp_regression¶
def mlp_regression(
x: _ArrayLike,
y: _ArrayLike,
hidden: Sequence[int] | int | None = ...,
activation: str = ...,
alpha: float = ...,
solver: str = ...,
learning_rate: float | None = ...,
batch_size: int | None = ...,
max_epochs: int = ...,
validation_fraction: float = ...,
patience: int | None = ...,
n_seeds: int = ...,
seed: int = ...,
standardize: bool = ...,
x_test: _ArrayLike | None = ...,
) -> dict[str, Any]:
Feed-forward neural regressor (one or two hidden layers) with a seed ensemble, early stopping on a TEMPORAL validation split, and scikit-learn MLPRegressor's exact objective — the "NN" of the macro forecasting horse races, and tsecon's only native neural net (no framework dependency; torch / foundation-model adapters are out of core by scope ruling).
`hidden`: tuple or list of layer widths, default `(16,)`; an int is
one layer; at most two layers by design. `activation`: "tanh"
(default), "relu", "logistic". `alpha` (1e-4): L2 penalty on the
weights, sklearn scale — the objective is
`(1/(2n)) sum (y - f(x))^2 + (alpha/(2n)) sum_l ||W_l||_F^2`,
intercepts unpenalized. `solver`: "adam" (default; sklearn's
constants; `learning_rate` None -> 1e-3; `batch_size` None -> one
full-batch step per epoch, an int -> seeded shuffled mini-batches;
`max_epochs` 500; early stopping with `patience` None -> 20 epochs
without a relative-1e-4 improvement of the validation loss, best
epoch's weights kept) or "lbfgs" (tsecon's L-BFGS on the full
objective with the analytic gradient, `max_epochs` capping its
iterations; passing `learning_rate`, `batch_size`, or `patience`
explicitly under lbfgs RAISES — they cannot apply — so leave them
None). `validation_fraction` (0.2; 0 disables early stopping; at most
0.5): the LAST `floor(validation_fraction * n)` rows are held out —
never a random split. `standardize` (True): the scaler for `x` and
`y` is fit on the TRAINING rows only and replayed on the validation
rows and on `x_test`. `n_seeds` (5) members from independent Philox
substreams of `seed` (0) are averaged. `x_test`: optional
`(n_test, p)` rows to predict.
Returns `fitted` (ensemble mean on every row of `x`, original y
scale), `predicted` / `member_predictions` (ensemble mean and the
`(n_seeds, n_test)` array on `x_test`; None without it),
`train_loss_path` / `validation_loss_path` (lists of per-member
per-epoch arrays; two entries — initial and final — under lbfgs; the
validation path is empty when validation_fraction=0), `best_epoch`
and `converged` (per member; True = early stopping fired / L-BFGS
converged, False = ran out of max_epochs), `n_parameters`, `weights`
(per member `{"coefs": [...], "intercepts": [...]}` in sklearn's
fan_in x fan_out layout, standardized scale), `n_train`,
`n_validation`, `x_mean`, `x_scale`, `y_mean`, `y_scale` (the
training-row scaler; identity when standardize=False), `solver`,
`activation`. Keys: fitted, predicted, member_predictions,
train_loss_path, validation_loss_path, best_epoch, converged,
n_parameters, weights, n_train, n_validation, x_mean, x_scale,
y_mean, y_scale, solver, activation.
Validated against scikit-learn 1.9.0 MLPRegressor (independent
package, fixtures/neural.json): forward pass = sklearn predict at its
fitted weights (1e-12), objective (1e-10), analytic gradient =
sklearn's own backprop (1e-10) and a central finite difference (1e-6
relative), gradient norm at sklearn's converged weights (1e-8). The
optimizer trajectory is deliberately not pinned. Estimator grade:
property / Monte Carlo — recovers y_t = sin(2 y_{t-1}) + e_t out of
sample (R^2 0.94 mini-batch Adam, 0.95 lbfgs, 0.80 all-defaults; 0.75
linear), the ensemble beats the mean member in every replication
(Jensen) and the median member in a majority (7/10 Rust draws, 10/10
NumPy draws) on a documented overfitting DGP, early stopping fires on
an easy problem and cannot at max_epochs=1. Reproducibility: single-threaded Rust, every draw a
pure function of `seed` — bit-identical on the same build; across
platforms only the last ulp of libm tanh/exp can differ, so the
cross-platform promise is statistical (seed-ensemble)
reproducibility. Errors name the array with NaN/inf (x, y, x_test),
list the accepted activation/solver names, name the two-layer limit,
and report `insufficient data: {got} observations, at least {needed}
required` with the validation split counted.
echo_state_network¶
def echo_state_network(
x: _ArrayLike,
y: _ArrayLike,
reservoir_size: int = ...,
spectral_radius: float = ...,
leak_rate: float = ...,
input_scaling: float = ...,
sparsity: float = ...,
washout: int = ...,
ridge_alpha: float = ...,
seed: int = ...,
x_test: _ArrayLike | None = ...,
) -> dict[str, Any]:
Echo state network (reservoir computing; Jaeger 2001; Lukosevicius 2012): a fixed sparse random recurrent reservoir, a leaky-integrator tanh state recursion, and a ridge-trained linear readout.
`reservoir_size` (200) units; `spectral_radius` (0.9) the reservoir
matrix is rescaled to (leading-eigenvalue modulus from a dense
eigenvalue decomposition; values above 1 accepted); `leak_rate` a in
(0, 1] (1.0 = plain ESN): `s_t = (1 - a) s_{t-1} + a tanh(W s_{t-1} +
W_in u_t)`, `s_0 = 0`, no reservoir bias; `input_scaling` (1.0):
W_in uniform on [-input_scaling, input_scaling]; `sparsity` (0.1):
the CONNECTIVITY, i.e. the fraction of nonzero reservoir entries
(standard normal values); `washout` (50): leading rows discarded
before the readout fit; `ridge_alpha` (1e-6): readout penalty,
minimizing `||y - Z b||^2 + ridge_alpha ||b||^2` on
Z_t = [1, u_t, s_t] (scikit-learn Ridge(fit_intercept=False) scale;
the constant column is penalized like every coefficient, Lukosevicius
eq. 9); `seed` (0); `x_test`: optional `(n_test, p)` rows treated as
the CONTINUATION of `x` (states carry on from the last training
state; no washout re-applied).
Returns `fitted` (readout on the rows that entered the fit, length
n - washout), `predicted` (on x_test, else None), `readout`
(coefficients on [1, u, s], length 1 + p + reservoir_size),
`spectral_radius_achieved` (recomputed on the scaled matrix),
`reservoir_size`, `n_washout`, `n_train` (n - washout). Keys: fitted,
predicted, readout, spectral_radius_achieved, reservoir_size,
n_washout, n_train.
Validation (fixtures/neural.json): the state recursion on an explicit
small reservoir is pinned at 1e-12 against a NumPy transcription that
reservoirpy 0.4.2's Reservoir (same explicit W / Win / lr) reproduced
exactly at generation time — a third-party pin of the mechanics; the
readout at 1e-10 against the closed-form ridge, cross-checked there
against scikit-learn Ridge; the spectral radius against
numpy.linalg.eigvals (1e-6). Estimator grade: property — NARMA-10
out-of-sample NRMSE 0.32 (mean over four data seeds) with
input_scaling=0.3 and otherwise default settings on 1000 training
rows, 0.19 with reservoir_size=400 on 2000 rows (the all-defaults
call averages 0.43: input_scaling=1 over-drives tanh for NARMA's u in
[0, 0.5]); achieved radius within 1e-6 of the target; same seed
bit-identical, different seeds differ. Reproducibility:
single-threaded, every draw a pure function of `seed`; last-ulp
libm/eigenvalue differences across platforms. Errors name the array
with NaN/inf (x, y, x_test); `washout >= n` names the fix; fewer than
two rows after the washout reports `insufficient data: {got}
observations, at least {needed} required` with the washout counted.
Threshold confidence sets (Hansen 1997/2000)¶
setar_threshold_ci¶
def setar_threshold_ci(
y: _ArrayLike,
p: int,
delay: int = ...,
trim: float = ...,
delays: Sequence[int] | None = ...,
constant: bool = ...,
level: float = ...,
het_robust: bool = ...,
slope_level: float | None = ...,
slope_region_level: float | None = ...,
null_threshold: float | None = ...,
) -> dict[str, Any]:
Hansen (1997/2000) likelihood-ratio confidence set for the threshold
of a two-regime SETAR(p), built on exactly the setar fit.
The fit is `setar(y, p, delay/delays, trim, constant)` itself — the
returned `threshold`, `delay`, `thresholds` and `ssr_path` are
bit-identical to its (`delays` overrides `delay`, as there). Over the
candidate grid the profile `LR_n(gamma) = nobs * (S(gamma) - S_min) /
S_min` (Hansen 2000; Hansen 1997 for the TAR) is inverted against the
closed-form critical value `c = -2 ln(1 - sqrt(level))` — the `level`
quantile of `P(xi <= x) = (1 - exp(-x/2))^2`, Hansen (2000) Table 1:
4.50 / 5.94 / 7.35 / 10.59 at 80 / 90 / 95 / 99% — giving the set
`{gamma : LR_n(gamma) <= eta2 * c}`. The set always contains the
estimate (`LR_n = 0` there), is typically asymmetric, and CAN BE
DISJOINT, so it is returned as a list of closed `[low, high]`
`intervals` (maximal runs of in-set candidates, grid endpoints as in
Hansen's own programs) with `is_connected`, `n_intervals`, and the
convex hull `ci_low` / `ci_high`; `in_set` flags each candidate and
`n_in_set` counts them. `LR_n` is a step function, constant between
adjacent candidates, so a `null_threshold` gamma_0 (any value inside
`[thresholds[0], thresholds[-1]]`) is evaluated at the largest
candidate `<= gamma_0` (`null_threshold_used`): `lr_at_null` and the
p-value `pvalue_at_threshold = p(lr_at_null / eta2)`, `p(x) = 1 - (1 -
exp(-x/2))^2` — test inversion at one point; all three are None when
no null is passed.
`het_robust=True` applies Hansen's (2000, section 3.4)
heteroskedasticity correction: the critical value is scaled by
`eta2 = E[e^2 (x'delta)^2 | q = gamma] / (sigma^2 E[(x'delta)^2 | q =
gamma])`, estimated as his programs do — at the threshold estimate,
regress `(x'delta_hat)^2` and `e_hat^2 (x'delta_hat)^2` each on a
quadratic polynomial in the threshold variable `y_{t-d}` (with
intercept), take the ratio of the fitted values at the estimate, and
divide by `sigma_hat^2 = S_min / nobs`. `eta2` is exactly 1 otherwise.
With no threshold effect the fitted ratio can be non-positive; that is
refused with a teaching error rather than returned as a negative
scale.
`slope_level=0.95` (say) adds Hansen's (2000, section 3.3)
CONSERVATIVE slope intervals: the union, over every candidate in the
`slope_region_level` threshold set (default 0.80, his applied
convention), of the conventional per-regime intervals `b_j(gamma) +/-
z se_j(gamma)` — classical per-regime SEs as `setar` reports (at the
estimate exactly `bse_low` / `bse_high`), or White HC0 SEs under
`het_robust`. They come back as `slope_ci_low` and `slope_ci_high`,
each a 2 x k nested list `[[low-regime coefficients], [high-regime
coefficients]]` (constant first, then lags 1..p), with the region used
in `slope_region_low` / `slope_region_high` / `slope_n_region`; all
are None without `slope_level`, and `slope_region_level` passed
without `slope_level` RAISES (it would be inert).
Validation: closed-form critical values and p-values pinned at 1e-14
(Table 1 reproduced to the printed decimals); the LR profile, the
`eta2` convention, the intervals and the slope unions pinned at 1e-10
against an independent NumPy transcription (fixtures/setar_ci.json —
no third-party threshold-CI code runs in the fixture container);
coverage MEASURED by seeded Monte Carlo in the crate's property tests
and quoted in the model card (asymptotically conservative for a fixed
threshold effect, as Hansen's theory says).
Further arguments, with defaults: `delay` (1), `trim` (0.15), `delays`
(None), `constant` (True), `level` (0.95), `het_robust` (False),
`slope_level` (None), `slope_region_level` (None: 0.80 when slope
intervals are requested), `null_threshold` (None).
Returned keys: `threshold`, `delay`, `nobs`, `k`, `level`, `lr_crit`,
`lr_crit_scaled`, `eta2`, `het_robust`, `thresholds`, `ssr_path`,
`lr_stat`, `in_set`, `intervals`, `n_intervals`, `is_connected`,
`ci_low`, `ci_high`, `n_in_set`, `null_threshold_used`, `lr_at_null`,
`pvalue_at_threshold`, `slope_level`, `slope_region_level`,
`slope_region_low`, `slope_region_high`, `slope_n_region`,
`slope_ci_low`, `slope_ci_high`.
Distributed-lag panel regressions (climate-impact specification)¶
panel_distributed_lag¶
def panel_distributed_lag(
outcome: _ArrayLike,
regressors: _ArrayLike,
lags: int,
powers: int = ...,
entity_effects: bool = ...,
time_effects: bool = ...,
entity_trends: bool = ...,
se_type: str = ...,
bandwidth: float | None = ...,
eval_points: _ArrayLike | None = ...,
mask: _ArrayLike | None = ...,
) -> dict[str, Any]:
Distributed-lag panel regression — the climate-impact specification of Dell-Jones-Olken (2012) and Burke-Hsiang-Miguel (2015):
`mask` (default None = a balanced panel) is an N x T array of 0/1 flags,
1 where the entity is observed in that period, for an UNBALANCED panel:
cells outside the mask are ignored and may hold NaN; a lagged row enters
only when the entity is observed in every period t - L ..= t, and the
default `eval_points` (the pooled regressor mean) runs over the observed
cells (validated against linearmodels PanelOLS on the Arellano-Bond
EmplUK panel and a seeded ragged panel, fixtures/panel_unbalanced.json).
COST: on an UNBALANCED panel with `time_effects=True` the time effects
are partialled out through the projected time dummies (one per observed
period) by a rank-revealing least-squares step, which is CUBIC in the
number of periods — measured at N = 6: 0.11 s at T = 400, 0.69 s at
T = 800, 5.2 s at T = 1600, 38 s at T = 3200, against 3 ms for the same
panel with no mask or with `time_effects=False` (both linear, and a
mask of all ones takes the balanced path bit-identically). Long
unbalanced panels are practical only without time effects, or by
trimming T.
y_it = sum_{l=0..L} beta_l x_{i,t-l} [+ sum_l gamma_l x^2_{i,t-l}]
+ alpha_i + delta_t [+ g_i t] + e_it
`outcome` is N x T; `regressors` is k x N x T (weather variables:
strictly exogenous, no lagged outcome — a lagged dependent variable
would put Nickell bias back into the within estimator; use `panel_lp`
with a bias correction for dynamic panels). `lags` is L: lags 0..L of
every regressor enter and the first L periods of each entity are
dropped (on an unbalanced panel, see `mask` above, a lagged row needs
every one of its lags observed). `powers=1` is the linear response, `powers=2` adds the lags
of the square (the BHM quadratic response). `entity_effects` (True),
`time_effects` (True) and `entity_trends` (False; requires entity
effects) choose the fixed effects; at least one effect is required.
`se_type` is "nonrobust", "cluster" (by entity — the DJO default) or
"driscoll_kraay" (the BHM robustness choice; needs a long T).
`bandwidth` is the Driscoll-Kraay lag truncation and acts ONLY under
`se_type="driscoll_kraay"` (4.0 when omitted there); passing it
explicitly with any other `se_type` raises instead of being silently
absorbed. `eval_points` (1-D, default None) are the points at which the
marginal effect of the cumulative quadratic response is evaluated for
every regressor; it acts ONLY under `powers=2` (None there means each
regressor's pooled sample mean) and passing it under `powers=1` raises
— the linear cumulative response has one constant marginal effect,
`cumulative_effect` itself.
Design columns are ordered regressor-major, then power, then lag
(`names` lists them, e.g. `x0_L0`, `x0_L1`, `x0^2_L0`, ...). Returned
keys: `params`, `names`, `bse`, `tvalues`, `cov` (K x K nested lists),
`lag_effects` and `lag_se` (`[regressor][power-1][lag]`),
`cumulative_effect`, `cumulative_se` (delta method, `sqrt(1' V 1)`),
`cumulative_ci_low`, `cumulative_ci_high` (normal 95%,
`[regressor][power-1]`), and under `powers=2` `eval_points`,
`marginal_effect`, `marginal_se` (`[regressor][point]`; the marginal
effect `B_1 + 2 B_2 x` of the cumulative response), `turning_point`,
`turning_point_se` (`[regressor]`; `-B_1/(2 B_2)`, NaN when `B_2` is
exactly zero) — all five are None under `powers=1`; plus `nobs`
(`N (T - L)`), `n_entities`, `n_periods_used` (`T - L`), `lags`,
`powers`, `df_resid`, `se_type`, `entity_effects`, `time_effects`,
`entity_trends`.
Validation (fixtures/panel_dl.json): nine cases x three covariance
estimators pinned at 1e-10 against linearmodels PanelOLS (slopes, SEs,
t-statistics, full covariance; the trends variant via explicit entity x
trend regressors), the delta-method quantities against the documented
NumPy transcription at 1e-10; the cumulative-effect interval's coverage
is measured in seeded Monte Carlo and quoted on the panel model card.
The `lags=0`, `time_effects=False` call is bit-identical to `panel_fe`.
JSZ affine term structure¶
jsz_fit¶
def jsz_fit(
yields: _ArrayLike,
maturities: Sequence[int],
n_factors: int = ...,
periods_per_year: float = ...,
w: _ArrayLike | None = ...,
n_starts: int = ...,
seed: int | None = ...,
) -> dict[str, Any]:
JSZ canonical Gaussian affine term-structure model (Joslin-Singleton-Zhu 2011) by maximum likelihood.
Risk-neutral dynamics in the JSZ canonical form — ordered eigenvalues
`lambda_q` (per period), one drift `k_inf_q`, short rate `r = iota'X` —
rotated onto `n_factors` yield portfolios `P_t = w y_t` (default `w`:
the first principal-component loadings of the panel) that are priced
WITHOUT error, the remaining `M - n_factors` yield directions with iid
error `sigma_e`. The P-measure VAR(1) of the portfolios is concentrated
out by OLS (exactly statsmodels `VAR(1)`), `k_inf_q` and `sigma_e` are
profiled analytically, and the numerical search runs over `lambda_q`
and the Cholesky factor of `sigma` only, from JSZ's recommended start
(the OLS eigenvalues) plus `n_starts - 1` seeded perturbations, the best
basin polished by BFGS then Nelder-Mead. The fit depends on `w` only
through its row space (any basis of the same portfolio space gives the
same `llf`, `fitted`, `lambda_q`, `k_inf_q`, `sigma_e`).
UNITS: `yields` is `T x M` ANNUALIZED continuously-compounded zero-coupon
log yields in DECIMAL (0.05, not 5.0); `maturities` are strictly
ascending integer PERIODS (months for monthly data; 1 need not be
present); `periods_per_year` (12.0) converts to the per-period quantities
the recursions price. `n_factors` (3); `w` (None: PCA loadings, an
`n_factors x M` array otherwise); `n_starts` (5); `seed` (None, meaning seed 0 — not fresh entropy; the returned
`seed` key is the value used) seeds the
perturbed starts through `tsecon_rng` and may only be passed when
`n_starts > 1` (with a single start it would be inert, so it raises).
Returns `lambda_q`, `k_inf_q` (per period), `sigma` (MLE innovation
covariance of the portfolio VAR, annualized units), `sigma_e`
(annualized), the OLS P-measure VAR `mu_p`, `phi_p` with statsmodels-
convention standard errors `mu_p_se`, `phi_p_se` and `sigma_ols`
(`sigma_u_mle`), the risk-neutral VAR in the portfolio rotation `k0_q_p`,
`k1_q_p`, the ACM-unit market prices of risk `lambda0 = mu_p - k0_q_p`,
`lambda1 = phi_p - k1_q_p`, loadings in the portfolio rotation `a_p`,
`b_p` (`fitted = a_p + b_p P`) and for the literal canonical latent state
`a_x`, `b_x`, the `fitted` yields (`T x M`), `risk_neutral` yields (the
recursion under the P-measure VAR, the acm_term_premium convention),
`term_premium = fitted - risk_neutral`, `rmse` per maturity, the
portfolios `factors` (`T x n_factors`) and weights `w`, `llf` (the log-
likelihood of the yield panel, basis-invariant), `converged`, `n_iter`,
and the echoed `maturities`, `n_factors`, `periods_per_year`, `n_starts`,
`seed`. No standard errors are reported for the Q parameters: a
numerical Hessian of a profile likelihood near a unit root is not an
honest asymptotic covariance; the P-VAR standard errors show where the
market-price-of-risk imprecision actually lives.
Validation (fixtures/jsz.json): the Riccati recursions against a
documented-formula NumPy transcription at 1e-12; the portfolio VAR(1)
against statsmodels at 1e-9; the likelihood against the documented
formula at 1e-7 absolute (~1e-11 relative); the MLE against a SciPy
multi-start optimum (cross-optimizer, lambda_q within 1e-5) on a
simulated canonical model — recovering lambda_q to 9e-5, k_inf_q to
1.3%, sigma_e to 0.3% at T = 500 — and on the real 1990-2007 GSW panel;
the AFNS special case lambda_q = (1, e^-lam, e^-lam) spans the Nelson-
Siegel loadings exactly and its convexity intercept converges at first
order in the period length to the CDR (2011) closed form of
afns_adjustment (gap 9.8e-6 at monthly, 6.1e-7 at 1/192 year).
Returned keys: `a_p`, `a_x`, `b_p`, `b_x`, `converged`, `factors`,
`fitted`, `k0_q_p`, `k1_q_p`, `k_inf_q`, `lambda0`, `lambda1`,
`lambda_q`, `llf`, `maturities`, `mu_p`, `mu_p_se`, `n_factors`,
`n_iter`, `n_starts`, `periods_per_year`, `phi_p`, `phi_p_se`,
`risk_neutral`, `rmse`, `seed`, `sigma`, `sigma_e`, `sigma_ols`,
`term_premium`, `w`.
jsz_loadings¶
def jsz_loadings(
lambda_q: _ArrayLike,
k_inf_q: float,
sigma_x: _ArrayLike,
maturities: Sequence[int],
periods_per_year: float = ...,
) -> dict[str, Any]:
The JSZ canonical bond-loading recursions at given parameters.
Evaluates `A_{n+1} = A_n + K0' B_n + 1/2 B_n' sigma_x B_n`, `B_{n+1} =
K1' B_n - iota` from `A_0 = B_0 = 0` for the literal canonical form
`K0 = (k_inf_q, 0, ..)'`, `K1 = J(lambda_q)` (diagonal; a Jordan block
with a 1 on the superdiagonal wherever two consecutive eigenvalues are
exactly equal — the AFNS pattern `(1, rho, rho)`), and returns the
per-maturity yield coefficients `a_x = -A_n/n * periods_per_year` and
`b_x = -B_n'/n` so that `y(n) = a_x + b_x X`. `lambda_q` must be ordered
non-increasing (per period); `sigma_x` is the per-period `N x N` state
innovation covariance in the canonical basis; `periods_per_year` (1.0)
only rescales the intercept.
Returned keys: `a_x`, `b_x`, `k0_q`, `k1_q`, `maturities`.
Generalized impulse responses (Koop-Pesaran-Potter)¶
var_girf¶
def var_girf(
data: _ArrayLike,
p: int,
shock_var: int = ...,
size: float = ...,
shock: str = ...,
horizon: int = ...,
n_draws: int = ...,
seed: int = ...,
trend: str = ...,
antithetic: bool = ...,
histories: int | None = ...,
bands: tuple[float, float] | None = ...,
) -> dict[str, Any]:
Generalized impulse responses (Koop-Pesaran-Potter 1996) of a linear VAR(p) by simulation — the engine's exact reduction to the closed-form impulse responses, so nonlinear GIRFs can be checked against a linear benchmark on the same footing.
For every lag window of `data` (or a seeded subsample of `histories`
of them) and every future-innovation draw, the fitted VAR is simulated
forward twice with the same innovations — once with the shock added to
the impact-period innovation, once without — and the paired difference
is averaged (common random numbers; `antithetic` (+z, -z) pairs when
set, needing an even `n_draws`). In a linear model the paired
difference is Psi_h delta for every draw and history, so the result
carries no Monte Carlo noise and does not depend on `n_draws` or
`seed` beyond rounding; `mc_se`/`draw_sd` are zero to rounding (below
1e-15), and `mc_se` is NaN below two effective draws — at the default
`n_draws=2` with `antithetic=True` there is exactly one, so `mc_se` is
NaN there (pass `n_draws=4` or `antithetic=False` for a finite value).
`shock`: "orthogonal" — `size` standard deviations of the `shock_var`-th
Cholesky-orthogonalized innovation in the variable ordering, so
`girf[h]` equals `var_irf(orth=True)[h][:, shock_var] * size` to 1e-12;
"generalized" — the Pesaran-Shin (1998) shock (innovation `shock_var`
moved by `size` standard deviations, the others by their conditional
expectation; no ordering), so `girf[h]` equals
Phi_h Sigma e_j / sqrt(sigma_jj) * size with Sigma the df-adjusted
residual covariance. `size` may be negative.
Keys: `girf` ([h][variable], h = 0..horizon, mean over histories),
`lower`/`upper` (the `bands` quantiles across histories — zero width
here), `per_history` ([history][h][variable]), `mc_se` (Monte Carlo
standard error of `girf`), `draw_sd` (across-draw spread of one
realized paired difference), `draw_lower`/`draw_upper` (mean over
histories of the per-history across-draw `bands` quantiles),
`n_histories`, `n_draws`, `n_effective_draws` (n_draws / 2 under
antithetic), `horizon`, `shock` (echo), `shock_var`, `shock_size_used`
(impact-period innovation to `shock_var` in raw units: size * P[j, j]
orthogonal, size * sqrt(sigma_jj) generalized), `shock_vector` (the
full raw innovation added at impact).
Validation: statsmodels VARResults.irf(orth=True) and the Pesaran-Shin
closed form, both at 1e-12 with a single draw (fixtures/girf.json).
Further arguments, with defaults: `shock_var` (0), `size` (1.0),
`shock` ("orthogonal"), `horizon` (10), `n_draws` (2), `seed` (0),
`trend` ("c"), `antithetic` (True), `histories` (None = every lag
window; an int draws a seeded subsample of that many — a count at or
above the number of available windows uses all of them, reported in
`n_histories`), `bands` (None = (0.16, 0.84)).
Memory: the engine refuses up front, as a `ValueError` naming `n_draws`,
`horizon` and the number of histories, any call whose draw buffers and
per-history results would exceed its fixed 2 GiB budget; below the budget
the buffers are allocated fallibly, so an allocator refusal is the same
error, never an abort.
threshold_var_girf¶
def threshold_var_girf(
data: _ArrayLike,
p: int,
threshold_index: int = ...,
delay: int = ...,
trim: float = ...,
delays: Sequence[int] | None = ...,
constant: bool = ...,
shock_var: int = ...,
size: float = ...,
shock: str = ...,
horizon: int = ...,
n_draws: int = ...,
seed: int = ...,
regime: str = ...,
histories: int | None = ...,
bands: tuple[float, float] | None = ...,
antithetic: bool = ...,
) -> dict[str, Any]:
Regime-dependent generalized impulse responses (Koop-Pesaran-Potter
1996) of the two-regime threshold VAR: fits threshold_var with the
same p/threshold_index/delay|delays/trim/constant, then
simulates the fitted nonlinear system forward from the sample's actual
lag windows, regime-switching period by period.
Histories are every lag window t >= max(p, delay) of `data` (the shock
hits period t; its regime is decided by data[t - delay, threshold_index]
<= threshold), in time order; `regime`="low"/"high" keeps only the
windows whose shock-date regime is that one, and `histories`=m a
seeded subsample of m of the selected windows. For each history and
each of `n_draws` future-innovation draws the model is simulated twice
with the same standard-normal draws — the shocked path adds the shock
at impact — and the paired difference is averaged. At every period
each path reads its own regime from its own simulated window and
scales the common draw by the Cholesky factor of THAT regime's ML
residual covariance (sigma_low/sigma_high), so a path that crosses the
threshold switches both coefficients and innovation covariance (Balke
2000's regime-by-regime draws; R tsDyn's GIRF pools residuals because
its TVAR fits one covariance). The impact shock is scaled by the
covariance of the regime the history is in at the shock date:
"orthogonal" — `size` standard deviations of the `shock_var`-th
Cholesky-orthogonalized innovation of that regime; "generalized" — the
Pesaran-Shin shock of that regime; the raw impact therefore differs
across regimes when their covariances do (see `shock_size_used`).
`antithetic` uses (+z, -z) pairs (even `n_draws`); the Monte Carlo
standard error is computed from the pair means.
Keys: `girf` ([h][variable], mean over the used histories),
`lower`/`upper` (the `bands` quantiles across the used histories — the
KPP history-conditional distribution), `per_history`
([history][h][variable]), `mc_se` (Monte Carlo standard error of
`girf`, histories fixed; NaN below two effective draws), `draw_sd`
(across-draw spread of one realized paired difference),
`draw_lower`/`draw_upper` (mean over histories of the per-history
across-draw `bands` quantiles), `girf_low_regime`/`girf_high_regime`
(means over the used histories of each regime; None when the selection
holds none of that regime), `history_regimes` (0 low / 1 high per used
history), `history_times` (the shock date t of each), `n_histories`,
`n_low_histories`, `n_high_histories`, `n_draws`, `n_effective_draws`,
`horizon`, `shock` (echo), `shock_var`, `shock_size_used` ([low, high]
impact-period innovation to `shock_var` in raw units), `shock_vector`
([regime][variable] raw innovation added at impact), `regime` (echo),
`threshold`, `delay`, `threshold_index`.
Reproducible: one Philox substream per (history, draw) spawned from
`seed`, bit-identical at any thread count and across processes.
Validation (honest grade): the engine's linear reduction is pinned at
1e-12 against statsmodels and the Pesaran-Shin closed form
(`var_girf`); the regime-switching simulation is pinned at 1e-10
against an independent NumPy transcription of the documented engine
reproducing its random streams (fixtures/girf.json); sign asymmetry,
size non-proportionality, regime dependence, 1/sqrt(n_draws)
convergence and antithetic variance reduction are measured by seeded
Monte Carlo property tests (see the model card). No third-party TVAR
GIRF runs in the build container.
Further arguments, with defaults: `threshold_index` (0), `delay` (1),
`trim` (0.1), `delays` (None; a list searches the delay and overrides
`delay`), `constant` (True), `shock_var` (0), `size` (1.0), `shock`
("orthogonal"), `horizon` (20), `n_draws` (500), `seed` (0), `regime`
("all"), `histories` (None = every selected window; an int draws a
seeded subsample of that many — a count at or above the number of
selected windows uses all of them, reported in `n_histories`), `bands`
(None = (0.16, 0.84)), `antithetic` (True).
Memory: the engine refuses up front, as a `ValueError` naming `n_draws`,
`horizon` and the number of histories, any call whose draw buffers and
per-history results would exceed its fixed 2 GiB budget; below the budget
the buffers are allocated fallibly, so an allocator refusal is the same
error, never an abort.
Cointegrating regressions (FM-OLS / DOLS / CCR)¶
fmols¶
def fmols(
y: _ArrayLike,
x: _ArrayLike,
trend: str = ...,
kernel: str = ...,
bandwidth: float | None = ...,
bandwidth_rule: str | None = ...,
force_int: bool = ...,
df_adjust: bool = ...,
diff: bool | None = ...,
x_trend: str | None = ...,
) -> dict[str, Any]:
Phillips-Hansen (1990) fully modified OLS (FM-OLS) of one cointegrating vector, with asymptotically valid (mixed-normal) inference.
`y` is the regressand (length T), `x` the (T, k) matrix of I(1)
regressors (do NOT add your own constant: deterministics come from
`trend`). The static OLS `y = x'beta + d'delta + e` is super-consistent
but its t-statistics are invalid — serial correlation in `e` and
correlation between `e` and the regressor innovations `dx` leave a
second-order bias and a nuisance-parameter limit. FM-OLS corrects the
regressand for endogeneity (`y+ = y - omega_12 Omega_22^-1 eta_2`) and
subtracts the serial-correlation bias `lambda+_12` from the moment
equations, using the kernel long-run covariance of the residual system
`eta = (OLS residual, detrended dx)`; the corrected estimator has
covariance `omega_1.2 (Z'Z)^-1` and standard-normal t-statistics.
`trend`: "n", "c" (default), "ct", "ctt" (constant, trend, quadratic
trend; the trend runs 1..T). `kernel`: "bartlett" (default), "parzen",
"quadratic-spectral". `bandwidth`: an explicit kernel bandwidth (>= 0;
Bartlett/Parzen weight lag j by k(j/(bandwidth+1)) for j <=
floor(bandwidth), quadratic spectral by k(j/bandwidth)), or None
(default) to select it by `bandwidth_rule`: "newey-west" (default;
arch's rule — the Newey-West 1994 plug-in on the unit-weighted sum of
the residual system with ceil(4 (T/100)^rate) pilot lags) or
"andrews" (the Andrews 1991 AR(1) parametric plug-in on the same
series). Passing `bandwidth_rule` together with an explicit `bandwidth`
RAISES (the rule would be inert). `force_int` (default True, as arch)
ceils the bandwidth — automatic or explicit; the automatic one is also
capped at T - 1. `df_adjust` (default False) scales the covariance by
(T-1)/(T-1-p), with T-1 the rows of the residual system and p the
number of ESTIMATED COEFFICIENTS — the k regressors AND the
deterministics of `trend`, i.e. `len(params)` (at the default
`trend="c"` and k = 2 that is 199/196, not 199/197). `x_trend`
(default None = `trend`; must carry at least the terms of `trend`)
sets the deterministics the regressors are detrended with before
differencing; `diff` (default None, which behaves as False) removes
the trend from the differences instead of the levels and RAISES when
the effective x_trend has no trend term (it would be inert) — and note
that `diff=False` passed EXPLICITLY raises there too, for the same
reason: the default is the `None` sentinel, not `False`.
Keys: `estimator`, `params` (x columns first, then the deterministics —
see `param_names`), `se`, `tvalues`, `pvalues` (two-sided normal),
`cov`, `param_names`, `resid` (length `nobs` = T, the full sample),
`nobs`, `n_x`, `n_det`, `trend`, `x_trend`, `kernel`, `bandwidth` (the
one actually used), `bandwidth_rule` (None when explicit), `force_int`,
`diff`, `df_adjust`, `long_run_variance` (`omega_1.2`, df-scaled),
`omega` / `lambda` / `sigma` (the (1+k)x(1+k) long-run, one-sided
long-run and short-run covariances of the residual system, nested
lists), `n_lags` (positive lags the window covered), `rsquared`,
`rsquared_adj`, `ols_params` and `ols_se` (the plain static OLS for
comparison — its SEs are NOT valid for inference).
Validation: arch 8.0 `FullyModifiedOLS` at 1e-10 across every trend,
kernel and option (fixtures/fmols.json); t-statistic size and
super-consistency measured by seeded Monte Carlo (model card).
Further arguments, with defaults: `trend` ("c"), `kernel` ("bartlett"),
`bandwidth` (None), `bandwidth_rule` (None: "newey-west"), `force_int`
(True), `df_adjust` (False), `diff` (None: False), `x_trend` (None:
`trend`).
ccr¶
def ccr(
y: _ArrayLike,
x: _ArrayLike,
trend: str = ...,
kernel: str = ...,
bandwidth: float | None = ...,
bandwidth_rule: str | None = ...,
force_int: bool = ...,
df_adjust: bool = ...,
diff: bool | None = ...,
x_trend: str | None = ...,
) -> dict[str, Any]:
Park (1992) canonical cointegrating regression (CCR) of one cointegrating vector, with asymptotically valid inference.
Same inputs, options and keys as `fmols`. Where FM-OLS corrects the
regressand and the moment equations, CCR transforms the DATA: with
`Sigma`, `Lambda`, `Omega` the short-run, one-sided and two-sided
long-run covariances of the residual system `eta` and `beta_ols` the
static OLS coefficients, `x* = x - (Sigma^-1 Lambda_2)'eta` and `y* =
y - (Sigma^-1 Lambda_2 beta_ols + kappa)'eta` with `kappa = (0,
Omega_22^-1 omega_21)`, and `params` is the OLS of `y*` on `[x*, d]`
over t = 2..T, with covariance `omega_1.2 (Z*'Z*)^-1`. Asymptotically
equivalent to FM-OLS; the two differ in finite samples.
`df_adjust` scales the covariance by (T-1)/(T-1-p) with p the
estimated coefficients (regressors and deterministics), as documented —
arch 8.0's `CanonicalCointegratingReg.fit` scales only `omega_11`
(an operator-precedence slip); everything else is arch-exact.
Keys: `estimator`, `params`, `se`, `tvalues`, `pvalues`, `cov`,
`param_names`, `resid`, `nobs`, `n_x`, `n_det`, `trend`, `x_trend`,
`kernel`, `bandwidth`, `bandwidth_rule`, `force_int`, `diff`,
`df_adjust`, `long_run_variance`, `omega`, `lambda`, `sigma`, `n_lags`,
`rsquared`, `rsquared_adj`, `ols_params`, `ols_se`.
Validation: arch 8.0 `CanonicalCointegratingReg` at 1e-10
(fixtures/fmols.json), the `df_adjust` scaling as documented.
Further arguments, with defaults: `trend` ("c"), `kernel` ("bartlett"),
`bandwidth` (None), `bandwidth_rule` (None: "newey-west"), `force_int`
(True), `df_adjust` (False), `diff` (None: False), `x_trend` (None:
`trend`).
dols¶
def dols(
y: _ArrayLike,
x: _ArrayLike,
trend: str = ...,
lags: int | None = ...,
leads: int | None = ...,
ic: str | None = ...,
common: bool | None = ...,
max_lag: int | None = ...,
max_lead: int | None = ...,
cov_type: str = ...,
kernel: str = ...,
bandwidth: float | None = ...,
bandwidth_rule: str | None = ...,
force_int: bool = ...,
df_adjust: bool = ...,
) -> dict[str, Any]:
Stock-Watson (1993) / Saikkonen (1991) dynamic OLS (DOLS) of one
cointegrating vector: the static regression augmented with lags lags
and leads leads of the regressor differences (the contemporaneous
difference is always included), so the augmented error is orthogonal
to the regressor innovations and OLS on the augmented design is
asymptotically mixed normal.
`y` (length T) and `x` (T, k) as for `fmols`; `trend` as there. The
regression runs over the T - 1 - lags - leads rows every term is
defined on, design `[x, deterministics, dx_{t-lags}, ..., dx_t, ...,
dx_{t+leads}]` (k columns per block; the trend runs 1..nobs over that
sample). `lags` / `leads` (default None) fix the counts; when either is
None it is chosen by minimising `ic` — "bic" (default), "aic" or
"hqic": `ln(RSS/nobs) + n_params c/nobs` — over 0..`max_lag` /
0..`max_lead` (default None = ceil(12 (T/100)^(1/4))) on the COMMON
sample of the largest candidate (ties to the smaller lag, then lead),
the chosen model then refit on its own sample; `common` (default None
= False) restricts the search to lags == leads. The sentinel rule:
`ic` passed with both `lags` and `leads` fixed RAISES, `max_lag` passed
with `lags` fixed RAISES, `max_lead` with `leads` RAISES, and `common`
with both fixed RAISES (each would be inert). A search whose largest
candidate has no residual degrees of freedom is refused (arch runs it
underdetermined).
`cov_type`: "unadjusted" (default) — `sigma2_HAC (Z'Z/n)^-1 / n` with
`sigma2_HAC` the kernel long-run variance of the residuals; "robust" —
the kernel-HAC sandwich `(Z'Z/n)^-1 S_HAC (Z'Z/n)^-1 / n` on the scores.
`kernel`, `bandwidth`, `bandwidth_rule` and `force_int` (default False,
arch's DOLS default) as for `fmols`, the automatic bandwidth chosen on
the residuals ("unadjusted") or the scores ("robust"); `df_adjust`
(default False) scales the covariance by nobs/(nobs - n_params).
Keys: `params` (the cointegrating vector: x columns then the
deterministics — `param_names`), `se`, `tvalues`, `pvalues`, `cov`,
`param_names`, `full_params` / `full_se` / `full_cov` /
`full_param_names` (every coefficient incl. the difference blocks),
`resid` (length `nobs`), `nobs` (the augmented regression's rows),
`n_total` (T), `n_x`, `n_det`, `n_params`, `trend`, `lags`, `leads`,
`selected` (False when both were fixed), `ic`, `ic_value` (the
minimised criterion; NaN when both were fixed), `max_lag` / `max_lead`
(the caps the search used), `common`, `cov_type`, `kernel`,
`bandwidth`, `bandwidth_rule`, `force_int`, `df_adjust`,
`long_run_variance` (kernel LRV of the residuals at `bandwidth`,
df-scaled — the `sigma2_HAC` of the unadjusted covariance),
`rsquared`, `rsquared_adj`, `ols_params`, `ols_se` (the plain static
OLS on the full sample — SEs NOT valid for inference).
Validation: arch 8.0 `DynamicOLS` at 1e-10 across every trend, both
covariance types, the three criteria, fixed/searched/common/capped
leads and lags (fixtures/fmols.json).
Further arguments, with defaults: `trend` ("c"), `lags` (None), `leads`
(None), `ic` (None: "bic"), `common` (None: False), `max_lag` (None),
`max_lead` (None), `cov_type` ("unadjusted"), `kernel` ("bartlett"),
`bandwidth` (None), `bandwidth_rule` (None: "newey-west"), `force_int`
(False), `df_adjust` (False).
Conditional VAR forecasts, residual diagnostics, lag-order selection¶
var_conditional_forecast¶
def var_conditional_forecast(
data: _ArrayLike,
conditions: Sequence[Sequence[float | None]] | _ArrayLike,
lags: int = ...,
trend: str = ...,
steps: int | None = ...,
alpha: float = ...,
) -> dict[str, Any]:
Conditional (hard-path) VAR forecast: the forecast of every series when some cells of the future path are pinned to given values.
`conditions` is a nested list with one row per horizon and one entry per
series: a number pins that (horizon, series) cell, None (or NaN) leaves
it free — a NumPy array with NaN for the free cells or a pandas DataFrame
with missing entries works too. Rows beyond `len(conditions)` up to
`steps` are free; `steps` defaults to `len(conditions)`. At least one
cell must be pinned (the all-free case is `var_forecast`).
COST: the guard on `steps` is a MEMORY budget (`steps * k` cells times
the constrained cells must stay inside 2^24 doubles), but the work is
QUADRATIC in `steps` — measured on a k = 3 VAR(2): 0.03 s at 1 000
steps, 0.46 s at 4 000, 7.6 s at 16 000, 31 s at 32 000 (x4 per
doubling), so a `steps` the budget admits can still run for hours.
Forecast horizons are tens of periods in practice; treat five figures
as a typo.
Method: Doan-Litterman-Sims (1984) / Waggoner-Zha (1999) — Gaussian
conditioning of the joint forecast-error distribution on the pinned
cells, equivalently the unconditional path plus the response to the
minimum-norm future shocks that deliver the conditions; identical to a
Kalman smoother with the free future cells set missing
(Bańbura-Giannone-Lenza 2015), the second golden leg.
Keys (each path `steps x k`, row h = horizon h + 1): `point` (pinned
cells hold their condition exactly), `unconditional` (bitwise
`var_forecast(...)["point"]`), `cov` ([h][i][j] conditional covariance
per horizon; pinned cells have a zero row/column), `se` (exactly 0 at
pinned cells), `unconditional_se`, `lower`/`upper` (innovation
uncertainty only, coefficients treated as known), `shocks` (implied
reduced-form innovations), `orth_shocks` (Cholesky-orthogonalised in the
column order — the only ordering-dependent key), `constrained`,
`n_constrained`, `mahalanobis` (squared Sigma-norm of the implied
shocks), `mahalanobis_pvalue` (chi2(n_constrained) tail: small means the
model is pushed hard), `steps`, `alpha`.
Further arguments, with defaults: `lags` (2), `trend` ("c"), `steps`
(None: `len(conditions)`), `alpha` (0.05).
var_diagnostics¶
def var_diagnostics(
data: _ArrayLike,
lags: int = ...,
trend: str = ...,
nlags: int = ...,
) -> dict[str, Any]:
Residual diagnostics of a fitted VAR(p): multivariate Portmanteau
(unadjusted and small-sample adjusted, chi2(k^2 (nlags - lags));
nlags must exceed lags), multivariate Jarque-Bera with its skewness
and kurtosis components (Cholesky orthogonalisation in the column order,
the statsmodels test_normality convention; Doornik-Hansen is not
provided), and the stability roots.
Keys: `portmanteau`, `portmanteau_adjusted`, `portmanteau_df`,
`portmanteau_pvalue`, `portmanteau_adjusted_pvalue`, `nlags`,
`jarque_bera`, `jarque_bera_pvalue`, `jarque_bera_df`, `skewness`,
`skewness_pvalue`, `kurtosis`, `kurtosis_pvalue`, `skewness_components`,
`kurtosis_components`, `roots` (reciprocal-root moduli, descending;
stable iff the last exceeds 1), `eigenvalue_moduli` (companion
eigenvalue moduli, descending; stable iff the first is below 1),
`is_stable`, `nobs`, `k`, `lags`.
Matches statsmodels `test_whiteness(adjusted=False/True)`,
`test_normality`, `roots`, `is_stable` at 1e-10.
Further arguments, with defaults: `lags` (2), `trend` ("c"), `nlags` (10).
var_select_order¶
VAR lag-order selection by AIC/BIC/HQIC/FPE on a common sample
(statsmodels VAR.select_order): every candidate p is fitted after
dropping the first max_lags - p rows so the criteria are comparable.
Candidates start at p = 0 with trend="c" (intercept-only baseline) and
at p = 1 with trend="n"; ties go to the smaller order.
Keys: `aic`, `bic`, `hqic`, `fpe` (selected orders), `candidates`,
`aic_values`, `bic_values`, `hqic_values`, `fpe_values`, `max_lags`,
`trend`.
Further arguments, with defaults: `max_lags` (8), `trend` ("c").
Multiple forecast comparisons (Reality Check / SPA, MCS, StepM)¶
spa_test¶
def spa_test(
benchmark_losses: _ArrayLike,
model_losses: _ArrayLike,
block_size: int | None = ...,
reps: int = ...,
bootstrap: str = ...,
studentize: bool = ...,
nested: bool = ...,
seed: int = ...,
) -> dict[str, Any]:
White's (2000) Reality Check and Hansen's (2005) test for Superior
Predictive Ability: does the BEST of m competing models beat the
benchmark once the search over all of them is accounted for?
`benchmark_losses` is the benchmark's loss series over the `n`
evaluation periods and `model_losses` a `T x m` array with one loss
column per competing model (a 1-D array is one model), index-aligned —
e.g. squared errors from `backtest` runs under the same scheme. With
`d_{t,k} = benchmark_t - model_{k,t}` (positive favours the model) the
null is `max_k E[d_k] <= 0`, and the statistic is `sqrt(n) max_k dbar_k /
omega_k` (`studentize=True`, Hansen's SPA) or `sqrt(n) max_k dbar_k`
(`studentize=False`, White's RC), where `omega_k^2` is the
stationary-bootstrap long-run variance of Hansen (2005, eq. 9) with
restart probability `1/block_size` (or, with `nested=True`, the bootstrap
variance of the resampled mean over the same resamples). The null
distribution is a block bootstrap of the whole loss-differential panel
(rows resampled together), re-centred three ways: `p_value_upper`
re-centres every model (White's original Reality Check, conservative when
poor models pad the comparison), `p_value_consistent` leaves models
significantly worse than the benchmark — by Hansen's `sqrt(2 log log n)`
threshold — un-centred (the recommended p-value, also returned as
`p_value`), and `p_value_lower` re-centres none of the models with a
negative sample mean (the liberal bound); always `lower <= consistent <=
upper`. Each p-value is the fraction of the `reps` replicate statistics
above the observed one.
`block_size=None` uses the Politis-White (2004) / Patton-Politis-White
(2009) optimal length of each loss-differential column, averaged over the
columns and rounded (reported in `block_size`, with `block_size_auto`
True). The schemes are the library's stationary (geometric blocks, mean
`block_size`), circular-block and moving-block bootstraps, one Philox
substream per replication, bit-identical at any thread count.
Validation: `arch.bootstrap.SPA`/`RealityCheck` reproduced EXACTLY
(means, variances, every replicate statistic, p-values, critical values)
when the Rust core is fed arch's own resample indices, with
`studentize=False` — arch 8.0's `studentize` flag is inert (measured:
identical output on/off), so arch computes the un-studentized statistic
with Hansen's re-centrings; the studentized path is pinned at 1e-12
against a NumPy transcription of Hansen's formulas on the same
resamples. The public seeded path lands within 0.05 of arch at 4000
replications. Size and power are measured by seeded Monte Carlo: the
un-studentized rejection rates under the least favourable null match
arch's own on the same design, and power against a dominated benchmark
is 0.985 at 5%. `studentize=True` (the default, Hansen's statistic)
divides the observed statistic and every bootstrap replicate by the
SAME estimated omega_k, which over-rejects in short samples — measured
0.123 at a nominal 0.05 with n=200 and AR(0.5) losses, shrinking to
0.093 at n=800. No package computes that statistic, so it is measured,
not validated; prefer `studentize=False` in a short evaluation sample
(see the forecasting model card for the full tables).
Further arguments, with defaults: `block_size` (None: Politis-White),
`reps` (1000), `bootstrap` ("stationary"; or "circular",
"moving_block"), `studentize` (True), `nested` (False), `seed` (0).
Returned keys: `statistic` (sqrt(n)-scaled), `best_model` (index
attaining the maximum), `p_value` (= `p_value_consistent`),
`p_value_lower`, `p_value_consistent`, `p_value_upper`, `crit_levels`
([0.90, 0.95, 0.99]), `crit_lower`, `crit_consistent`, `crit_upper`
(critical values at those levels, same scale as `statistic`),
`mean_loss_diff` (`dbar_k`), `loss_diff_var` (`omega_k^2`),
`recentered` (bool per model: re-centred under the consistent p-value),
`boot_lower`, `boot_consistent`, `boot_upper` (the `reps` replicate
statistics), `n`, `m`, `block_size`, `block_size_auto`, `reps`,
`bootstrap`, `studentize`, `nested`.
model_confidence_set¶
def model_confidence_set(
losses: _ArrayLike,
size: float = ...,
method: str = ...,
block_size: int | None = ...,
reps: int = ...,
bootstrap: str = ...,
seed: int = ...,
) -> dict[str, Any]:
The Hansen-Lunde-Nason (2011) Model Confidence Set: which of m models
are statistically indistinguishable from the best?
`losses` is a `T x m` array with one loss column per model (`m >= 2`),
index-aligned over the same evaluation periods. Starting from all
models, each step tests equal predictive ability across the models still
in the set with the range statistic `T_R = max_{i,j} |dbar_ij| /
sqrt(var*(dbar_ij))` (`method="R"`, HLN's recommended default) or the max
statistic `T_max = max_i dbar_i. / sqrt(var*(dbar_i.))` (`method="max"`),
against a block bootstrap of the loss panel (the same resamples reused at
every step, re-centred at the sample means); the worst model — the row of
the maximizing pair under `T_R`, every model attaining the maximum under
`T_max` — is eliminated with the step's p-value, until one model remains.
A model's MCS p-value is the running maximum of the step p-values along
the elimination path, so the set at any `size` is `{k : p_MCS(k) >
size}` (`included`) and the sets are nested in `size`; the p-values do
not depend on `size`. Hansen-Lunde-Nason's guarantee — the set contains
the best model(s) with probability at least `1 - size` — is ASYMPTOTIC
and about the whole best set. Measured on a design with two
exactly-equally-best models it holds at about 0.87 against a nominal 0.90
and does not improve from n=150 to n=600, matching `arch`'s own rates on
the same design; the easier event "a best model is in the set" does hold
at the nominal level. The model card has the table.
`block_size=None` uses the Politis-White optimal length of each loss
column, averaged and rounded (`block_size`, `block_size_auto`). Identical
loss columns are degenerate under `method="R"` only — their pairwise
bootstrap variance is exactly zero, so the panel is refused by name —
while `method="max"` standardizes against the cross-sectional mean, gives
the duplicates one statistic and eliminates them together in one step.
`arch` 8.0.0 handles neither: measured, its `method="R"` raises after
warning about the 0/0 division and its `method="max"` does not return.
Validation: `arch.bootstrap.MCS` reproduced EXACTLY (mean losses,
elimination order, included/excluded sets, MCS p-values, the pairwise
variance matrix) when fed arch's own resample indices; the public seeded
path reproduces arch's set and lands within 0.05 of its p-values at 4000
replications; coverage of the best set is measured by seeded Monte
Carlo and cross-checked against arch's own frequencies on the same
design (see the forecasting model card).
Further arguments, with defaults: `size` (0.10), `method` ("R"; or
"max"), `block_size` (None: Politis-White), `reps` (1000), `bootstrap`
("stationary"; or "circular", "moving_block"), `seed` (0).
Returned keys: `included` (model indices in the set, ascending),
`excluded`, `mcs_p_values` (per model), `elimination_order` (every model
in the order eliminated, survivor last), `step_p_values` (the raw
p-value of the step that eliminated each model, aligned with
`elimination_order`; 1.0 for the survivor), `statistics` (observed `T_R`
/ `T_max` per step), `n_steps`, `mean_losses`, `n`, `m`, `size`,
`method`, `block_size`, `block_size_auto`, `reps`, `bootstrap`.
stepm_test¶
def stepm_test(
benchmark_losses: _ArrayLike,
model_losses: _ArrayLike,
size: float = ...,
block_size: int | None = ...,
reps: int = ...,
bootstrap: str = ...,
studentize: bool = ...,
nested: bool = ...,
seed: int = ...,
) -> dict[str, Any]:
The Romano-Wolf (2005) StepM procedure: WHICH models beat the benchmark,
controlling the family-wise error rate at size, on the SPA bootstrap.
Arguments as `spa_test`. Step 1 declares superior every model whose
statistic (`sqrt(n) dbar_k / omega_k`, or un-studentized) exceeds the
`1 - size` quantile of the bootstrap maximum over ALL models under the
consistent re-centring; each later step recomputes that quantile over
the models not yet declared superior and adds those now exceeding it,
until a step adds nothing (or every model is superior — `arch` raises
there; this stops). The full-set SPA result is returned alongside.
Validation: reproduces `arch.bootstrap.StepM`'s superior set exactly on
arch's own resamples (`studentize=False`), the studentized rule against
the NumPy transcription.
Further arguments, with defaults: `size` (0.05), `block_size` (None:
Politis-White), `reps` (1000), `bootstrap` ("stationary"), `studentize`
(True), `nested` (False), `seed` (0).
Returned keys: `superior_models` (indices, ascending), `n_superior`,
`steps` (models declared superior at each step; the last entry is empty
when the procedure stopped because a step added nothing), `n_steps`,
`step_crit_values` (the `1 - size` bootstrap quantile compared against at
each step, on the `statistic` scale), `size`, plus every key of
`spa_test` for the full-set test: `statistic`, `best_model`, `p_value`,
`p_value_lower`, `p_value_consistent`, `p_value_upper`, `crit_levels`,
`crit_lower`, `crit_consistent`, `crit_upper`, `mean_loss_diff`,
`loss_diff_var`, `recentered`, `boot_lower`, `boot_consistent`,
`boot_upper`, `n`, `m`, `block_size`, `block_size_auto`, `reps`,
`bootstrap`, `studentize`, `nested`.
Exponential smoothing (ETS)¶
ets_fit¶
def ets_fit(
y: _ArrayLike,
error: str = ...,
trend: str | None = ...,
damped: bool = ...,
seasonal: str | None = ...,
seasonal_periods: int | None = ...,
initialization: str = ...,
horizon: int = ...,
level: float | None = ...,
n_sim: int | None = ...,
seed: int | None = ...,
optimizer: str | None = ...,
smoothing_params: Sequence[float] | None = ...,
initial_states: Sequence[float] | None = ...,
max_iter: int | None = ...,
) -> dict[str, Any]:
Innovations state-space exponential smoothing — one member of the ETS(Error, Trend, Seasonal) taxonomy of Hyndman, Koehler, Snyder & Grose (2002) / Hyndman et al. (2008), fitted by maximum likelihood.
`error` is "add" or "mul"; `trend` and `seasonal` are None, "add" or
"mul"; `damped=True` damps the trend (estimates `phi`; refused without
a trend, where it would be inert); `seasonal_periods` is the period m
(12 monthly, 4 quarterly), required with a seasonal component and
refused without one. ETS(A,N,N) is simple exponential smoothing,
(A,A,N) Holt, (A,Ad,N) the damped trend, (A,A,A) / (M,A,M) the
additive / multiplicative Holt-Winters. Any multiplicative component
needs strictly positive `y` (refused otherwise, naming the offending
observation). NaN is refused: the innovations form conditions every
state on the observed error and has no missing-value mechanism
(R's ets refuses NaN too) — interpolate first, or use a Kalman-filter
model.
The smoothing parameters are Hyndman's `alpha`, `beta`, `gamma`, `phi`
(not beta* = beta/alpha or gamma* = gamma/(1 - alpha)), searched in the
traditional box 0 < alpha < 1, 0 < beta < alpha, 0 < gamma < 1 - alpha,
0.8 <= phi <= 0.98 (R's and statsmodels' default bounds).
`initialization`: "estimated" (default: the initial states are free
parameters, started from the heuristic; the seasonal indices are
normalised to sum to zero / average one, R's convention, and count
m - 1 free parameters), "heuristic" (the Hyndman 2008 section 2.6.1
heuristic — a centred moving average over the first cycles and a
linear regression of its first ten values — held fixed; needs at
least 10 observations and, with a seasonal, 2m and 10 + 2 floor(m/2)),
or "known" (fixed at `initial_states`). `initial_states` is
`[level, trend?, seasonal[0], ..., seasonal[m-1]]` where
`seasonal[j]` is the index in force for observation j; it is required
with "known" and refused with the other two. `smoothing_params`
(`[alpha, beta?, gamma?, phi?]`, the components present, in that
order) evaluates the model at FIXED parameters with no optimisation —
statsmodels' `smooth(params)` — and needs initialization "heuristic"
or "known" (with "estimated" nothing would be estimated: refused);
`optimizer` and `max_iter` are then inert and refused if passed.
`optimizer` is "auto" (the effective default: L-BFGS and Nelder-Mead
from a staged start — the smoothing parameters alone at the heuristic
states first — then a BFGS polish of whichever did better; the
returned `optimizer` key reads "nelder_mead+bfgs"), "nelder_mead",
"bfgs" or "lbfgs"; `max_iter` caps each stage's iterations.
`horizon=h` adds h-step forecasts with `level` (0.95 when omitted)
prediction intervals: for the class-1 models — additive error with
additive or no trend and seasonal — the exact Gaussian intervals from
the closed-form variances of Hyndman et al. (2008, Table 6.1)
(`interval_method="exact"`); for every other model `n_sim` (5000 when
omitted) seeded innovation paths through the fitted recursion, the
bounds being empirical quantiles (`"simulated"`; `seed` 0 when
omitted). `level`, `n_sim` and `seed` are refused with `horizon=0`,
and `n_sim` / `seed` are refused for a class-1 model, where nothing is
simulated. The point forecast is always the zero-innovation path (R's
and statsmodels' convention). Allocation guards, not modelling limits:
`horizon` may not exceed 1000000, and `n_sim * horizon` (the simulated
values held at once) may not exceed 2^28; both are refused by name.
Returned keys: `spec` (e.g. "ETS(A,Ad,N)"), `short_name` ("AAdN"),
`error`, `trend`, `damped`, `seasonal`, `seasonal_periods` (None
without a seasonal), `alpha`, `beta`, `gamma`, `phi` (None when the
component is absent), `params` and `param_names` (the packed
smoothing vector), `initial_level`, `initial_trend`,
`initial_seasonal` (None when absent), `initial_states` and
`initial_state_names` (packed), `initialization`, `fitted`
(one-step-ahead), `resid` (`y - fitted`, or `(y - fitted) / fitted`
under multiplicative errors), `level_path`, `trend_path`,
`seasonal_path` (the states after each update; None when absent),
`final_level`, `final_trend`, `final_seasonal` (the forecast anchor;
`final_seasonal[j]` is the index for forecast step j), `final_states`,
`loglik` (the concentrated Gaussian log-likelihood, statsmodels'
convention; R's `ets` omits the constant -(n/2)(ln(2 pi / n) + 1)),
`sigma2` (`mean(resid^2)`), `nobs`, `k_params` (smoothing parameters +
free initial states under "estimated" + sigma2), `aic`, `aicc`, `bic`,
`converged`, `n_iterations`, `n_fevals`, `optimizer` ("none" at fixed
parameters), `class1`, `horizon`, and — None when `horizon=0` —
`forecast`, `forecast_variance`, `forecast_lower`, `forecast_upper`,
`interval_level`, `interval_method`, `n_sim`, `seed`.
Validation (fixtures/ets.json): the twenty models without a
multiplicative seasonal are pinned at fixed parameters to statsmodels
`ETSModel` (log-likelihood, fitted values, states, forecasts,
simulations) at 1e-10, the six class-1 forecast variances to
statsmodels' exact `get_prediction` and to the Table 6.1 closed forms
at 1e-10 / 1e-12, the heuristic initialisation to
`holtwinters.ExponentialSmoothing` at 1e-10, and the maximum likelihood
to statsmodels' L-BFGS-B fit (match-or-beat, two optimizers); the ten
multiplicative-seasonal models are pinned at 1e-12 to an independent
transcription of the published recursion and their simulator to
statsmodels `simulate` at 1e-10 — statsmodels' own smoother uses the
classical Holt-Winters seasonal update there (a stated, measured
convention gap). Interval coverage and parameter recovery are measured
by seeded Monte Carlo and quoted on the model card.
Further arguments, with defaults: `error` ("add"), `trend` (None),
`damped` (False), `seasonal` (None), `seasonal_periods` (None),
`initialization` ("estimated"), `horizon` (0), `level` (None: 0.95
when forecasting), `n_sim` (None: 5000 when simulating), `seed` (None:
0 when simulating), `optimizer` (None: "auto"), `smoothing_params`
(None: estimated), `initial_states` (None), `max_iter` (None).
auto_ets¶
def auto_ets(
y: _ArrayLike,
seasonal_periods: int | None = ...,
ic: str = ...,
allow_multiplicative_trend: bool = ...,
restrict: bool = ...,
damped: bool | None = ...,
initialization: str = ...,
horizon: int = ...,
level: float | None = ...,
n_sim: int | None = ...,
seed: int | None = ...,
optimizer: str | None = ...,
) -> dict[str, Any]:
Automatic ETS model selection — the candidate-set search of Hyndman
et al. (2008, section 7.2) as R's forecast::ets runs it: every
admissible member of the taxonomy is fitted by maximum likelihood
(ets_fit) and the one with the smallest information criterion
(ic: "aicc" default, "aic", "bic") is returned, fitted, with the
ranked candidate table.
Candidates: error "add"/"mul", trend None/"add" (plus "mul" with
`allow_multiplicative_trend=True`; R's default excludes it), damped
and undamped trends (`damped=None`; `True`/`False` restricts to one),
seasonal None/"add"/"mul" when `seasonal_periods` >= 2 (None: non-
seasonal candidates only). Multiplicative errors and seasonals are
tried only when every `y` > 0. `restrict=True` (R's default) drops the
combinations with infinite forecast variance or a mis-scaled error —
additive error with any multiplicative component, and (M,M,A) — so a
positive seasonal series has 15 candidates (6 additive-error, 9
multiplicative-error), a non-positive seasonal one 6, a non-seasonal
positive series 6, and a non-seasonal non-positive series 3.
`initialization` is "estimated" or "heuristic" for every
candidate; `optimizer` and the forecast options (`horizon`, `level`,
`n_sim`, `seed`) are those of `ets_fit` — `n_sim` and `seed` act only
if the selected model is not class 1 (the winner is not known in
advance, so they are accepted regardless; with `horizon=0` they are
refused as inert), including their allocation guards (`horizon`
at most 1000000, `n_sim * horizon` at most 2^28). NaN is refused.
Returned keys: every key of `ets_fit` for the selected model (its very
fit from the search, not a refit — refitting reproduces it exactly),
plus `ic`, `ic_value`, `candidates` — a list of dicts with `spec`,
`short_name`, `loglik`, `aic`, `aicc`, `bic`, `ic_value`, `k_params`,
`converged`, `status` ("ok" or "error") and `error` (the message when
a candidate failed; failures never abort the search) ranked by the
criterion, failures last — `n_candidates` and `n_fitted`. Read the
table: candidates within ~2 of the best criterion are near-ties the
data do not distinguish.
Validation (honest grade, as for `auto_arima`): every candidate's
likelihood is the golden-pinned `ets_fit` likelihood; the candidate
set reproduces R's `forecast::ets` enumeration exactly
(fixtures/ets.json); the selection loop itself has no runnable
third-party reference (the M3 forecast-competition parity of the
method is R-only), so it is graded by seeded Monte-Carlo recovery of
the generating component form, quoted on the model card.
Further arguments, with defaults: `seasonal_periods` (None), `ic`
("aicc"), `allow_multiplicative_trend` (False), `restrict` (True),
`damped` (None: both), `initialization` ("estimated"), `horizon` (0),
`level` (None: 0.95), `n_sim` (None: 5000), `seed` (None: 0),
`optimizer` (None: "auto").
Returned keys: `aic`, `aicc`, `alpha`, `beta`, `bic`, `candidates`,
`class1`, `converged`, `damped`, `error`, `final_level`,
`final_seasonal`, `final_states`, `final_trend`, `fitted`,
`forecast`, `forecast_lower`, `forecast_upper`, `forecast_variance`,
`gamma`, `horizon`, `ic`, `ic_value`, `initial_level`,
`initial_seasonal`, `initial_state_names`, `initial_states`,
`initial_trend`, `initialization`, `interval_level`,
`interval_method`, `k_params`, `level_path`, `loglik`,
`n_candidates`, `n_fevals`, `n_fitted`, `n_iterations`, `n_sim`,
`nobs`, `optimizer`, `param_names`, `params`, `phi`, `resid`,
`seasonal`, `seasonal_path`, `seasonal_periods`, `seed`,
`short_name`, `sigma2`, `spec`, `trend`, `trend_path` — every
`ets_fit` key for the selected model, read there, plus the five
selection extras above.
structural time-series models (unobserved components) and TVP regression¶
unobserved_components¶
def unobserved_components(
y: np.ndarray,
level: str = "llevel",
seasonal: int | None = None,
stochastic_seasonal: bool | None = None,
freq_seasonal: list[float] | None = None,
freq_seasonal_harmonics: list[int] | None = None,
stochastic_freq_seasonal: list[bool] | None = None,
cycle: bool = False,
damped_cycle: bool | None = None,
stochastic_cycle: bool | None = None,
cycle_period_bounds: list[float] | None = None,
exog: np.ndarray | None = None,
forecast_steps: int = 0,
forecast_exog: np.ndarray | None = None,
fixed_params: list[float] | None = None,
n_starts: int = 3,
) -> dict[str, Any]:
Harvey's structural time-series ("unobserved components") models by exact-diffuse maximum likelihood — level/trend, dummy and trigonometric seasonals, a (damped) stochastic cycle, and regressors:
y_t = mu_t + gamma_t + c_t + beta' x_t + eps_t
with every state initialized exactly diffuse (Koopman 1997), NaN in `y`
treated as missing, and the components assembled exactly as statsmodels'
`UnobservedComponents(..., use_exact_diffuse=True)` enumerates them (the
log-likelihoods are directly comparable).
`level` picks the level/trend block by its statsmodels name (long or
short form): "irregular"/"ntrend", "fixed intercept", "deterministic
constant"/"dconstant", "local level"/"llevel" (default), "random
walk"/"rwalk", "fixed slope", "deterministic trend"/"dtrend", "local
linear deterministic trend"/"lldtrend", "random walk with
drift"/"rwdrift", "local linear trend"/"lltrend", "smooth
trend"/"strend", "random trend"/"rtrend". `seasonal=s` adds an
`s-1`-state dummy seasonal (`stochastic_seasonal`, default True, gives
it a variance; False makes it fixed dummies). `freq_seasonal=[p, ...]`
adds trigonometric seasonals with `freq_seasonal_harmonics` harmonics
each (default `floor(p/2)`) and `stochastic_freq_seasonal` flags
(default True each). `cycle=True` adds the stochastic cycle:
`damped_cycle` (default False) estimates a damping in (0, 1),
`stochastic_cycle` (default False) gives it a variance, and
`cycle_period_bounds=[min, max]` confines its frequency to
`(2 pi/max, 2 pi/min)`. Its default is `[2, len(y)]`, and an infinite
`max` is read the same way: under exact-diffuse initialization the
log-likelihood of a stochastic cycle DIVERGES as the frequency goes to
zero (the second cycle state becomes weakly observable and its diffuse
resolution contributes `-ln(lambda)`), so an unbounded period is not a
safe search region — and a cycle longer than the sample is not
identified in any case. statsmodels leaves that bound at infinity.
`exog` (T x k) enters with
time-invariant coefficients estimated jointly (statsmodels
`mle_regression=True`); `forecast_steps=h` returns h-step forecasts and
needs `forecast_exog` (h x k) when `exog` is given. `fixed_params`
(statsmodels order: `sigma2.irregular`, the state variances in
component order, `frequency.cycle`, `damping.cycle`, `beta.x1`...)
evaluates the model there instead of estimating. `n_starts` (default 3)
is the deterministic start ladder. Counts that size an allocation are
bounded and refuse rather than abort: `seasonal` and each
`freq_seasonal` period at most `len(y)` (they cost states, and a period
the sample never completes is not identified), `forecast_steps` at most
100000, `n_starts` at most 64.
Options that act only under a component RAISE when passed without it:
`stochastic_seasonal` without `seasonal`; `freq_seasonal_harmonics` /
`stochastic_freq_seasonal` without `freq_seasonal`; `damped_cycle` /
`stochastic_cycle` / `cycle_period_bounds` without `cycle=True`;
`forecast_exog` without `exog` or without `forecast_steps`.
Estimation: BFGS + Nelder-Mead on the exact prediction-error
log-likelihood in statsmodels' square-root/logistic working space,
scale-adaptive (y standardized, mapped back exactly). A variance whose
estimate cannot be told from zero (zeroing it costs < 1e-4
log-likelihood — the pile-up) is flagged in `at_boundary` with a NaN
standard error; `se` are observed-information (numerical Hessian)
CONDITIONAL on the flagged parameters sitting exactly at their boundary,
i.e. the information matrix is inverted over the free parameters only.
statsmodels' `cov_type="approx"` inverts the full matrix instead,
including the boundary directions where it is indefinite, so the two
agree exactly when nothing is flagged and differ by definition when
something is.
`aic`/`bic` use statsmodels' `k_params + k_diffuse` degrees of freedom.
The component keys without a prefix are the SMOOTHED (two-sided) paths;
`filtered_*` are the one-sided ones. Variances inside the diffuse period
are the finite part; `std_resid` is NaN there and at missing periods.
Validation (honest grade): fixed-parameter log-likelihood, filtered
states and variances, smoothed states, residuals, forecasts and forecast
variances pinned at 1e-8 against statsmodels for 26 component
combinations and NaN-inserted series (fixtures/uc.json); the SMOOTHED
variances at 1e-8 too, except inside the diffuse period of the hardest
combinations, where the exact-diffuse smoother is ill-conditioned and
the tolerance is the distance between statsmodels' own two smoother
implementations (up to 4.5e-3 there, against filters that agree to
2.9e-11); the MLE pinned
to the better of statsmodels' own fit and a SciPy re-optimization of the
identical criterion (two optimizers); the Durbin-Koopman (2012) Nile
local level reproduced to the book's printed precision; the
Harvey-Durbin (1986) UK seat-belt BSM re-estimated to that same optimum,
with its slope and seasonal variance pile-ups flagged (the series is
fetched from Rdatasets when the fixture is generated and never
redistributed, so only its derived optimum is stored); parameter
recovery, forecast-interval coverage and scale invariance measured by
seeded Monte Carlo (see the model card).
Returned keys: `trend_specification`, `param_names`, `params`, `se`,
`at_boundary`, `loglik`, `aic`, `bic`, `nobs`, `nobs_observed`,
`nobs_diffuse`, `k_states`, `k_diffuse`, `k_params`, `estimated`,
`converged`, `n_iter`, `n_fevals`, `state_names`, `filtered_state`,
`filtered_state_var`, `smoothed_state`, `smoothed_state_var` (nested
lists, nobs x k_states), `level`, `level_var`, `filtered_level`,
`filtered_level_var`, `slope`, `slope_var`, `filtered_slope`,
`filtered_slope_var`, `seasonal`, `seasonal_var`, `filtered_seasonal`,
`filtered_seasonal_var`, `cycle`, `cycle_var`, `filtered_cycle`,
`filtered_cycle_var` (None when the component is absent),
`freq_seasonal`, `freq_seasonal_var`, `filtered_freq_seasonal`,
`filtered_freq_seasonal_var` (lists with one array per block), `fitted`,
`resid`, `std_resid`, `forecast`, `forecast_var`.
tvp_regression¶
def tvp_regression(
y: np.ndarray,
x: np.ndarray,
constant: bool = True,
fixed_params: list[float] | None = None,
n_starts: int = 3,
) -> dict[str, Any]:
Regression with random-walk (time-varying) coefficients by exact-diffuse maximum likelihood, with the pile-up check:
y_t = x_t' beta_t + eps_t, beta_{t+1} = beta_t + eta_t,
eps_t ~ N(0, sigma2_eps), eta_t ~ N(0, diag(sigma2_beta)),
beta_1 diffuse.
`x` is T x k (`constant=True`, the default, prepends a random-walk
intercept). NaN in `y` is a missing period; `x` must be finite. The
observation variance and the k coefficient-innovation variances are
estimated by BFGS + Nelder-Mead on the exact prediction-error
log-likelihood (square-root working space, scale-adaptive, `n_starts`
deterministic starts, default 3, at most 64); `fixed_params=[sigma2_eps,
sigma2_beta_1, ..., sigma2_beta_k]` evaluates the filter there instead —
with every state variance 0 it is recursive least squares (statsmodels
`RecursiveLS`), the expanding-window OLS path.
The pile-up problem (Shephard-Harvey 1990; Stock-Watson 1998): a
random-walk-coefficient variance whose MLE cannot be told from zero
(zeroing it costs < 1e-4 log-likelihood) is flagged in `pile_up` (and
`at_boundary`) and gets a NaN standard error, so a coefficient the data
cannot show moving is not reported as moving by a tiny amount. `se` are
observed-information (numerical Hessian) conditional on the flagged
variances being exactly zero; `aic`/`bic` use statsmodels'
`k_params + k_diffuse` degrees of freedom; `std_resid` is NaN inside the
diffuse period (the first k informative observations) and at missing
periods.
Validation (honest grade): fixed-parameter log-likelihood, filtered and
smoothed coefficient paths and variances, residuals pinned at 1e-8
against a statsmodels `MLEModel` transcription of the documented
state-space form and against `RecursiveLS` in the zero-variance limit
(its filtered coefficients and concentrated log-likelihood); the MLE
pinned to the better of statsmodels' fit and a SciPy re-optimization of
the same criterion (two optimizers) with the true-zero variance's
pile-up flagged; the pile-up frequency on constant versus moving
coefficients measured by seeded Monte Carlo (model card).
Returned keys: `coef_names`, `k`, `param_names`, `params`, `se`,
`at_boundary`, `sigma2_eps`, `sigma2_beta`, `pile_up`, `loglik`, `aic`,
`bic`, `nobs`, `nobs_observed`, `nobs_diffuse`, `k_params`, `estimated`,
`converged`, `n_iter`, `n_fevals`, `beta_filtered`, `beta_filtered_var`,
`beta_smoothed`, `beta_smoothed_var` (nested lists, nobs x k), `fitted`,
`resid`, `std_resid`.