statspai.timeseries

timeseries

Time series methods: models for causal inference contexts and forecasting.

Provides VAR (vector autoregression), structural break tests, Granger causality, unit-root tests and cointegration analysis; and a forecasting toolkit -- ARIMA, exponential smoothing (ETS), benchmark methods, STL decomposition, forecast accuracy, rolling-origin cross-validation and hierarchical reconciliation.

StateDependentLPResult dataclass

Bases: ResultProtocolMixin

Regime-specific impulse responses from a state-dependent projection.

*_state1 is the regime F = 1 and *_state0 the regime F = 0; diff = irf_state1 - irf_state0 with its standard error and two-sided p-value (missing for the Cholesky projection, whose bands are not standard errors).

to_frame

to_frame() -> DataFrame

One row per regime and horizon, regime state1 first.

difference_frame

difference_frame() -> DataFrame

Response in F = 1 minus response in F = 0, by horizon.

ARDLResult

Bases: ResultProtocolMixin

A fitted AR(p) or ADL(p, q) forecasting regression.

Attributes:

Name	Type	Description
model	EconometricResults	The OLS fit. Pass it to sp.test, sp.lincom and so on; lagged regressors are named <var>_L<k>.
params, std_errors	Series	Coefficients and standard errors of model.
lags	int	Autoregressive order p.
x_lags	dict	Lag order of each additional regressor.
omitted, omitted_x	list of str	Terms left out of the fit because of collinearity, and the regressors that lost every term that way (they are not in x or x_lags). Both are empty for a full-rank design.
nobs	int	Observations in the estimation sample.
ser	float	Standard error of the regression, sqrt(SSR / (T - K)).
bic, aic	float	ln(SSR / T) + K ln(T) / T and ln(SSR / T) + 2 K / T, with K the number of coefficients including the intercept.
ic_table	DataFrame or None	The criteria for every candidate order, when lags was selected.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> y = np.zeros(200)
>>> for t in range(1, 200):
...     y[t] = 1.0 + 0.5 * y[t - 1] + rng.normal()
>>> res = sp.ardl(pd.DataFrame({"y": y}), "y", lags=1)
>>> type(res).__name__
'ARDLResult'
>>> list(res.params.index)
['Intercept', 'y_L1']
>>> res.forecast().shape
(1, 4)

forecast

forecast(alpha: Optional[float] = None, *, steps: int = 1, exog: Optional[DataFrame] = None) -> DataFrame

Forecasts of the periods after the last estimation observation.

Beyond one step the forecast is dynamic: each forecast of y stands in for the lag the next one needs.

Parameters:

Name	Type	Description	Default
alpha	float	Level of the interval. Default: the model's.	None
steps	int	Number of periods ahead.	1
exog	DataFrame	Future values of the x series, one row per period after the sample in order, with the columns named as in the fit. Needed whenever a forecast uses an x dated after the sample: from the first step with contemporaneous=True, from the second otherwise.	None

Returns:

Type	Description
DataFrame	One row per period with forecast, rmsfe, lower, upper.

Notes

The interval is forecast +/- z * RMSFE. At one step the RMSFE is the standard error of the regression. At step h it is that times sqrt(1 + psi_1**2 + ... + psi_{h-1}**2), where psi are the moving-average weights of the autoregressive part. It ignores the estimation error in the coefficients, treats the future x as known and assumes normal errors; a pseudo out-of-sample RMSFE from :meth:poos is the more honest width at one step.

long_run

long_run(alpha: Optional[float] = None) -> DataFrame

Long-run effect of each x on y.

A permanent unit rise in x moves y in the end by the sum of the coefficients on x and its lags divided by one minus the sum of the autoregressive coefficients. In the geometric-lag model y[t] = a + b x[t] + lam y[t-1] this is b / (1 - lam).

Returns:

Type

Description

DataFrame

One row per x in the model, and Intercept for the long-run mean of y at x = 0 when the model has a constant, with coef, std err, z, pvalue and the confidence limits. The standard errors come from the delta method on the model's covariance matrix, so they are robust when the fit is. attrs['ar_sum'] is the sum of the autoregressive coefficients and attrs['adjustment'] that sum minus one, the error-correction coefficient.

Notes

The ratio has no finite-sample moments and its normal approximation is poor when the sum of the autoregressive coefficients is near one. With that sum at or above one there is no long-run level to converge to and the method refuses.

granger

granger(x: Optional[str] = None) -> Dict[str, Any]

F test that every lag of x has a zero coefficient.

The Granger-causality statistic: whether lags of x help predict y given the lags of y. It uses the model's own covariance estimator, so it is heteroskedasticity- or autocorrelation-robust when the fit is. "Causality" here means predictive content.

poos

poos(start: Any, window: str = 'expanding') -> DataFrame

Pseudo out-of-sample one-step forecasts.

For each date from start to the end of the data the model is re-estimated on observations before that date and used to forecast it, as a forecaster standing at the previous date would have done.

Parameters:

Name	Type	Description	Default
start	label or int	First date to forecast: a label of the data's index (or of the time column), or a row position.	required
window	('expanding', 'rolling')	"expanding" re-estimates on every observation from the start of the estimation sample; "rolling" keeps the window length fixed at its initial size.	"expanding"

Returns:

Type	Description
DataFrame	actual, forecast and error (actual minus forecast) per forecast date. attrs holds rmsfe (root mean squared forecast error), bias (mean error), bias_se and n_forecasts. A bias several standard errors from zero says the forecasts were systematically off, which a break in the relationship would produce.

ARIMAResult dataclass

Bases: ResultProtocolMixin

Fitted ARIMA(p,d,q) / SARIMAX model returned by :func:statspai.arima.

Carries the estimated params / se (indexed by parameter name), information criteria (aic / bic / aicc), the log-likelihood, residuals and fitted values, plus inference accessors (:attr:tvalues, :attr:pvalues, :meth:conf_int) and a .summary() / .forecast() interface.

Examples:

>>> import statspai as sp
>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> n = 120
>>> y = np.zeros(n)
>>> for t in range(1, n):
...     y[t] = 0.6 * y[t - 1] + rng.normal(0, 1)
>>> res = sp.arima(y, order=(1, 0, 0))
>>> type(res).__name__
'ARIMAResult'
>>> bool(np.isfinite(res.aic))
True
>>> list(res.params.index)
['const', 'ar.L1', 'sigma2']

period property

period: int

Seasonal period (1 for a non-seasonal model).

n_arma_params property

n_arma_params: int

Number of estimated ARMA coefficients, p + q + P + Q: the degrees of freedom a portmanteau test of the residuals gives up.

sigma2 property

sigma2: float

Maximum likelihood estimate of the innovation variance.

std_errors property

std_errors: Series

Alias for :attr:se (regression-style naming).

tvalues property

tvalues: Series

z-statistics params / se (SARIMAX uses a normal reference).

pvalues property

pvalues: Series

Two-sided p-values from the normal reference distribution.

conf_int

conf_int(alpha: float = 0.05) -> DataFrame

Confidence intervals for each parameter.

Parameters:

Name	Type	Description	Default
alpha	float	1 - alpha is the coverage (0.05 → 95% CI).	0.05

Returns:

Type	Description
DataFrame	Indexed by parameter name with lower / upper columns.

forecast

forecast(horizon: int = 10, alpha: float = 0.05, *, level: Optional[Any] = None, exog: Optional[Any] = None, dof_adjust: bool = False) -> DataFrame

Forecasts with prediction intervals.

Parameters:

Name	Type	Description	Default
horizon	int	Number of periods ahead.	10
alpha	float	One minus the coverage of the interval, used when level is not given; the columns are then forecast, lower, upper.	0.05
level	float or sequence of float	Coverage in percent, e.g. (80, 95); the columns are then forecast, lower_80, upper_80, lower_95, upper_95, the layout of :func:statspai.ets.	None
exog	array - like or DataFrame	Values of the regressors over the forecast periods, one row per period, columns in the order of the fit (a DataFrame is matched by column name). Required when the model has regressors: the forecast is conditional on them.	None
dof_adjust	bool	Scale the innovation variance by n / (n - k), n the number of observations after differencing and k the number of estimated coefficients, as R's forecast::Arima and fable::ARIMA do. The default uses the maximum likelihood variance, as stats::arima and Stata do.	False

Returns:

Type	Description
DataFrame	Indexed by the forecast periods when the series had a date or period index.

Notes

The intervals treat the coefficients, and the future regressors, as known.

BaggedForecastResult dataclass

Bases: ResultProtocolMixin

Bagged forecasts returned by :func:statspai.bagged_forecast.

Attributes:

Name	Type	Description
forecast	DataFrame	forecast (the average over the bootstrapped series) with ensemble_min and ensemble_max, the range of the members. The range describes how much the forecast depends on the sample; it is not a prediction interval.
members	DataFrame	The forecast from each bootstrapped series, one column each.
forecaster	str
n_boot	int

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(1)
>>> y = 100 + np.cumsum(rng.normal(0.5, 1.0, 60))
>>> bag = sp.bagged_forecast(y, "drift", horizon=4, n_boot=20)
>>> bag.forecast.shape, bag.members.shape
((4, 3), (4, 20))

BeveridgeNelsonResult dataclass

Bases: ResultProtocolMixin

Decomposition returned by :func:statspai.beveridge_nelson.

Attributes:

Name	Type	Description
trend, cycle	Series	Permanent and transitory components, aligned with the input; trend + cycle is the series. With method == 'ols' the first order values are missing, with 'mle' the first one.
drift	float	Mean of the first difference, intercept / (1 - sum(ar_coefs)).
order	int	Autoregressive order of the model for the first difference.
ma_order	int	Moving-average order; zero for an autoregression.
intercept	float	Constant of the autoregressive form, drift * (1 - sum(ar_coefs)).
ar_coefs, ma_coefs	ndarray	phi_1..phi_p and theta_1..theta_q of phi(L) (dy - drift) = theta(L) e with theta(L) = 1 + theta_1 L + ....
sigma2	float	Innovation variance: RSS / (N - order - 1) under 'ols', the maximum-likelihood estimate under 'mle'.
long_run_multiplier	float	psi(1) = theta(1) / phi(1): the long-run response of the level to a unit innovation. The trend moves by drift + psi(1) * innovation each period.
variance_ratio	float	Variance of the change in the trend over the variance of the change in the series implied by the fitted model, psi(1)^2 / sum(psi_j^2).
residuals	Series	Innovations of the model, aligned with the input: least-squares residuals under 'ols', one-step prediction errors of the exact likelihood under 'mle'.
selection	DataFrame or None	Information criteria by order, when the order was chosen.
n_obs	int	Number of first differences the model was fitted to.
method	str	'ols' (autoregression by least squares) or 'mle' (ARMA by exact Gaussian maximum likelihood).
loglik	float or None	Log-likelihood under 'mle'.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> y = np.cumsum(0.5 + rng.normal(size=200))    # random walk with drift
>>> bn = sp.beveridge_nelson(y, order=0)
>>> float(abs(bn.cycle).max())
0.0

summary

summary() -> str

Plain-text report of the fitted model and the decomposition.

plot

plot(ax: Any = None) -> Any

The series with its trend (top) and the cycle (bottom).

BVARResult dataclass

Bases: ResultProtocolMixin

Posterior summary of a Bayesian VAR with Minnesota prior.

Produced by :func:bvar. Holds the posterior-mean coefficient matrix and residual covariance, and exposes .forecast(), .irf() and .summary().

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> T = 120
>>> y1 = np.zeros(T); y2 = np.zeros(T)
>>> for t in range(1, T):
...     y1[t] = 0.5 * y1[t - 1] + 0.2 * y2[t - 1] + rng.normal()
...     y2[t] = 0.3 * y2[t - 1] + 0.1 * y1[t - 1] + rng.normal()
>>> df = pd.DataFrame({"gdp": y1, "inflation": y2})
>>> res = sp.bvar(df, lags=2)
>>> type(res).__name__
'BVARResult'
>>> res.coef.shape
(5, 2)
>>> res.lags
2
>>> isinstance(res.summary(), str)
True

coef_names property

coef_names: list

Row labels of coef: L<lag>.<variable> then _cons.

irf

irf(shock_var: int = 0, horizon: int = 20) -> ndarray

Orthogonalised impulse responses (Cholesky decomposition).

credible_interval

credible_interval(level: float = 0.9) -> tuple[ndarray, ndarray]

Posterior credible interval for every coefficient.

Returns (lower, upper) matrices the same shape as coef. With the error covariance fixed at Sigma0 the coefficient posterior is exactly normal; coef_sd holds its marginal standard deviations.

coef_table

coef_table() -> DataFrame

Posterior means, one column per equation, rows labelled.

CointegrationResult

Bases: ResultProtocolMixin

Results from cointegration test.

Produced by :func:engle_granger (and the Johansen routine). Exposes the test statistic(s), critical values, estimated cointegration rank and a formatted .summary().

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> T = 200
>>> x = np.cumsum(rng.normal(size=T))      # random walk
>>> y = 2.0 * x + rng.normal(size=T)       # cointegrated with x
>>> df = pd.DataFrame({"y": y, "x": x})
>>> res = sp.engle_granger(df, variables=["y", "x"])
>>> type(res).__name__
'CointegrationResult'
>>> res.test_type
'Engle-Granger'
>>> res.n_vars
2
>>> isinstance(res.summary(), str)
True

References

[@engle1987integration]

DLMResult dataclass

Bases: ResultProtocolMixin

Result of :func:dlm.

Attributes:

Name	Type	Description
variances	DataFrame	The observation variance (obs) and the state variance of every coefficient (state:<term>): the estimate, and for method='gibbs' the posterior sd, interval and effective sample size.
smoothed	DataFrame	The coefficient paths given the whole sample, one column per term, with <term>_sd (and <term>_lower / <term>_upper).
filtered	DataFrame	The coefficient paths given the data up to each date.
fitted	Series	One-step-ahead forecasts of y.
loglik	float	Log likelihood at the variances (maximum likelihood), or at their posterior means (Gibbs).
params	Series	The smoothed coefficients at the last date.
draws	DataFrame or None	Gibbs draws of the variances.
method, n_obs, level

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> T = 150
>>> x = rng.normal(size=T)
>>> slope = 1 + np.cumsum(0.1 * rng.normal(size=T))
>>> df = pd.DataFrame({"x": x, "y": 0.5 + slope * x + 0.3 * rng.normal(size=T)})
>>> fit = sp.dlm("y ~ x", df)
>>> isinstance(fit, sp.DLMResult)
True
>>> list(fit.smoothed.columns)[:2]
['Intercept', 'x']
>>> fit.forecast(2, pd.DataFrame({"x": [0.0, 1.0]})).shape
(2, 3)

forecast

forecast(steps: int = 1, data: Optional[DataFrame] = None) -> DataFrame

Forecast y beyond the sample.

Parameters:

Name	Type	Description	Default
steps	int	Horizon.	1
data	DataFrame	Regressor values for the forecast dates (steps rows). Not needed for a model with only an intercept.	None

Returns:

Type	Description
DataFrame	mean, lower and upper: the coefficients stay at their last filtered mean and their uncertainty grows by the state variances each step. The interval is for the outcome (it includes the observation variance). With method='gibbs' the variances are at their posterior means.

plot

plot(terms: Optional[Sequence[str]] = None) -> Any

Smoothed coefficient paths with their intervals.

TSCVResult dataclass

Bases: ResultProtocolMixin

Rolling-origin forecast errors returned by :func:statspai.tscv.

Attributes:

Name	Type	Description
errors	DataFrame	actual - forecast; one row per forecast origin (labelled by the last observation of its training window), one column per horizon h=1..H. nan where the outcome lies beyond the sample or the forecaster failed.
forecasts, actuals	DataFrame	Laid out like errors.
failures	list of (origin, message)	Origins where the forecaster raised.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> y = 100 + np.cumsum(rng.normal(0.2, 1.0, 60))
>>> cv = sp.tscv(y, "drift", horizon=2, initial=20)
>>> cv.errors.shape
(40, 2)
>>> list(cv.accuracy().columns)
['RMSE', 'MAE', 'MAPE', 'n']

accuracy

accuracy() -> DataFrame

RMSE, MAE and MAPE by forecast horizon, over the origins.

GARCHResult dataclass

Bases: ResultProtocolMixin

Fitted GARCH(p,q) model returned by :func:garch.

Holds the conditional-variance parameters, the volatility path, and standardised residuals, plus :meth:forecast for multi-step variance.

Examples:

>>> import statspai as sp
>>> import numpy as np
>>> rng = np.random.default_rng(0)
>>> T = 400
>>> eps = np.zeros(T)
>>> s2 = np.ones(T)
>>> omega, a1, b1 = 0.05, 0.1, 0.85
>>> for t in range(1, T):
...     s2[t] = omega + a1 * eps[t - 1] ** 2 + b1 * s2[t - 1]
...     eps[t] = np.sqrt(s2[t]) * rng.standard_normal()
>>> res = sp.garch(eps, p=1, q=1)
>>> bool(res.persistence < 1.0)   # alpha + beta < 1 => stationary
True
>>> res.forecast(horizon=3).shape
(3,)

persistence property

persistence: float

Sum of the coefficients that carry the variance forward: sum(alpha) + sum(beta), plus half the threshold terms for 'gjr' (a symmetric innovation is negative half the time); sum(beta) for 'egarch'.

forecast

forecast(horizon: int = 1) -> ndarray

Multi-step ahead variance forecast (analytic recursion).

E[eps^2_{T+h}] = sigma^2_{T+h|T} for h >= 1, so future squared shocks are replaced by their forecasts; every ARCH and GARCH lag is used.

forecast_mean

forecast_mean(horizon: int = 1) -> ndarray

Multi-step forecast of the series: the constant mean plus the autoregressive forecast of the disturbance (ar=).

value_at_risk

value_at_risk(alpha: float = 0.01) -> float

One-step-ahead value at risk: the alpha quantile of the forecast distribution of the next observation.

mu_{T+1} + sigma_{T+1} * q_alpha with q_alpha the quantile of the standardised innovation (standard normal, or Student t scaled to unit variance). A return below it has probability alpha under the fitted model.

ITSResult dataclass

Bases: ResultProtocolMixin

Result container for :func:its.

Examples:

>>> import statspai as sp
>>> import pandas as pd
>>> rng = np.random.default_rng(0)
>>> t = np.arange(60)
>>> y = 1 + 0.1 * t + (t >= 30) * 3.0 + rng.normal(scale=0.5, size=60)
>>> df = pd.DataFrame({"y": y, "t": t})
>>> res = sp.its(df, y="y", time="t", intervention=30)
>>> type(res).__name__
'ITSResult'
>>> float(res.level_change) > 0
True

to_dict

to_dict() -> dict

JSON-safe dict of every field (agent-native serialization).

JohansenLRTest dataclass

Bases: ResultProtocolMixin

Result of :func:johansen_lrtest.

Attributes:

Name	Type	Description
statistic, df, pvalue	(float, int, float)	Likelihood-ratio statistic and its chi-squared reference.
hypothesis	str	'beta', 'beta_known' or 'loading'.
rank	int	Cointegration rank held fixed under both hypotheses.
eigenvalues, eigenvalues_restricted	ndarray	Squared canonical correlations without and with the restriction.
beta, loading	DataFrame	Restricted estimates (columns _ce1 ...); beta is normalised so the first non-zero element of each column is one when that is possible.
n_used	int	Observations entering the statistic.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> trend = np.cumsum(rng.normal(size=300))
>>> df = pd.DataFrame({
...     "c": trend + rng.normal(size=300),
...     "y": trend + rng.normal(size=300),
... })
>>> test = sp.johansen_lrtest(df, rank=1, beta_known=[1, -1])
>>> type(test).__name__
'JohansenLRTest'
>>> list(test.beta.index)
['c', 'y']

LocalProjectionsResult dataclass

Bases: ResultProtocolMixin

Impulse-response container returned by :func:local_projections.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(3)
>>> T = 200
>>> shock = rng.normal(size=T)
>>> y = np.zeros(T)
>>> for t in range(1, T):
...     y[t] = 0.5 * y[t - 1] + 0.8 * shock[t] + rng.normal(0, 0.5)
>>> df = pd.DataFrame({"y": y, "shock": shock})
>>> res = sp.local_projections(df, outcome="y", shock="shock", horizons=8)
>>> type(res).__name__
'LocalProjectionsResult'
>>> len(res.horizons) == 9
True
>>> bool(res.irf[0] > 0)
True

LongRunVariance dataclass

Bases: ResultProtocolMixin

Long-run variance returned by :func:statspai.lrvar.

Attributes:

Name	Type	Description
lrvar	float or ndarray	The long-run variance; a K x K matrix for K > 1 series.
var_mean	float or ndarray	Variance of the sample mean, lrvar / n_obs.
se_mean	float or ndarray	Its square root (of the diagonal, for several series).
variance	float or ndarray	The ordinary variance gamma(0) with divisor n_obs.
kernel	str
bandwidth	float	Scale S of the kernel: lag j has weight k(j / S).
bandwidth_rule	str	'fixed' or the automatic rule that produced bandwidth.
prewhite	int
adjust	bool
n_obs	int
names	list of str

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> x = np.array([1.0, 2.0, 4.0, 3.0, 5.0, 4.0, 6.0, 5.0])
>>> fit = sp.lrvar(x, kernel="bartlett", bandwidth=2)
>>> round(float(fit), 4)
3.3047

to_frame

to_frame() -> DataFrame

The long-run covariance matrix as a labelled DataFrame.

summary

summary() -> str

Plain-text report of the estimate and how it was computed.

MarkovSwitchingResult dataclass

Bases: ResultProtocolMixin

Result of :func:mswitch.

Attributes:

Name	Type	Description
params	DataFrame	One row per coefficient: state ('all' for a coefficient common to the states), term, estimate, se, z, pvalue, ci_lower, ci_upper. The interval of a standard deviation is built on the log scale.
transition	DataFrame	transition.loc[i, j] is the probability of moving from state i to state j; rows sum to one.
transition_table	DataFrame	The same probabilities with delta-method standard errors and intervals built on the logit scale.
durations	DataFrame	Expected duration 1 / (1 - p_ii) of each state.
filtered, smoothed, predicted	DataFrame	Probability of each state given the data up to the date, the whole sample (Kim smoother), and the previous date.
fitted, resid	Series	One-step predictions E[y_t | y_{t-1}, ...] (the state means weighted by the predicted probabilities) and their errors.
theta, vcov	(Series, DataFrame)	Estimates and covariance on the estimation scale (log standard deviations, transition logits q as in Stata).
loglik, aic, bic, hqic, n_obs, n_params, states, model, ar, vce
converged	bool
starts	DataFrame	Log likelihood reached from each starting value; degenerate marks an end point where one state is never visited.
model_info	dict	order_rule, n_distinct_maxima, ergodic and notes.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(1)
>>> s = np.zeros(300, dtype=int)
>>> for t in range(1, 300):
...     stay = 0.95 if s[t - 1] == 0 else 0.9
...     s[t] = s[t - 1] if rng.uniform() < stay else 1 - s[t - 1]
>>> y = np.where(s == 0, -1.0, 2.0) + rng.normal(size=300)
>>> fit = sp.mswitch(y, states=2, seed=0)
>>> type(fit).__name__
'MarkovSwitchingResult'
>>> list(fit.smoothed.columns)
['state1', 'state2']

plot

plot(which: str = 'smoothed') -> Any

Probability of each state over time ('smoothed'/'filtered').

forecast

forecast(steps: int = 1) -> DataFrame

Forecasts of y and of the state probabilities.

Exact conditional expectations given the sample and the estimates; available without exogenous regressors, for model='dr' with ar=0 and for model='ar' with common AR coefficients.

MarkovSwitchingLRTest dataclass

Bases: ResultProtocolMixin

Result of :func:mswitch_lrtest.

Attributes:

Name	Type	Description
statistic	float	2 * (loglik_alt - loglik_null), not below zero.
pvalue	float	(1 + #{LR* >= LR}) / (n_valid + 1).
critical_values	Series	Bootstrap critical values at 10%, 5% and 1%: the ceil((n_valid + 1)(1 - level))-th smallest LR*; inf when there are too few replicates for the level.
reject	bool	pvalue <= alpha.
loglik_null, loglik_alt	float
states, null_states, n_obs, reps, n_valid	int
naive_df, naive_pvalue	(int, float)	Difference in the number of parameters and the chi-squared tail probability with it. Not a valid reference; for comparison only (see the Notes of :func:mswitch_lrtest).
replicates	DataFrame	One row per simulated series: lr, both log likelihoods, and the flags that diagnostics counts.
diagnostics	dict	n_failed (no statistic; left out of the p-value), n_floored (alternative fit below the null fit; statistic set to zero), n_multiple_maxima (the starts of the alternative ended at different local maxima), n_alt_not_converged, n_variance_floor, the same flags for the observed series, and the layout of the starting values.
null_params	Series	Estimates of the null model on the estimation scale; the series are simulated from them.
fit_null	MarkovSwitchingResult or None	None for a single regime (see null_params).
fit_alt	MarkovSwitchingResult
interpretation	str	What the test does and does not establish.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(1)
>>> s = np.zeros(120, dtype=int)
>>> for t in range(1, 120):
...     s[t] = s[t - 1] if rng.uniform() < 0.95 else 1 - s[t - 1]
>>> y = np.where(s == 0, -1.5, 1.5) + rng.normal(size=120)
>>> test = sp.mswitch_lrtest(y, states=2, reps=19, starts=3, seed=0)
>>> type(test).__name__
'MarkovSwitchingLRTest'
>>> bool(test.reject)
True

Hierarchy dataclass

An aggregation structure returned by :func:statspai.hierarchy.

Attributes:

Name	Type	Description
Y	DataFrame	Every series of the structure, aggregates first and the bottom level last; rows are time periods, columns series identifiers (the key values joined by sep).
S	DataFrame	Summing matrix: one row per series, one column per bottom-level series; Y.T = S @ Y_bottom.T.
tags	dict	Level name (the key columns joined by sep) to the list of its series.

Examples:

>>> import pandas as pd
>>> import statspai as sp
>>> df = pd.DataFrame({
...     "country": "AU", "state": ["A", "A", "B", "B"],
...     "t": [1, 2, 1, 2], "y": [1.0, 2.0, 3.0, 4.0]})
>>> h = sp.hierarchy(df, [["country"], ["country", "state"]], time="t", value="y")
>>> h.Y.columns.tolist()
['AU', 'AU/A', 'AU/B']
>>> h.S.values.tolist()
[[1.0, 1.0], [1.0, 0.0], [0.0, 1.0]]

bottom property

bottom: List[str]

Identifiers of the bottom-level series.

ReconcileResult dataclass

Bases: ResultProtocolMixin

Reconciled forecasts returned by :func:statspai.reconcile.

Attributes:

Name	Type	Description
forecasts	DataFrame	Coherent forecasts, laid out like the base forecasts.
method	str
G	DataFrame	The matrix that maps base forecasts of all series to bottom-level forecasts; forecasts = base @ (S @ G).T.
weights	ndarray or None	The matrix W of the generalised least squares projection.
shrinkage	float or None	Shrinkage intensity of mint_shrink.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> S = pd.DataFrame([[1, 1], [1, 0], [0, 1]], index=["T", "A", "B"],
...                  columns=["A", "B"], dtype=float)
>>> base = pd.DataFrame([[10.0, 4.0, 5.0]], columns=["T", "A", "B"])
>>> rec = sp.reconcile(base, S, method="ols")
>>> rec.forecasts.round(4).values.tolist()
[[9.6667, 4.3333, 5.3333]]

intervals

intervals(level: float = 95) -> Tuple[DataFrame, DataFrame]

Lower and upper bounds of the reconciled normal prediction intervals, laid out like forecasts.

Available when sd= was given to :func:statspai.reconcile. level is the coverage in percent.

SimpleForecastResult dataclass

Bases: ResultProtocolMixin

A fitted benchmark method returned by :func:statspai.simple_forecast.

Attributes:

Name	Type	Description
method	{'mean', 'naive', 'snaive', 'drift'}
period	int	Seasonal period of the seasonal naive method (1 otherwise).
fitted_values, residuals	ndarray	One-step fitted values and residuals; nan for the first observations a method cannot fit.
sigma	float	Residual standard deviation, sqrt(SSE / (T - K - M)) with K estimated parameters and M missing residuals.
drift	float or None	Average change per period (drift method).

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> y = np.array([10.0, 12.0, 11.0, 13.0, 14.0, 13.0, 15.0, 16.0])
>>> fit = sp.simple_forecast(y, "drift")
>>> round(fit.drift, 4)
0.8571
>>> fit.forecast(2, level=95).round(2)["forecast"].tolist()
[16.86, 17.71]

forecast

forecast(horizon: int = 10, level: Levels = (80, 95), *, bootstrap: bool = False, n_paths: int = 5000, seed: Optional[int] = 0) -> DataFrame

Point forecasts and prediction intervals.

Parameters:

Name	Type	Description	Default
horizon	int		10
level	float or sequence of float	Coverage in percent.	(80, 95)
bootstrap	bool	Take the intervals from future paths built by resampling the residuals instead of assuming them normal.	False
n_paths	int		5000
seed	int or None		0

Returns:

Type	Description
DataFrame	forecast, lower_<level>, upper_<level>.

simulate

simulate(horizon: int = 10, n_paths: int = 5, *, bootstrap: bool = True, seed: Optional[int] = 0) -> DataFrame

Future sample paths, one per column.

CumulativePeriodogramResult dataclass

Bases: ResultProtocolMixin

Bartlett's test returned by :func:statspai.cumulative_periodogram_test.

Attributes:

Name	Type	Description
statistic	float	sqrt(q) times the largest distance between the cumulative periodogram and the straight line of white noise.
pvalue	float
n	int
table	DataFrame	freq (cycles per observation), cumulative (share of the variance at frequencies up to freq) and expected (the same under white noise, 2 * freq).
alpha	float
band	float	Half-width of the 1 - alpha band around expected.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(3)
>>> res = sp.cumulative_periodogram_test(rng.normal(size=200))
>>> bool(res.pvalue > 0.05)
True

reject property

reject: bool

Whether white noise is rejected at level alpha.

summary

summary() -> str

Plain-text report of the test.

plot

plot(ax: Any = None) -> Any

Cumulative periodogram with the white-noise band.

SpectrumResult dataclass

Bases: ResultProtocolMixin

Estimated spectral density returned by :func:statspai.periodogram.

Attributes:

Name	Type	Description
table	DataFrame	One row per frequency: freq (cycles per observation), omega (radians per observation, 2 pi freq), cycle_length (observations per cycle, 1 / freq), spectrum and, for the periodogram estimates, lower and upper.
method	{'raw', 'smoothed', 'ar'}
scale	{'cycles', 'radians'}
df	float or None	Equivalent degrees of freedom of each ordinate (None for the autoregressive spectrum).
bandwidth	float or None	Equivalent bandwidth of the smoothing kernel, in cycles per observation (None for the autoregressive spectrum).
n	int	Observations used.
order	int or None	Autoregressive order (method='ar').
ar_coef	ndarray or None	Yule-Walker coefficients (method='ar').
var_pred	float or None	Innovation variance (method='ar').
settings	dict

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> res = sp.periodogram(rng.normal(size=120), method="smoothed", spans=5)
>>> list(res.table.columns)
['freq', 'omega', 'cycle_length', 'spectrum', 'lower', 'upper']

peak property

peak: Dict[str, float]

Frequency at which the estimate is largest.

summary

summary() -> str

Plain-text description of the estimate.

plot

plot(ax: Any = None, log: bool = True, omega: bool = False) -> Any

Plot the estimate against frequency.

Parameters:

Name	Type	Description	Default
ax	matplotlib Axes		None
log	bool	Logarithmic vertical axis, on which the confidence band has constant width.	True
omega	bool	Put radians rather than cycles per observation on the horizontal axis.	False

KalmanResult dataclass

Bases: ResultProtocolMixin

Output of :func:kalman_filter.

Attributes:

Name	Type	Description
predicted_state, predicted_cov	ndarray	E[X_t | Y_1..Y_{t-1}] and its mean squared error, shapes (T, m) and (T, m, m).
filtered_state, filtered_cov	ndarray	The same given Y_1..Y_t.
smoothed_state, smoothed_cov	ndarray or None	The same given the whole sample.
innovations	ndarray	One-step prediction errors Y_t - A_t - G_t X_{t|t-1}, (T, n), NaN where Y is missing.
innovations_cov	ndarray	Their covariance, (T, n, n), NaN in the rows and columns of missing elements.
std_innovations	ndarray	L_t^{-1} times the prediction errors, L_t the lower Cholesky factor of their covariance among the observed elements: uncorrelated with unit variance under the model.
loglik	float	Gaussian log-likelihood, the sum of loglik_obs from date burn + 1 on.
loglik_obs	ndarray	Contribution of each date; zero when the whole row is missing.
n_obs	int	Scalar observations that enter loglik.
n_dates	int	Dates with at least one of them.
x0, P0	ndarray	Moments of the initial state that were used.
init	str	'user', 'stationary', 'diffuse' or 'exact'.
predicted_cov_inf, filtered_cov_inf	ndarray or None	Under the exact diffuse initial state, the part of the covariance that is multiplied by infinity, (T, m, m); zero once the data have pinned the initial state down. predicted_cov and filtered_cov are then the finite part. None otherwise.
n_diffuse	int	Scalar observations absorbed by the diffuse initial state.
P0_inf	ndarray or None	Diffuse part of the initial covariance.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> y = np.cumsum(rng.normal(size=60)) + rng.normal(size=60)
>>> out = sp.kalman_filter(y, F=1.0, G=1.0, Q=1.0, R=1.0)
>>> isinstance(out, sp.KalmanResult)
True
>>> out.states("filtered").shape
(60, 2)

states

states(which: str = 'smoothed') -> DataFrame

State estimates with their standard errors as a DataFrame.

which is 'predicted', 'filtered' or 'smoothed'; the columns are the state names and <name>_se. Under the exact diffuse initial state a standard error is infinite as long as the data have not pinned that state down.

forecast

forecast(steps: int = 1, **future: Any) -> Dict[str, Any]

Forecast states and observations beyond the sample.

Parameters:

Name	Type	Description	Default
steps	int	Horizon.	1
**future	Any	Values of F, G, Q, R or A over the forecast dates, with a leading axis of length steps (or constant). A matrix that varied in the sample must be given here.	{}

Returns:

Type	Description
dict	state and obs (DataFrames of forecasts), state_se and obs_se (root mean squared errors), state_cov and obs_cov (arrays (steps, ., .)). The mean squared error of an observation includes R; the system matrices are taken as known.

plot

plot(which: str = 'smoothed', alpha: float = 0.05) -> Any

State paths with pointwise 1 - alpha bands.

StateSpaceResult dataclass

Bases: ResultProtocolMixin

Maximum-likelihood fit returned by :func:statespace.

Attributes:

Name	Type	Description
params	Series	Estimates of the parameter vector passed to build.
se	Series	Standard errors.
table	DataFrame	estimate, se, z, pvalue, lower, upper.
cov	DataFrame	Covariance of the estimates.
transformed	DataFrame or None	The same table for transform(params), by the delta method.
loglik, aic, bic	float
n_obs, n_dates	int	Scalar observations in the likelihood and dates carrying them; bic uses n_dates.
converged	bool	Scaled gradient near zero and a positive definite Hessian.
filter	KalmanResult	Filter and smoother output at the estimates.
model_info	dict	gradient, scaled_gradient, hessian_pd, steps (what each optimiser returned), vce.

Examples:

A local level model with both variances estimated:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(1)
>>> level = np.cumsum(0.3 * rng.normal(size=150))
>>> y = level + rng.normal(size=150)
>>> def build(th):
...     return {"F": 1.0, "G": 1.0, "Q": np.exp(th[0]), "R": np.exp(th[1])}
>>> fit = sp.statespace(y, build, [0.0, 0.0], init="diffuse", burn=1,
...                     param_names=["log_q", "log_r"])
>>> type(fit).__name__
'StateSpaceResult'
>>> list(fit.params.index)
['log_q', 'log_r']

forecast

forecast(steps: int = 1, **future: Any) -> Dict[str, Any]

Forecast at the estimates; see :meth:KalmanResult.forecast.

Parameter uncertainty is not included in the mean squared errors.

plot

plot(which: str = 'smoothed') -> Any

State paths with pointwise bands; see :meth:KalmanResult.plot.

DecompositionResult dataclass

Bases: ResultProtocolMixin

Decomposition of a series returned by :func:statspai.stl and :func:statspai.classical_decompose.

Attributes:

Name	Type	Description
observed, trend, remainder	ndarray
seasonal	ndarray	The seasonal component (the sum of them with several periods; the separate ones are in seasonal_components).
periods	tuple of int
multiplicative	bool	Whether observed = trend * seasonal * remainder rather than their sum (classical decomposition only).
method	str

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> t = np.arange(120)
>>> y = 0.1 * t + 3 * np.sin(2 * np.pi * t / 12) + np.cos(t)
>>> dec = sp.stl(y, period=12)
>>> dec.strength["seasonal"] > 0.9
True
>>> list(dec.to_frame().columns)
['observed', 'trend', 'seasonal', 'remainder', 'seasonally_adjusted']

seasonally_adjusted property

seasonally_adjusted: ndarray

The series with the seasonal component removed.

strength property

strength: Dict[str, float]

Strength of trend and of seasonality, each between 0 and 1: max(0, 1 - Var(R) / Var(T + R)) and max(0, 1 - Var(R) / Var(S + R)) (Wang, Smith and Hyndman, 2006). With several seasonal periods there is one seasonal_<period> entry per period as well.

to_frame

to_frame() -> DataFrame

Components as columns, on the index of the series.

forecast

forecast(horizon: int = 10, level: Levels = (80, 95), *, method: str = 'ets', **kwargs: Any) -> DataFrame

Forecast by modelling the seasonally adjusted series and adding back a seasonal naive forecast of the seasonal component.

Parameters:

Name	Type	Description	Default
horizon	int		10
level	float or sequence of float		(80, 95)
method	('ets', 'naive', 'drift', 'arima')	Model of the seasonally adjusted series: a non-seasonal ETS model chosen by AICc, the naive or drift method, or a non-seasonal ARIMA model chosen automatically. Further keyword arguments go to that model.	"ets"

Notes

The intervals are those of the seasonally adjusted series shifted by the seasonal forecast: the uncertainty of the seasonal component is ignored, as in R's forecast::stlf.

plot

plot(figsize: Tuple[float, float] = (9, 7)) -> Any

Stacked panels: the series and each component.

StructuralBreakResult

Bases: ResultProtocolMixin

Results from structural break tests.

Returned by :func:structural_break. Holds the detected break point(s), the sup-F statistic(s) with their Andrews (1993) asymptotic p-value(s), and segment goodness-of-fit (RSS, BIC).

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> # mean shift halfway through the sample
>>> y = np.concatenate([rng.normal(0.0, 1.0, 100),
...                     rng.normal(3.0, 1.0, 100)])
>>> df = pd.DataFrame({"y": y})
>>> res = sp.structural_break(df, y="y", method="sup-f")
>>> isinstance(res, sp.StructuralBreakResult)
True
>>> res.test_type
'Sup-F'
>>> res.n_obs
200
>>> res.n_breaks >= 1
True

plot

plot(ax: Any = None, **kwargs: Any) -> Any

Plot with break dates marked.

SVARResult

Bases: ResultProtocolMixin

Outcome of :func:svar.

Attributes:

Name	Type	Description
identification	str	"short-run", "long-run", "sign" or "proxy" (one shock: impact has one column, fevd shares are of the full forecast-error variance).
impact	DataFrame	The structural impact matrix P (rows: variables, columns: shocks): the response in the same period to a one-standard-deviation shock. With sign restrictions, the pointwise median over the admissible set (which is not itself an admissible matrix).
A, B	DataFrame or None	Short-run identification: the estimated matrices of A e = B u.
C	DataFrame or None	Long-run identification: the long-run impact matrix.
table	DataFrame or None	The free parameters with standard error, z and p-value.
log_likelihood	float or None
overid	dict or None	Over-identified models: the likelihood-ratio test of the over-identifying restrictions (statistic, df, pvalue).
n_accepted, n_tried	int or None	Sign restrictions: draws kept, and rotations tried to get them.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> T = 300
>>> u = rng.normal(size=(T, 2))
>>> e = u @ np.array([[1.0, 0.0], [0.5, 0.8]]).T
>>> y = np.zeros((T, 2))
>>> for t in range(1, T):
...     y[t] = 0.5 * y[t - 1] + e[t]
>>> fit = sp.var(pd.DataFrame(y, columns=["x", "y"]), lags=1)
>>> res = sp.svar(fit, B=[[np.nan, 0], [np.nan, np.nan]])
>>> res.identification
'short-run'
>>> res.impact.shape
(2, 2)
>>> res.irf(4).shape   # (periods + 1) x variables x shocks, long format
(20, 4)

irf

irf(periods: int = 20, cumulative: bool = False, *, ci: Optional[str] = None, reps: int = 500, seed: Optional[int] = None, boot: str = 'kilian', alpha: Optional[float] = None) -> DataFrame

Structural impulse responses, one row per shock, response and period.

Column irf is the response to a one-standard-deviation shock. With sign restrictions it is the pointwise median over the admissible set and lower / upper are its alpha / 2 and 1 - alpha / 2 quantiles: the spread of the identified set at the estimated reduced form, not a confidence interval.

Parameters:

Name	Type	Description	Default
periods	int		20
cumulative	bool	Accumulated responses (the response of the level when the variable enters the VAR in differences).	False
ci	(None, 'bootstrap')	Short- and long-run identification: add se, lower and upper from a residual bootstrap of the VAR, the identification repeated on every bootstrap sample. Each replicate's shocks are signed to agree with the estimate.	None
reps	int		500
seed	int		None
boot	('kilian', 'efron', 'hall')	The bias-corrected bootstrap, the plain percentile band, or the band reflected about the estimate (see :func:statspai.irf).	'kilian'
alpha	float	1 - alpha is the pointwise coverage; default the alpha given to :func:svar.	None

fevd

fevd(periods: int = 20, *, ci: Optional[str] = None, reps: int = 500, seed: Optional[int] = None, boot: str = 'kilian', alpha: Optional[float] = None) -> DataFrame

Forecast-error variance decomposition: the share of the period-step forecast-error variance of each response due to each shock (zero at period = 0, as in Stata's irf table fevd). With sign restrictions, the pointwise median over the admissible set, with its quantile band.

ci='bootstrap' adds se, lower and upper for a short- or long-run identification, by the bootstrap of :meth:irf (same reps, seed, boot and alpha). The percentile bands stay inside the unit interval; boot='hall' reflects them about the estimate and need not.

historical_decomposition

historical_decomposition() -> DataFrame

Historical decomposition of every variable into shock contributions.

Returns one row per period and variable: actual, baseline (the path from the initial values and deterministic terms with all shocks set to zero), one column per identified shock with its accumulated contribution sum_s Theta_s[:, j] u_{j, t-s} from the first residual on, and with a proxy identification other for the shocks it leaves unidentified. actual = baseline + sum of the shock columns (+ other) holds exactly in every row. Not defined for sign restrictions, whose median impact matrix is not itself an admissible model.

FilterResult dataclass

Bases: ResultProtocolMixin

Trend and cycle returned by :func:statspai.tsfilter.

Attributes:

Name	Type	Description
observed, trend, cycle	Series	On the index of the input. trend and cycle are NaN where the filter is not defined (see :func:statspai.tsfilter).
method	str
params	dict	The settings that define the filter. For 'hamilton' also the regression coefficients, constant first.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> y = np.cumsum(np.random.default_rng(0).normal(0.5, 1, 120))
>>> res = sp.tsfilter(y, method="hp")
>>> bool(np.allclose(res.trend + res.cycle, res.observed))
True

to_frame

to_frame() -> DataFrame

Observed series, trend and cycle as columns.

gain

gain(omega: Union[Sequence[float], ndarray, None] = None) -> Series

Gain of the cycle filter: by how much it scales a wave of frequency omega.

Parameters:

Name	Type	Description	Default
omega	array - like	Frequencies in radians per observation, between 0 and pi. Default: 201 equally spaced points on that interval.	None

Returns:

Type	Description
Series	Gain, indexed by omega. The squared gain multiplies the spectral density of the input.

Raises:

Type	Description
MethodIncompatibility	For method='cf' without sma_order, whose weights differ at every date, so that no single gain describes it, and for the one-sided Hodrick-Prescott filter.

Notes

'hp': 4 s (1 - cos w)**2 / (1 + 4 s (1 - cos w)**2) with s = smooth. 'bw': s (1 - cos w)**m / ((1 + cos w)**m + s (1 - cos w)**m) with m = order and s = tan(pi / high)**(-2 m). Both describe the filter away from the ends of the sample. 'bk' and 'cf' with sma_order: |b_0 + 2 sum_j b_j cos(j w)| from the weights used, exact at every date where the cycle is defined. 'hamilton': |1 - sum_j beta_j exp(-i w (h + j))| from the estimated coefficients.

summary

summary() -> str

Plain-text description of the decomposition.

plot

plot(axes: Any = None) -> Any

Series with its trend (top) and the cycle (bottom).

TVPVARResult dataclass

Bases: ResultProtocolMixin

Result of :func:tvp_var.

Attributes:

Name	Type	Description
coef_filtered, se_filtered	ndarray, shape (T, K, K p + 1)	Coefficients given the data up to each date, and their standard errors. Axis 1 is the equation, axis 2 the regressor in the order of terms (lags first, _cons last, as in :func:statspai.var). T counts the dates after the first lags.
coef_smoothed, se_smoothed	ndarray or None	The same given the whole sample (method='kalman' only).
sigma	DataFrame	Error covariance. 'kalman': constant, the covariance of the one-step forecast errors (see the notes of :func:tvp_var). 'forgetting': the value at the last date.
sigma_t	ndarray(T, K, K) or None	The path of the error covariance ('forgetting').
variances	DataFrame or None	'kalman': one row per equation, the error variance obs and the innovation variance of every coefficient.
loglik	float	'kalman': the sum of the equations' Kalman-filter log likelihoods (the equations are estimated separately). 'forgetting': the sum of the log one-step predictive densities.
index	Index	Dates of the T rows.
var_names, terms	list of str
lags, n_obs, method, alpha

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> y = np.zeros((160, 2))
>>> for t in range(1, 160):
...     a = 0.2 + 0.6 * t / 160
...     y[t, 0] = a * y[t - 1, 0] + rng.normal()
...     y[t, 1] = 0.3 * y[t - 1, 1] + rng.normal()
>>> df = pd.DataFrame(y, columns=["x", "z"])
>>> fit = sp.tvp_var(df, lags=1, method="forgetting", lam=0.98)
>>> fit.coef_filtered.shape
(159, 2, 3)
>>> fit.terms
['L1.x', 'L1.z', '_cons']
>>> list(fit.irf(at=-1, periods=4).columns)
['date', 'shock', 'response', 'period', 'irf']

coefficients

coefficients(kind: Optional[str] = None, alpha: Optional[float] = None) -> DataFrame

Coefficient paths, one row per date, equation and regressor.

Parameters:

Name	Type	Description	Default
kind	('smoothed', 'filtered')	Default: smoothed when the method has them, else filtered.	'smoothed'
alpha	float	1 - alpha is the pointwise coverage of lower / upper (normal quantiles); default the alpha of the fit.	None

Returns:

Type	Description
DataFrame	date, equation, term, coef, se, lower, upper.

irf

irf(at: Any = -1, periods: int = 10, orthogonal: bool = True, cumulative: bool = False, kind: Optional[str] = None) -> DataFrame

Impulse responses of the VAR frozen at one or several dates.

The coefficients and the error covariance of date at are held fixed and the responses are those of that constant VAR. It answers "how would a shock propagate if the economy stayed as it was at that date", and ignores that coefficients keep moving afterwards.

Parameters:

Name	Type	Description	Default
at	label, int or list of them	Dates: labels of index, or integers, which are always positions among the estimation dates (negative from the end). The default is the last date.	-1
periods	int		10
orthogonal	bool	Responses to one-standard-deviation shocks from the Cholesky factor of the error covariance in the order of the variables, as in :func:statspai.irf. False: unit innovations.	True
cumulative	bool		False
kind	('smoothed', 'filtered')		'smoothed'

Returns:

Type	Description
DataFrame	date, shock, response, period, irf.

stability

stability(kind: Optional[str] = None) -> DataFrame

Largest companion root of the frozen VAR at every date.

Returns:

Type	Description
DataFrame	Indexed by date: max_root (modulus) and explosive (max_root >= 1). attrs['n_explosive'] counts the flagged dates.

forecast

forecast(steps: int = 1, alpha: Optional[float] = None) -> DataFrame

Forecasts from the last filtered coefficients.

The coefficients stay at their last filtered value (the forecast of a random walk) and are treated as known: the standard errors come from the error covariance of the last date only, as in VARResult.forecast, and understate the uncertainty.

Returns:

Type	Description
DataFrame	Indexed by horizon 1 .. steps; for each variable <name>, <name>_se, <name>_lower, <name>_upper.

plot

plot(equation: Optional[str] = None, terms: Optional[Sequence[str]] = None, kind: Optional[str] = None) -> Any

Coefficient paths with pointwise 1 - alpha bands.

Parameters:

Name	Type	Description	Default
equation	str	One equation; default all of them (one column each).	None
terms	list of str	Regressors to show; default all.	None
kind	('smoothed', 'filtered')		'smoothed'

UnitRootResult

Bases: ResultProtocolMixin

Outcome of :func:unitroot.

Attributes:

Name	Type	Description
test	str	"ADF", "DF-GLS", "Phillips-Perron" or "KPSS".
statistic	float	The t statistic on the lagged level (ADF, DF-GLS), the adjusted Z(t) (Phillips-Perron), or the KPSS statistic.
pvalue	float or None	MacKinnon (1994) approximate p-value for ADF and Phillips-Perron. None for DF-GLS and KPSS, whose null distributions are only tabulated at a few points.
null	str	"unit root", or "stationarity" for KPSS -- there reject is evidence of a unit root.
z_rho	float or None	Phillips-Perron only: the adjusted n (rho - 1) statistic.
critical_values	dict	{"1%": ..., "5%": ..., "10%": ...} for this sample size (and, for DF-GLS, this lag order).
reject	bool	Whether the null is rejected at alpha.
lags	int	Lagged differences in the test regression; for Phillips-Perron and KPSS, the Newey-West truncation lag.
lag_selection	str	"fixed" or the rule that chose lags.
n_obs	int	Observations in the test regression.
rho, se	float	Coefficient on the lagged level and its standard error. The regression differs by test: for ADF and DF-GLS the dependent variable is the difference, dy[t] = rho * y[t-1] + ..., so rho is near 0 under a unit root; for Phillips-Perron it is the level, y[t] = rho * y[t-1] + ..., so rho is near 1 (the number Stata's pperron prints). The two are one apart and share the standard error; rho_level and rho_diff name them without ambiguity. nan for KPSS.
rho_level, rho_diff	float	The coefficient on y[t-1] when the dependent variable is y[t] (rho_level, 1 under a unit root) and when it is dy[t] (rho_diff = rho_level - 1, 0 under a unit root), whichever regression the test ran. nan for KPSS.
trend	str	"n", "c" or "ct".

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> walk = pd.DataFrame({"y": np.cumsum(rng.normal(size=200))})
>>> res = sp.unitroot(walk, "y", lags=1)
>>> type(res).__name__
'UnitRootResult'
>>> sorted(res.critical_values)
['1%', '10%', '5%']
>>> res.lags
1

VARResult

Bases: ResultProtocolMixin

Results from VAR estimation.

Returned by :func:var; carries coefficient tables, the residual covariance, information criteria, and helper methods .irf(), .granger_test() and .plot_irf().

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(3)
>>> T = 200
>>> x = np.zeros(T); y = np.zeros(T)
>>> for t in range(1, T):
...     x[t] = 0.5 * x[t - 1] + rng.normal()
...     y[t] = 0.3 * y[t - 1] + 0.4 * x[t - 1] + rng.normal()
>>> df = pd.DataFrame({"y": y, "x": x})
>>> vr = sp.var(df, variables=["y", "x"], lags=2)
>>> type(vr).__name__
'VARResult'
>>> vr.var_names
['y', 'x']
>>> vr.lags
2
>>> print(vr.summary().splitlines()[0])
Vector Autoregression (VAR)

References

lutkepohl2005new

irf

irf(periods: int = 20, impulse: Optional[str] = None, response: Optional[str] = None, orthogonal: bool = True, **kwargs: Any) -> Dict[str, Any]

Compute impulse response functions (see :func:irf).

fevd

fevd(periods: int = 20, *, ci: Optional[str] = None, alpha: float = 0.05, reps: int = 1000, seed: Optional[int] = None, boot: str = 'kilian') -> DataFrame

Forecast-error variance decomposition under the Cholesky ordering of the variables.

One row per shock, response and period: the share of the period-step forecast-error variance of response due to the orthogonalised innovation of shock (zero at period = 0, as in Stata's irf table fevd). The orthogonalisation uses the same residual covariance as :meth:irf.

Parameters:

Name	Type	Description	Default
periods	int		20
ci	(None, asymptotic, bootstrap)	Add se, lower and upper. 'asymptotic' is the delta method (Stata's irf table fevd, stderr), with the band fevd -/+ z se, which can leave the unit interval. 'bootstrap' is the residual bootstrap of :func:irf with percentile bands, which cannot.	None
alpha	float		0.05
reps	int		1000
seed	int		None
boot	(kilian, efron)	The bias-corrected bootstrap, or the plain percentile one (see :func:irf).	'kilian'

Notes

A share of exactly zero or one (the first-ordered variable one step ahead) is on the boundary: its asymptotic standard error is zero and no normal band applies to it.

granger_test

granger_test(caused: str, causing: str) -> Dict[str, Any]

Test Granger causality.

granger_table

granger_table() -> DataFrame

Granger-causality Wald tests for every equation (Stata vargranger).

lag_exclusion

lag_exclusion() -> DataFrame

Wald tests that a whole lag can be dropped (Stata varwle).

lm_test

lm_test(lags: int = 2) -> DataFrame

LM test for residual autocorrelation at each lag order up to lags (Stata varlmar).

stability

stability() -> DataFrame

Eigenvalues of the companion matrix (Stata varstable); attrs['stable'] is True when all lie inside the unit circle.

equation_table

equation_table() -> DataFrame

Root MSE, R-squared and Wald chi-squared of each equation.

forecast

forecast(steps: int = 1, alpha: float = 0.05) -> DataFrame

Dynamic forecasts steps periods past the sample, with standard errors that take the coefficients as known.

plot_irf

plot_irf(periods: int = 20, orthogonal: bool = True, **kwargs: Any) -> Any

Plot impulse response functions.

VECResult

Bases: ResultProtocolMixin

Fitted vector error-correction model.

Attributes:

Name	Type	Description
alpha	DataFrame	Adjustment coefficients, one row per equation and one column per cointegrating relation.
beta	DataFrame	Cointegrating vectors (columns _ce1 ...), normalised so that the first rank rows are the identity. Rows _cons / _trend hold the deterministic terms inside the relation.
beta_se	DataFrame	Standard errors of beta; NaN for the normalised entries and for a constant that is identified only through the equations.
coefs	dict of DataFrame	For each equation D_<name>: coef and se of the lagged cointegrating relation(s) (L._ce1), the lagged differences (LD.<name>, L2D.<name> ...) and the deterministic terms.
sigma_ml	DataFrame	Maximum-likelihood residual covariance.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> common = np.cumsum(rng.normal(size=300))
>>> df = pd.DataFrame({
...     "y": common + rng.normal(scale=0.5, size=300),
...     "x": 0.5 * common + rng.normal(scale=0.5, size=300),
... })
>>> fit = sp.vec(df, ["y", "x"], lags=1, rank=1)
>>> type(fit).__name__
'VECResult'
>>> list(fit.beta.index)
['y', 'x', '_cons']
>>> list(fit.coefs)
['D_y', 'D_x']
>>> fit.lm_test(2).shape
(2, 4)

equation_table

equation_table() -> DataFrame

Per equation: number of coefficients, root MSE (divisor T - d), R-squared about zero, and the Wald chi-squared that every coefficient of the equation is zero.

companion

companion() -> ndarray

Companion matrix of the VAR in levels the model implies.

irf

irf(periods: int = 20, impulse: Optional[str] = None, response: Optional[str] = None, orthogonal: bool = True) -> Dict[str, Any]

Impulse responses of the levels.

Because the system has K - r unit roots, a shock can move the levels for good: the responses settle at a constant that need not be zero, where those of a stationary VAR die out.

Parameters:

Name	Type	Description	Default
periods	int	Horizon; responses are returned for periods 0 .. periods.	20
impulse	str	Restrict to one shock or one responding variable.	None
response	str	Restrict to one shock or one responding variable.	None
orthogonal	bool	Responses to one-standard-deviation orthogonalised shocks, by the Cholesky factor of the maximum-likelihood residual covariance in the order of the variables (Stata irf create after vec). With False, responses to a unit forecast error in one equation.	True

Returns:

Type	Description
dict	'irf': {"<impulse> -> <response>": array}; 'periods': 0 .. periods. The layout of :func:statspai.irf.

Notes

No confidence bands are computed. The ordering matters for the orthogonalised responses exactly as in a VAR.

fevd

fevd(periods: int = 20) -> DataFrame

Forecast-error variance decomposition of the levels under the Cholesky ordering of the variables.

One row per shock, response and period: the share of the period-step forecast-error variance of response due to the orthogonalised innovation of shock (zero at period = 0, as in Stata's irf table fevd).

forecast

forecast(steps: int = 1, alpha: float = 0.05) -> DataFrame

Dynamic forecasts of the levels steps periods past the sample.

Each forecast feeds the next through the error-correction equations, deterministic terms included. The standard error is the root of sum_i Phi_i Sigma Phi_i' with the maximum-likelihood residual covariance; it treats the coefficients as known, so it understates the uncertainty in short samples, and it grows without bound with the horizon because the levels are integrated.

Returns:

Type	Description
DataFrame	Indexed by horizon 1 .. steps with, for each variable, the columns <name>, <name>_se, <name>_lower and <name>_upper.

stability

stability() -> DataFrame

Eigenvalues of the companion matrix (Stata vecstable).

A VECM with K series and rank r has K - r eigenvalues equal to one by construction (attrs['unit_moduli']); the model is stable when the others lie inside the unit circle (attrs['stable']).

lm_test

lm_test(lags: int = 2) -> DataFrame

LM test for residual autocorrelation at each lag order (Stata veclmar): (T - d - 0.5) log(|Sigma| / |Sigma_s|) from the equations re-estimated, at the fitted cointegrating relation, with the residuals lagged s periods added (zeros where the lag does not exist); d counts the coefficients of one augmented equation. Chi-squared with K^2 degrees of freedom.

CrossCorrelogram dataclass

Bases: ResultProtocolMixin

Cross-correlogram returned by :func:statspai.xcorr.

Attributes:

Name	Type	Description
table	DataFrame	Indexed by lag (-lags..lags) with columns xcorr, lower, upper and outside.
haugh	dict	Portmanteau test of no cross-correlation at any lag in the table (prewhiten='ar' only; empty otherwise): statistic, pvalue, statistic_adj, pvalue_adj, df.
n_obs	int	Observations the correlations are computed from.
prewhiten	str or None
ar_orders	dict	Order of the AR filter of each series.
ar_coefs	dict	Its coefficients.
residuals	DataFrame	The two series that were cross-correlated, about their means.
alpha	float
names	tuple of str

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(1)
>>> y = rng.normal(size=200)
>>> x = np.r_[0.0, 0.0, y[:-2]] + 0.3 * rng.normal(size=200)
>>> cc = sp.xcorr(x, y, lags=4)
>>> int(cc.table["xcorr"].idxmax())    # x follows y by two periods
2

summary

summary() -> str

Plain-text table of the cross-correlations and the test.

plot

plot(ax: Any = None) -> Any

Stem plot of the cross-correlations with the band.

ZivotAndrewsResult dataclass

Bases: ResultProtocolMixin

Result of :func:zivot_andrews.

Attributes:

Name	Type	Description
statistic	float	Minimum over the candidate break dates of the t statistic of the lagged level.
break_index	int	Number of observations in the first regime: the break dummies switch on at observation break_index + 1 (1-based).
break_label	object	Index label (or time value) of the last observation before the break.
critical_values	dict	Asymptotic 1%, 5% and 10% values.
lags	int	Lagged differences in the regression at the chosen break.
model	str
path	Series	The t statistic at every candidate break, indexed by break_index.
coefficients	DataFrame	The regression at the chosen break.
n_obs	int	Observations in the regression at the chosen break.
lag_method	str	'fixed', or the criterion with the scope of the choice, e.g. 'bic, once' or 'ttest, break'.
lags_path	Series or None	With lag_selection='break', the order chosen at every candidate break, indexed like path.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> y = np.cumsum(0.1 + rng.normal(size=150))   # random walk with drift
>>> res = sp.zivot_andrews(y, model="both")
>>> sorted(res.critical_values)
['1%', '10%', '5%']
>>> len(res.path) == 150 - 2 * 22 + 1
True

reject property

reject: Dict[str, bool]

Whether the unit root is rejected at each tabulated level.

plot

plot(ax: Any = None) -> Any

The t statistic over the candidate break dates.

bagged_forecast

bagged_forecast(y: Any, forecaster: Union[str, Callable[..., Any]] = 'ets', *, horizon: int = 10, n_boot: int = 100, period: int = 1, block_size: Optional[int] = None, boxcox: Any = 'auto', seed: Optional[int] = 0, data: Optional[DataFrame] = None, **kwargs: Any) -> BaggedForecastResult

Forecast every bootstrapped version of a series and average.

Parameters:

Name	Type	Description	Default
y	(array - like, Series or str)		required
forecaster	str or callable	"ets", "arima", "naive", "snaive", "drift", "mean" (extra keyword arguments are passed on), or a function f(series, horizon) as in :func:statspai.tscv.	"ets"
horizon	int		10
n_boot	int	Number of series forecast, the original included.	100
period	int	As in :func:statspai.bootstrap_series.	1
block_size	int	As in :func:statspai.bootstrap_series.	1
boxcox	int	As in :func:statspai.bootstrap_series.	1
seed	int	As in :func:statspai.bootstrap_series.	1
data	DataFrame		None

Returns:

Type	Description
BaggedForecastResult	forecast (average, minimum and maximum over the members) and members.

Notes

With forecaster="ets" and the model chosen separately for each series this is the bagged ETS of Bergmeir, Hyndman and Benitez (2016), R's forecast::baggedETS. Bagging averages over model choice and parameter estimates, the two sources of error a single fit ignores.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> t = np.arange(72)
>>> y = (40 + 0.5 * t) * (1 + 0.15 * np.sin(2 * np.pi * t / 4))
>>> y = y + rng.normal(0, 1.5, 72)
>>> bag = sp.bagged_forecast(y, "snaive", horizon=4, n_boot=10, period=4)
>>> list(bag.forecast.columns)
['forecast', 'ensemble_min', 'ensemble_max']

References

bergmeir2016bagging, hyndman2026fpppy

bootstrap_series

bootstrap_series(y: Any, n_boot: int = 100, *, period: int = 1, block_size: Optional[int] = None, boxcox: Any = 'auto', seed: Optional[int] = 0, data: Optional[DataFrame] = None) -> DataFrame

Bootstrapped versions of a series that keep its trend, seasonal pattern and the autocorrelation of what is left.

Parameters:

Name	Type	Description	Default
y	(array - like, Series or str)	The series, in time order; a column name when data is given.	required
n_boot	int	Number of series returned. The first is the original.	100
period	int	Seasonal period. With period > 1 the series is decomposed by STL with a periodic seasonal component; otherwise a loess trend is removed.	1
block_size	int	Length of the resampled blocks of the remainder. Default twice the period, or min(8, n // 2) for non-seasonal data.	None
boxcox	(float, 'auto' or None)	Box-Cox parameter applied before decomposing, so that the components add up. "auto" is Guerrero's choice restricted to [0, 1]; None decomposes the series as it is, which is also what happens when the series is not strictly positive.	"auto"
seed	int or None		0
data	DataFrame		None

Returns:

Type	Description
DataFrame	One column per series, boot_0 (the original) to boot_<n_boot - 1>, on the index of the series.

Notes

The procedure is that of Bergmeir, Hyndman and Benitez (2016) and of R's forecast::bld.mbb.bootstrap. The draws are random, so the series are not those R returns for the same seed.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> t = np.arange(80)
>>> y = (50 + t) * (1 + 0.2 * np.sin(2 * np.pi * t / 4)) + rng.normal(0, 2, 80)
>>> boot = sp.bootstrap_series(y, 10, period=4)
>>> boot.shape
(80, 10)
>>> bool(np.allclose(boot["boot_0"], y))
True

References

bergmeir2016bagging

chow_test

chow_test(data: DataFrame, y: str, x: Optional[Sequence[str]] = None, *, break_point: Any, time: Optional[str] = None, break_vars: Optional[Sequence[str]] = None, vce: str = 'nonrobust', alpha: float = 0.05) -> Dict[str, Any]

Chow test for a structural break at a known date.

Tests that the regression coefficients are the same before and after break_point. The model is fitted once with every coefficient in break_vars allowed to differ between the regimes, and the restrictions that the differences are zero are tested. With every coefficient free and classical errors this is Chow's statistic [(RSS - RSS_1 - RSS_2) / k] / [(RSS_1 + RSS_2) / (n - 2k)].

Parameters:

Name	Type	Description	Default
data	DataFrame	Observations in time order.	required
y	str	Dependent variable.	required
x	sequence of str	Regressors. A constant is always included; with x=None the test is for a shift in the mean.	None
break_point	int, label or list of these	Where the second regime starts. An int is the number of observations in the first regime, so break_point=50 puts rows 0..49 before the break. With time= it is a value of that column, and anything else is looked up in the index of data; in both cases the row carrying the value is the first of the new regime (Stata estat sbknown, break()). A list tests several breaks jointly.	required
time	str	Column that dates the observations. The rows are sorted by it before the regimes are formed.	None
break_vars	sequence of str	Coefficients allowed to change, as names from x plus "const". Default: all of them.	None
vce	('nonrobust', 'hc0', 'hc1')	Covariance estimator of the Wald statistic. With 'nonrobust' the F statistic is exact under normal errors with one variance in every regime; 'hc1' allows the variance to differ between regimes and observations.	'nonrobust'
alpha	float	Level used for reject.	0.05

Returns:

Type	Description
dict	statistic (F), df1, df2 and pvalue; chi2 (the Wald statistic df1 * F) and chi2_pvalue (Stata reports this pair); reject; break_point (sizes of the regimes that end at each break); n_obs; rss_restricted and rss_unrestricted; regimes (a DataFrame of the coefficients in each regime); vce and test.

Notes

A regime needs at least as many observations as it has free coefficients. Chow's second test, for a regime too short to fit, is not computed.

The F reference holds for a date fixed in advance. Picking the date with the largest statistic and then reading its p-value here rejects a stable relation far more often than alpha.

References

chow1960tests

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> n = 120
>>> x = rng.normal(size=n)
>>> shift = np.where(np.arange(n) < 60, 0.0, 1.5)
>>> df = pd.DataFrame({"y": 1 + 0.5 * x + shift + rng.normal(size=n), "x": x})
>>> res = sp.chow_test(df, "y", ["x"], break_point=60)
>>> res["df1"], res["df2"]
(2, 116)
>>> bool(res["reject"])
True

engle_granger

engle_granger(data: DataFrame, variables: Optional[List[str]] = None, lags: Optional[int] = None, trend: str = 'c', alpha: float = 0.05) -> CointegrationResult

Engle-Granger (1987) two-step cointegration test.

Step 1: OLS regression of y (first variable) on the remaining variables and the deterministic terms selected by trend. Step 2: augmented Dickey-Fuller regression on the step-1 residuals without deterministic terms,

.. math:: \Delta e_t = \rho e_{t-1} + \sum_{j=1}^{L} \gamma_j \Delta e_{t-j} + u_t ,

whose t-ratio on rho is the test statistic. This is the specification of Stata's egranger (Schaffer, SSC) and the one the MacKinnon critical values are tabulated for; the deterministic terms live in step 1 only.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
variables	list of str	Variables to test (first is dependent).	None
lags	int	Number of lagged differences L in the step-2 regression (Stata egranger, lags()). If None, uses the fixed rule int(4 * (n / 100) ** 0.25).	None
trend	(c, ct, ctt)	Deterministic terms in step 1: constant; constant + linear trend (egranger, trend); constant + linear + quadratic trend (egranger, qtrend). The trend is t = 0, 1, ....	'c'
alpha	(0.01, 0.05, 0.1)	Level at which rank (1 = cointegrated) is decided.	0.01

Returns:

Type

Description

CointegrationResult

test_stats is the ADF t-ratio; critical_values the (1%, 5%, 10%) MacKinnon (2010) response-surface critical values for N = len(variables) series and T = n - 1 (the egranger convention); eigenvectors holds the step-1 coefficients (constant, regressors, trend terms); pvalue is MacKinnon's (1994) asymptotic p-value, the one statsmodels' coint reports, or None with more than six series, which his surfaces do not cover. The decision in rank uses the finite-sample critical value, so near a boundary it can disagree with pvalue < alpha.

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> T = 200
>>> x = np.cumsum(rng.normal(size=T))      # random walk
>>> y = 2.0 * x + rng.normal(size=T)       # cointegrated with x
>>> df = pd.DataFrame({"y": y, "x": x})
>>> res = sp.engle_granger(df, variables=["y", "x"])
>>> res.test_type
'Engle-Granger'
>>> res.n_vars
2
>>> bool(res.rank in (0, 1))
True
>>> isinstance(res.summary(), str)
True

References

[@engle1987integration], [@mackinnon1994approximate]

johansen

johansen(data: DataFrame, variables: Optional[List[str]] = None, lags: int = 1, trend: str = 'c', test: str = 'trace', alpha: float = 0.05) -> CointegrationResult

Johansen (1991) cointegration test.

Tests for the cointegration rank using the trace or maximum eigenvalue test statistic.

Equivalent to Stata's vecrank and R's urca::ca.jo().

Parameters:

Name	Type	Description	Default
data	DataFrame		required
variables	list of str	Variables to test.	None
lags	int	Number of lagged differences in the VECM. This is Stata's vecrank, lags(p) minus one and ca.jo(K = lags + 1).	1
trend	str	Deterministic specification (Stata vecrank, trend() names in parentheses): 'n' (none); 'rc' (rconstant, constant restricted to the cointegrating space = ca.jo(ecdet="const")); 'c' (constant, unrestricted constant = ca.jo(ecdet="none")); 'rt' (rtrend, trend restricted to the cointegrating space, unrestricted constant = ca.jo(ecdet="trend")); 'ct' (trend, unrestricted constant and trend). Stata's long names are accepted as aliases.	'c'
test	(trace, maxeig)		'trace'
alpha	(0.05, 0.01)	Level of the Osterwald-Lenum critical values (the table Stata's vecrank uses) that decide rank.	0.05

Returns:

Type	Description
CointegrationResult	test_stats[r] tests H0: rank <= r (r = 0..k-1); eigenvalues are the k largest squared canonical correlations; eigenvectors are the cointegrating vectors (columns), normalised so the first element is 1 (ca.jo's @V convention; restricted cases carry the deterministic coefficient as the last row).

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> trend = np.cumsum(rng.normal(size=200))  # common stochastic trend
>>> df = pd.DataFrame({
...     'gdp': trend + rng.normal(scale=0.5, size=200),
...     'consumption': 0.8 * trend + rng.normal(scale=0.5, size=200),
...     'investment': 0.3 * trend + rng.normal(scale=0.5, size=200),
... })
>>> result = sp.johansen(
...     df, variables=['gdp', 'consumption', 'investment'], lags=2)
>>> type(result).__name__
'CointegrationResult'
>>> isinstance(result.summary(), str)
True

tscv

tscv(y: Any, forecaster: Forecaster = 'naive', *, horizon: int = 1, initial: Optional[int] = None, step: int = 1, window: Optional[int] = None, period: int = 1, data: Optional[DataFrame] = None, **kwargs: Any) -> TSCVResult

Time series cross-validation on a rolling forecast origin.

For each origin the forecaster sees only the observations up to it and forecasts the next horizon periods; the errors against what then happened estimate the out-of-sample accuracy at each horizon.

Parameters:

Name	Type	Description	Default
y	(array - like, Series or str)	The series, in time order; a column name when data is given.	required
forecaster	str or callable	A built-in method -- "naive", "snaive", "drift", "mean", "ets", "arima" (extra keyword arguments are passed to it; the four benchmark methods take boxcox and biasadj) -- or a function f(train, horizon). The function receives the training window as a pd.Series (with the index of y) and returns the horizon forecasts, a table with a forecast column, or a fitted model with a .forecast(horizon) method.	"naive"
horizon	int	Longest forecast horizon.	1
initial	int	Number of observations in the first training window. The default is the larger of 10 and two seasonal periods plus one, the least the built-in methods can be fitted on; choose it so that the first window is long enough for the model you evaluate.	None
step	int	Number of periods between successive origins.	1
window	int	Length of a rolling training window. None (default) keeps every past observation, an expanding window.	None
period	int	Seasonal period handed to the built-in methods.	1
data	DataFrame		None

Returns:

Type	Description
TSCVResult	errors (origins by horizons), forecasts, actuals, accuracy() by horizon and failures.

Notes

The error matrix has the layout of R's forecast::tsCV: the row of an origin holds y[origin + h] - forecast. Origins where the forecaster raises are reported in failures and warned about; they are not dropped silently.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(1)
>>> y = 50 + np.cumsum(rng.normal(0.3, 1.0, 80))
>>> cv = sp.tscv(y, "drift", horizon=4, initial=30)
>>> cv.accuracy()["RMSE"].is_monotonic_increasing
True
>>> ses = sp.tscv(y, lambda train, h: sp.ets(train, "ANN"), horizon=1, initial=30)
>>> ses.errors.shape
(50, 1)

References

hyndman2026fpppy

hierarchy

hierarchy(data: DataFrame, spec: Sequence[Sequence[str]], *, time: str, value: str, sep: str = '/') -> Hierarchy

Aggregate bottom-level series into a hierarchical or grouped structure and build its summing matrix.

Parameters:

Name	Type	Description	Default
data	DataFrame	Long table: one row per bottom-level series and period, with the key columns that define the aggregation.	required
spec	list of list of str	One list of key columns per level, e.g. [["country"], ["country", "state"], ["country", "state", "region"]]. The last entry is the bottom level and must identify the most disaggregated series; the others may nest (a hierarchy) or cross (a grouped structure, e.g. by state and by purpose).	required
time	str	Period and value columns.	required
value	str	Period and value columns.	required
sep	str	Separator of the key values in the series identifiers.	"/"

Returns:

Type	Description
Hierarchy	Y (all series, wide), S (summing matrix), tags.

Notes

The bottom level must be balanced: every bottom series observed in every period. Missing cells would silently bias every aggregate, so they are refused.

Examples:

>>> import pandas as pd
>>> import statspai as sp
>>> df = pd.DataFrame({
...     "all": "Total",
...     "state": ["A", "A", "A", "A", "B", "B", "B", "B"],
...     "kind": ["x", "x", "y", "y", "x", "x", "y", "y"],
...     "t": [1, 2] * 4,
...     "y": [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]})
>>> spec = [["all"], ["all", "state"], ["all", "kind"], ["all", "state", "kind"]]
>>> h = sp.hierarchy(df, spec, time="t", value="y")
>>> h.S.shape
(9, 4)
>>> h.Y["Total"].tolist()
[16.0, 20.0]

References

hyndman2011optimal, hyndman2026fpppy

cumulative_periodogram_test

cumulative_periodogram_test(data: ArrayLike, y: Optional[str] = None, *, alpha: float = 0.05) -> CumulativePeriodogramResult

Bartlett's cumulative-periodogram test that a series is white noise.

Parameters:

Name	Type	Description	Default
data	(DataFrame, Series or array)	The series, in time order. With a DataFrame, y names the column.	required
y	str	Column of data.	None
alpha	float	Level of the band reported in .band and drawn by .plot().	0.05

Returns:

Type	Description
CumulativePeriodogramResult	.statistic, .pvalue, .table, .summary(), .plot().

Raises:

Type	Description
DataInsufficient	Fewer than four observations or a constant series.
MethodIncompatibility	A gap in the series, or alpha outside (0, 1).

Notes

The periodogram of the demeaned series is taken at the frequencies k / n, k = 0, ..., q - 1 with q = floor(n / 2) + 1, and accumulated into F_k, the share of the total up to k / n. Under white noise F_k is close to 2 k / n. The statistic is B = sqrt(q) * max_k |F_k - 2 k / n| and is referred to the Kolmogorov distribution, P(B > b) = 2 sum_j (-1)**(j-1) exp(-2 j**2 b**2). This is the statistic of Stata wntestb; R cpgram draws the same curve. The p-value sums the series to machine precision; Stata stops once a term falls below about 1e-8, so its r(p) agrees to that accuracy only.

A series that was estimated (residuals of a fitted model) is closer to the line than white noise would be, and the test is then conservative.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(3)
>>> e = rng.normal(size=200)
>>> x = np.cumsum(e) * 0.2 + e
>>> sp.cumulative_periodogram_test(x).reject
True

periodogram

periodogram(data: ArrayLike, y: Optional[str] = None, *, method: str = 'raw', spans: Union[int, Sequence[int], None] = None, taper: float = 0.1, detrend: str = 'linear', pad: float = 0.0, fast: bool = False, order: Optional[int] = None, max_order: Optional[int] = None, n_freq: int = 500, scale: str = 'cycles', alpha: float = 0.05) -> SpectrumResult

Estimate the spectral density of a series.

Parameters:

Name	Type	Description	Default
data	(DataFrame, Series or array)	The series, in time order. With a DataFrame, y names the column.	required
y	str	Column of data.	None
method	('raw', 'smoothed', 'ar')	'raw' is the periodogram at the Fourier frequencies j / n, j = 1, ..., floor(n / 2). 'smoothed' averages neighbouring periodogram ordinates with the modified Daniell kernel(s) given by spans. 'ar' fits an autoregression by Yule-Walker and returns its spectral density on n_freq equally spaced frequencies from 0 to 1/2.	'raw'
spans	int or sequence of int	Widths of the modified Daniell smoothers, applied one after the other (method='smoothed'). A span s uses s // 2 ordinates on each side, the outermost with half weight.	None
taper	float	Share of the sample at each end that is tapered by a split cosine bell before the transform, between 0 and 0.5. Ignored by 'ar'.	0.1
detrend	('linear', 'mean', 'none')	What is removed before the transform: a least-squares line, the mean, or nothing. 'ar' always removes the mean only.	'linear'
pad	float	Zeros appended before the transform, as a share of the sample length; pad=1 doubles the number of frequencies.	0
fast	bool	Pad with zeros to the next length whose prime factors are 2, 3 and 5 (the default of R spec.pgram). The frequencies are then j / n_padded and no longer the Fourier frequencies of the sample.	False
order	int	Autoregressive order (method='ar'). Default: the order with the smallest AIC among 0, ..., max_order.	None
max_order	int	Largest order considered. Default min(n - 1, floor(10 * log10(n))).	None
n_freq	int	Number of frequencies for method='ar'.	500
scale	('cycles', 'radians')	'cycles': density per cycle per observation, integrating to the variance over (-1/2, 1/2]. 'radians': density per radian, the former divided by 2 pi.	'cycles'
alpha	float	lower and upper are a pointwise 1 - alpha interval.	0.05

Returns:

Type	Description
SpectrumResult	.table has one row per frequency; .df and .bandwidth describe the smoother; .summary() and .plot().

Raises:

Type	Description
MethodIncompatibility	Unknown option values, spans missing or given with the wrong method, a gap in the series.
DataInsufficient	Fewer than four observations, a constant series, or an autoregressive order the sample cannot support.

Notes

With taper=0 and scale='radians' the raw estimate at w_j = 2 pi j / n is the textbook periodogram I(w_j) = |sum_t x_t exp(-i w_j t)|**2 / (2 pi n); neither the mean nor a constant changes it at these frequencies, so detrend='mean' and detrend='none' agree. scale='cycles' returns |sum_t x_t exp(-2 pi i f t)|**2 / n.

The defaults taper=0.1 and detrend='linear' are those of R spec.pgram; fast=True completes the match when n has a prime factor above 5. The estimate is divided by the mean square of the taper, 1 - 5 taper / 4.

Each ordinate is treated as spectrum * chi2(df) / df with df = 2 n / (n_padded * c * sum_j K_j**2), where K is the kernel and c = (1 - 93 taper / 64) / (1 - 5 taper / 4)**2 allows for the taper; the interval is spectrum * df / chi2 at the two quantiles. The raw periodogram has df = 2 and is not consistent: the interval does not shrink as the sample grows. bandwidth is the standard deviation of the kernel, with a continuity correction of 1/12, in cycles per observation.

For 'ar' the density is s2 / |1 - sum_k a_k exp(-2 pi i f k)|**2 with s2 the Yule-Walker innovation variance times n / (n - order - 1), and the order minimises n log(s2_k) + 2 k (R ar and spec.ar). No interval is reported: the chi-squared approximation does not apply.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> t = np.arange(200)
>>> rng = np.random.default_rng(1)
>>> x = np.cos(2 * np.pi * t / 20) + 0.5 * rng.normal(size=200)
>>> res = sp.periodogram(x)
>>> float(res.peak["cycle_length"])
20.0
>>> smooth = sp.periodogram(x, method="smoothed", spans=(3, 5))
>>> round(smooth.df, 2)
9.36
>>> sp.periodogram(x, method="ar", order=2).order
2

References

[@brockwell1991time], [@neusser2016time]

kalman_filter

kalman_filter(y: Any, *, F: Any, G: Any, Q: Any, R: Any, A: Optional[Any] = None, x0: Optional[Any] = None, P0: Optional[Any] = None, init: str = 'auto', kappa: float = 10000000.0, diffuse: Optional[Any] = None, smooth: bool = True, burn: int = 0, data: Optional[DataFrame] = None, state_names: Optional[Sequence[str]] = None) -> KalmanResult

Kalman filter and smoother for a linear state space model.

The model is X_t = F_t X_{t-1} + V_t with Var(V_t) = Q_t and Y_t = A_t + G_t X_t + W_t with Var(W_t) = R_t. All system matrices are taken as known; :func:statespace estimates them.

Parameters:

Name	Type	Description	Default
y	array-like, Series, DataFrame, or column name(s) with ``data``	Observations, T rows and n columns, in time order. NaN is a missing value: a missing element drops its row of the observation equation at that date, and a date with nothing observed is a pure prediction step.	required
F	array-like, ``(m, m)`` or ``(T, m, m)``	Transition matrix. With a time axis, entry t moves the state from date t - 1 to date t.	required
G	array-like, ``(n, m)`` or ``(T, n, m)``	Loadings of the observables on the state. With one observable, (m,) and (T, m) are accepted, so a regression with time-varying coefficients passes its regressor matrix.	required
Q	array-like, ``(m, m)`` or ``(T, m, m)``	Covariance of the state disturbance. May be singular.	required
R	array-like, ``(n, n)`` or ``(T, n, n)``	Covariance of the measurement error. May be singular or zero.	required
A	array-like, ``(n,)`` or ``(T, n)``	Intercept of the observation equation. Default zero.	None
x0	array - like	Mean and covariance of the state one step before the first observation. x0 defaults to zero; P0 follows init.	None
P0	array - like	Mean and covariance of the state one step before the first observation. x0 defaults to zero; P0 follows init.	None
init	('auto', 'stationary', 'diffuse', 'exact')	Rule for P0 when it is not given. 'stationary' solves P0 = F P0 F' + Q and needs constant F and Q with every eigenvalue of F inside the unit circle. 'diffuse' sets P0 = kappa I, an approximation. 'auto' is stationary when F is stable and 'diffuse' otherwise. 'exact' is the exact diffuse filter and smoother: the states marked in diffuse (all of them by default) start with infinite variance and loglik is the diffuse log-likelihood; see Notes.	'auto'
kappa	float	Variance of the initial state under init='diffuse'.	1e7
diffuse	sequence of bool, or symmetric matrix	One flag per state: True gives that element of X_0 an infinite variance. Passing it selects the exact diffuse filter (init must be 'auto' or 'exact'). The other states start from the matching block of P0 when it is given and from their stationary distribution otherwise, which needs constant F and Q and those states not to depend on diffuse ones. A positive semi-definite matrix is P0_inf itself, for diffuse directions that are combinations of the states; the finite part is then all of P0 (zero when P0 is not given).	None
smooth	bool	Also run the fixed-interval smoother.	True
burn	int	Leading dates left out of loglik (the filter still uses them).	0
data	DataFrame	Source of the columns named by y.	None
state_names	sequence of str	Names of the states. Default x1, x2, ....	None

Returns:

Type	Description
KalmanResult	Predicted, filtered and smoothed states with their covariances, prediction errors, the log-likelihood, .states(), .forecast(), .summary() and .plot(). Under the exact diffuse filter also predicted_cov_inf, filtered_cov_inf, n_diffuse and P0_inf.

Raises:

Type	Description
MethodIncompatibility	System matrices of the wrong shape, an unstable F under init='stationary', or a prediction-error covariance that is not positive definite.
DataInsufficient	y is empty.

Notes

The log-likelihood is the Gaussian prediction-error decomposition, summed over every observed element of y from date burn + 1: date t contributes -0.5 [n_t log(2 pi) + log det S_t + v_t' S_t^{-1} v_t] with n_t the number of elements observed and S_t the covariance of their prediction errors. Missing elements contribute nothing.

init='diffuse' is the large-variance approximation: the first prediction errors have variance of order kappa and enter the likelihood with it, and state moments differ from their limit by a term of order 1 / kappa at every date (about 1e-6 for the default kappa on series of unit scale). Use burn to leave the first dates out of the likelihood.

init='exact' is that limit itself. The initial covariance is P0_* + kappa P0_inf with P0_inf = diag(diffuse) and kappa sent to infinity; the recursions carry the pair (P_*, P_inf), reported as predicted_cov and predicted_cov_inf (likewise filtered), until the data have used P_inf up, and are the ordinary ones afterwards. During those dates the elements of Y_t enter one at a time (after a rotation by the eigenvectors of R_t when it is not diagonal), so several observables, missing values and time-varying matrices need no special case. n_diffuse scalar observations are absorbed by the initial state. Each contributes -0.5 [log(2 pi) + log F_inf] to loglik, F_inf being the coefficient of kappa in its prediction-error variance, so that

loglik = lim [loglik(kappa) + 0.5 n_diffuse log(kappa)].

The normal constant is kept for every observation, as in Durbin and Koopman and in statsmodels; R KFAS leaves it out for the absorbed ones, which makes its logLik larger by 0.5 n_diffuse log(2 pi). The value also depends on which vector is called diffuse: here it is X_0, so the first predicted state has P_inf = F_1 P0_inf F_1'; software that puts P_inf = I on the first predicted state differs by the constant 0.5 log pdet(F_1 P0_inf F_1'), zero for random-walk states. Smoothed moments, and predicted and filtered moments from the end of the diffuse period, do not depend on that choice; predicted and filtered states inside it do, in the directions that are still unidentified, where .states() reports an infinite standard error. On the diffuse dates innovations and innovations_cov hold the joint prediction error and the finite part of its covariance, and std_innovations is NaN.

The covariance update uses the Joseph form, and the smoother the backward recursion that inverts only prediction-error covariances, so singular Q and R need no special treatment. It gives the same moments as the Rauch-Tung-Striebel recursion J_t = P_{t|t} F' P_{t+1|t}^{-1}, P_{t|T} = P_{t|t} + J_t (P_{t+1|T} - P_{t+1|t}) J_t'.

Examples:

A local level model:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(1)
>>> level = np.cumsum(0.5 * rng.normal(size=100))
>>> y = level + rng.normal(size=100)
>>> out = sp.kalman_filter(y, F=1.0, G=1.0, Q=0.25, R=1.0)
>>> out.smoothed_state.shape
(100, 1)

An AR(2) in companion form; the likelihood is the exact Gaussian one:

>>> F = np.array([[0.5, 0.3], [1.0, 0.0]])
>>> Q = np.diag([1.0, 0.0])
>>> out = sp.kalman_filter(y, F=F, G=[1.0, 0.0], Q=Q, R=0.0)
>>> out.init
'stationary'
>>> out.forecast(2)["obs"].shape
(2, 1)

A local linear trend with an exact diffuse initial state:

>>> out = sp.kalman_filter(y, F=[[1.0, 1.0], [0.0, 1.0]], G=[1.0, 0.0],
...                        Q=np.diag([0.25, 0.01]), R=1.0, init="exact")
>>> out.n_diffuse
2

References

[@kalman1960new], [@koopman1997exact], [@neusser2016time]

classical_decompose

classical_decompose(y: Any, period: int, *, model: str = 'additive', data: Optional[DataFrame] = None) -> DecompositionResult

Classical decomposition by moving averages.

The trend is a centred moving average of one full cycle (a 2 x m average when the period is even), the seasonal figure the average detrended value at each position of the cycle. Kept as the textbook baseline; :func:statspai.stl is the better tool: it has no missing ends, lets the seasonal pattern change and can be made robust to outliers.

Parameters:

Name	Type	Description	Default
y	(array - like, Series or str)		required
period	int		required
model	('additive', 'multiplicative')	observed = trend + seasonal + remainder or their product.	"additive"
data	DataFrame		None

Returns:

Type	Description
DecompositionResult	trend and remainder are nan for the first and last half cycle, where the moving average is not defined.

Notes

Equal to R's stats::decompose and to statsmodels.tsa.seasonal.seasonal_decompose with its default two-sided filter.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> y = np.tile([8.0, 12.0, 10.0, 14.0], 6) * (1 + 0.02 * np.arange(24))
>>> dec = sp.classical_decompose(y, period=4, model="multiplicative")
>>> round(float(dec.seasonal[:4].mean()), 6)
1.0

References

hyndman2026fpppy

cusum_test

cusum_test(data: DataFrame, y: str, x: Optional[List[str]] = None, alpha: float = 0.05, method: str = 'recursive') -> Dict[str, Any]

CUSUM test for parameter stability.

Tests H0: parameters are stable vs H1: parameter shift.

Parameters:

Name	Type	Description	Default
data	DataFrame		required
y	str	Dependent variable.	required
x	list of str	Regressors.	None
alpha	float		0.05
method	(recursive, ols)	Which residuals are cumulated. 'recursive' is the test of Brown, Durbin and Evans on one-step-ahead prediction errors. 'ols' is the test of Ploberger and Kramer on the residuals of the full-sample fit, cumsum(e) / (s sqrt(n)) with s**2 = RSS / (n - k); its limit is a Brownian bridge, so the boundary is a constant (1.358 at 5%). This is strucchange::efp(type="OLS-CUSUM") and Stata's estat sbcusum, ols. statsmodels' breaks_cusumolsresid divides the residual sum of squares by n unless ddof=k is passed.	'recursive'

Returns:

Type

Description

dict

For method='ols': 'cusum' (the scaled partial sums, one per observation), 'statistic' and 'max_cusum' (their largest absolute value), 'p_value' (from the law of the supremum of a Brownian bridge), 'critical_value' (a constant), 'boundary_coef' (the same number), 'reject', 'n_obs' and 'method'. For method='recursive': Keys: 'cusum' (standardised CUSUM path of the recursive residuals, strucchange::efp(type = "Rec-CUSUM")$process without its leading 0), 'max_cusum' (its supremum in absolute value), 'statistic' (max_s |cusum_s| / (1 + 2 s / m), the statistic of strucchange::sctest and Stata estat sbcusum), 'p_value' (Brownian-motion crossing probability of that statistic), 'boundary_coef' (a solving crossing probability = alpha; 0.9479 at 5%), 'critical_value' (the Brown-Durbin-Evans crossing boundary, an array a * [1 + 2 s / (n - k)] that widens from a to 3a across the sample -- not a constant), 'reject' (True if the path crosses that boundary anywhere), and 'n_obs'.

Notes

This is the recursive-residual CUSUM of Brown, Durbin & Evans (1975), as in R's strucchange::efp(type="Rec-CUSUM"). Its boundary is linear in the recursion index with coefficient a (0.948 at the 5% level), not the constant sup|Brownian-bridge| value (1.358) that belongs to the OLS-CUSUM of Ploberger & Kramer (1992).

References

brown1975techniques, ploberger1992cusum

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(0)
>>> n = 120
>>> x = rng.normal(size=n)
>>> y = 1.0 + 0.5 * x + rng.normal(scale=0.5, size=n)  # stable relation
>>> df = pd.DataFrame({"y": y, "x": x})
>>> res = sp.cusum_test(df, y="y", x=["x"])
>>> sorted(res.keys())
['boundary_coef', 'critical_value', 'cusum', 'max_cusum', 'method',
 'n_obs', 'p_value', 'reject', 'statistic']
>>> ols = sp.cusum_test(df, y="y", x=["x"], method="ols")
>>> round(ols["critical_value"], 3)
1.358
>>> res["n_obs"]
120
>>> bool(res["reject"])
False

boxcox_lambda

boxcox_lambda(y: Any, period: int = 1, *, lower: float = -1.0, upper: float = 2.0, data: Optional[DataFrame] = None) -> float

Box-Cox parameter that makes the variation of a series the same across its level, by Guerrero's method.

The series is cut into consecutive cycles; lambda minimises the coefficient of variation, across cycles, of sd / mean^(1 - lambda).

Parameters:

Name	Type	Description	Default
y	(array - like, Series or str)	A strictly positive series.	required
period	int	Seasonal period; cycles of length 2 are used for non-seasonal data.	1
lower	float	Search interval.	-1 and 2
upper	float	Search interval.	-1 and 2
data	DataFrame		None

Returns:

Type	Description
float	Apply it with scipy.special.boxcox(y, lam) and invert with scipy.special.inv_boxcox. lambda = 0 is the logarithm, 1 no transformation.

Notes

The criterion is that of R's forecast::BoxCox.lambda and fabletools::guerrero. It is minimised to full precision, where R's optimize stops at a tolerance of about 1e-4; expect agreement to four decimals.

A back-transformed point forecast is the median of the forecast distribution, not its mean. For a logarithm the mean is exp(mu) * (1 + sigma_h^2 / 2) to second order, with sigma_h the forecast standard deviation on the log scale.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> t = np.arange(96)
>>> y = np.exp(0.03 * t) * (1 + 0.2 * np.sin(2 * np.pi * t / 4))
>>> abs(sp.boxcox_lambda(y, period=4)) < 0.05      # the logarithm
True

References

guerrero1993time

fourier_terms

fourier_terms(n: Any, period: float, K: int, *, start: int = 1) -> DataFrame

Sine and cosine regressors for a seasonal pattern.

K pairs sin(2 pi k t / period), cos(2 pi k t / period), k = 1..K, describe a smooth seasonal pattern with 2K coefficients instead of period - 1 dummies. They handle long and non-integer periods (52.18 weeks in a year) and several seasonal patterns at once, as regressors of :func:statspai.arima (dynamic harmonic regression) or of :func:statspai.regress.

Parameters:

Name	Type	Description	Default
n	(int, Series, DataFrame or Index)	Number of rows; or an object whose length and index are used.	required
period	float	Length of the seasonal cycle, in observations.	required
K	int	Number of harmonics, at most period / 2.	required
start	int	Time of the first row. The regressors of the h periods after a sample of length T are fourier_terms(h, period, K, start=T + 1).	1

Returns:

Type	Description
DataFrame	Columns sin1_<period>, cos1_<period>, ... A sine that is identically zero (2k = period) is left out.

Notes

The values equal R's forecast::fourier(x, K) and, for the future rows, fourier(x, K, h). Choose K by AICc: every extra pair costs two parameters.

Examples:

>>> import statspai as sp
>>> X = sp.fourier_terms(8, period=4, K=2)
>>> list(X.columns)
['sin1_4', 'cos1_4', 'cos2_4']
>>> sp.fourier_terms(2, period=4, K=1, start=9).round(6).values.tolist()
[[1.0, 0.0], [0.0, -1.0]]

References

hyndman2026fpppy

ljungbox

ljungbox(x: Any, lags: Union[int, Sequence[int], None] = None, *, model_df: Optional[int] = None, method: str = 'ljung-box', period: int = 1, data: Optional[DataFrame] = None) -> DataFrame

Ljung-Box (or Box-Pierce) portmanteau test of no autocorrelation.

Parameters:

Name	Type	Description	Default
x	array-like, pd.Series, str or fitted model	The series, usually residuals; a column name when data is given. A fitted sp.arima / sp.ets / sp.simple_forecast result can be passed directly: its innovation residuals are tested.	required
lags	int or sequence of int	Number of autocorrelations in the statistic. An integer gives one test at that lag, a sequence one test per entry. Default: 10 for non-seasonal data and twice the period for seasonal data, capped at a fifth of the sample.	None
model_df	int	Number of ARMA parameters estimated in the model the residuals come from; the reference distribution is chi-squared with lags - model_df degrees of freedom. Default 0, or p + q + P + Q when x is a fitted ARIMA model. A lag at or below model_df has no reference distribution: its p-value is missing and a warning says so.	None
method	('ljung-box', 'box-pierce')		"ljung-box"
period	int	Seasonal period, used only for the default number of lags.	1
data	DataFrame		None

Returns:

Type	Description
DataFrame	Indexed by lag: statistic, df, p_value.

Notes

Q = n (n + 2) sum_{k=1}^{L} r_k^2 / (n - k) (Ljung and Box, 1978) or Q = n sum r_k^2 (Box and Pierce, 1970). Equal to R's Box.test(x, lag, type, fitdf), Stata's wntestq and statsmodels' acorr_ljungbox(x, lags, model_df=).

A small p-value says the residuals are autocorrelated: the model leaves forecastable structure behind, and its prediction intervals are not to be trusted.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> e = rng.normal(size=200)
>>> out = sp.ljungbox(e, lags=10)
>>> out.index.tolist(), list(out.columns)
([10], ['statistic', 'df', 'p_value'])
>>> walk = np.cumsum(e)
>>> bool(sp.ljungbox(walk, lags=10)["p_value"].iloc[0] < 0.001)
True

References

ljung1978measure, box1970distribution

ndiffs

ndiffs(y: Any, *, alpha: float = 0.05, max_d: int = 2, data: Optional[DataFrame] = None) -> int

Number of first differences that make a series stationary.

The series is differenced until a KPSS test no longer rejects level stationarity at alpha, at most max_d times.

Parameters:

Name	Type	Description	Default
y	(array - like, Series or str)		required
alpha	(0.1, 0.05, 0.025, 0.01)	Level of each KPSS test (the levels its critical values are tabulated for).	0.10
max_d	int		2
data	DataFrame		None

Returns:

Type	Description
int

Notes

The rule of Hyndman and Khandakar (2008) and of R's forecast::ndiffs(test="kpss"), including its lag truncation trunc(3 sqrt(n) / 13). :func:statspai.arima with auto=True uses it to fix d before comparing models by AICc, because likelihoods of differently differenced series are not comparable. For a seasonal series take the seasonal differences first, see :func:statspai.nsdiffs.

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> e = rng.normal(size=300)
>>> sp.ndiffs(e), sp.ndiffs(np.cumsum(e))
(0, 1)

References

hyndman2008automatic, kwiatkowski1992testing

nsdiffs

nsdiffs(y: Any, period: int, *, threshold: float = 0.64, max_D: int = 1, data: Optional[DataFrame] = None) -> int

Number of seasonal differences a series needs.

A seasonal difference is taken while the strength of seasonality of an STL decomposition exceeds threshold.

Parameters:

Name	Type	Description	Default
y	(array - like, Series or str)		required
period	int	Seasonal period.	required
threshold	float	Strength of seasonality above which a seasonal difference is taken; 0.64 minimised forecast errors on the M3 competition data.	0.64
max_D	int		1
data	DataFrame		None

Returns:

Type	Description
int

Notes

The strength is max(0, 1 - Var(R) / Var(S + R)) from :func:statspai.stl with its defaults, the measure of Wang, Smith and Hyndman (2006); the rule is that of R's forecast::nsdiffs(test="seas").

Examples:

>>> import numpy as np
>>> import statspai as sp
>>> t = np.arange(120)
>>> rng = np.random.default_rng(0)
>>> y = 10 * np.sin(2 * np.pi * t / 12) + rng.normal(size=120)
>>> sp.nsdiffs(y, 12), sp.nsdiffs(rng.normal(size=120), 12)
(1, 0)

References

wang2006characteristic, hyndman2008automatic

seasonal_dummies

seasonal_dummies(n: Any, period: int, *, start: int = 1, drop_first: bool = True) -> DataFrame

Seasonal indicator regressors for a regression with a seasonal pattern.

Parameters:

Name	Type	Description	Default
n	(int, Series, DataFrame or Index)	Number of rows; or an object whose length and index are used.	required
period	int	Seasonal period.	required
start	int	Time of the first row; row t is in season (t - 1) mod period + 1. The rows of the h periods after a sample of length T are seasonal_dummies(h, period, start=T + 1).	1
drop_first	bool	Leave out the indicator of season 1, the reference season of a regression with an intercept (the "dummy variable trap").	True

Returns:

Type	Description
DataFrame	Columns season_2 ... season_<period> of zeros and ones.

Notes

R's forecast::seasonaldummy and the season term of tslm build the same columns (with the last season as the reference in seasonaldummy). With a long period prefer :func:statspai.fourier_terms.

Examples:

>>> import statspai as sp
>>> sp.seasonal_dummies(5, 4).values.tolist()
[[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1], [0, 0, 0]]

References

hyndman2026fpppy

unitroot

unitroot(data: Union[DataFrame, Series, ndarray], y: Optional[str] = None, *, test: str = 'adf', trend: str = 'c', lags: Union[int, str] = 'aic', max_lags: Optional[int] = None, time: Optional[str] = None, alpha: float = 0.05) -> UnitRootResult

Unit-root test for one series: ADF, DF-GLS, Phillips-Perron or KPSS.

Equivalent to Stata's dfuller / dfgls / pperron / kpss, R's urca::ur.df / urca::ur.ers / urca::ur.pp / urca::ur.kpss and statsmodels.tsa.stattools.adfuller / kpss.

Parameters:

Name	Type	Description	Default
data	(DataFrame, Series or array)	The series, in time order (or pass time=). Leading and trailing missing values are trimmed; a gap inside the series is an error.	required
y	str	Column to test when data is a DataFrame.	None
test	('adf', 'dfgls', 'pp', 'kpss')	"dfgls" is the Elliott-Rothenberg-Stock test, more powerful when the largest root is near one. "pp" is the Phillips-Perron test: the Dickey-Fuller regression without lagged differences, corrected for serial correlation by a Newey-West long-run variance. "kpss" is the Kwiatkowski-Phillips-Schmidt-Shin test, whose null is stationarity (around a constant with trend="c", around a linear trend with "ct"): there reject=True is evidence of a unit root.	"adf"
trend	('c', 'ct', 'n')	Deterministic terms under the alternative: a constant, a constant and a linear trend, or nothing ("n", ADF only). Use "ct" for a series that visibly trends.	"c"
lags	int or {'aic', 'bic'}	Lagged differences in the test regression. An information criterion searches 0 .. max_lags on the common sample those regressions share, then the test is run on every observation the chosen order allows. Stata's dfuller default is lags=0. For "pp" and "kpss" this is the Newey-West truncation lag instead: an integer, or the default ("aic" as left unset, or "auto") for the rule floor(4 (T/100)^(2/9)) for Phillips-Perron (Stata's pperron), floor(4 (T/100)^(1/4)) for KPSS (the l4 rule of the original paper and of R's tseries::kpss.test).	"aic"
max_lags	int	Upper end of the search. Default: Schwert's (1989) rule floor(12 * (T / 100) ** 0.25).	None
time	str	Column to sort a DataFrame by first.	None
alpha	float	Level for result.reject. DF-GLS supports 0.01, 0.05 and 0.10; KPSS also 0.025.	0.05

Returns:

Type	Description
UnitRootResult

Raises:

Type	Description
MethodIncompatibility	Unknown test / trend / lag rule, trend="n" with DF-GLS, a gap or an infinite value in the series, a missing value in time, or a DF-GLS level other than 1%, 5%, 10%.
DataInsufficient	Too few observations for the requested lag order, a constant series, or a series the test regression fits exactly.

Notes

Rejecting means the data are unlikely under a unit root. Not rejecting is not evidence of a unit root: these tests have little power against a root of, say, 0.95 in a few decades of quarterly data.

The ADF critical values are MacKinnon's response surface, not the interpolated Fuller table Stata prints, so they differ in the second decimal; the statistic is identical.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> e = rng.normal(size=300)
>>> df = pd.DataFrame({"walk": np.cumsum(e)})
>>> df["ar"] = 0.0
>>> for t in range(1, 300):
...     df.loc[t, "ar"] = 0.5 * df.loc[t - 1, "ar"] + e[t]
>>> bool(sp.unitroot(df, "walk").reject)
False
>>> bool(sp.unitroot(df, "ar").reject)
True
>>> gls = sp.unitroot(df, "ar", test="dfgls", lags=1)
>>> gls.test
'DF-GLS'

References

[@dickey1979distribution], [@said1984testing], [@elliott1996efficient], [@mackinnon1994approximate], [@schwert1989tests], [@phillips1988testing], [@kwiatkowski1992testing]

granger_causality

granger_causality(var_result: Optional[VARResult] = None, data: Optional[DataFrame] = None, caused: Optional[str] = None, causing: Optional[str] = None, lags: Optional[int] = None, id: Optional[str] = None, time: Optional[str] = None) -> Dict[str, Any]

Granger causality test.

Tests whether causing variable Granger-causes caused variable.

Parameters:

Name	Type	Description	Default
var_result	VARResult	Pre-estimated VAR model.	None
data	DataFrame	Data (if var_result not provided).	None
caused	str	Variable being tested for causation.	None
causing	str or list of str	Variable(s) hypothesised to cause; a list tests all their lags jointly (Stata vargranger's ALL row when it lists every other variable).	None
lags	int	Number of lags (if fitting new VAR).	None
id	str	Unit and period columns of a panel in data: the Dumitrescu-Hurlin test (plm::pgrangertest, Stata xtgcause), one regression per unit with its own coefficients, averaging the unit Wald statistics. Returns W_bar, Z_bar, Z_tilde (the fixed-T version, whose two-sided p-value is p_value) with their p-values and the per-unit statistics in individual. lags defaults to 1.	None
time	str	Unit and period columns of a panel in data: the Dumitrescu-Hurlin test (plm::pgrangertest, Stata xtgcause), one regression per unit with its own coefficients, averaging the unit Wald statistics. Returns W_bar, Z_bar, Z_tilde (the fixed-T version, whose two-sided p-value is p_value) with their p-values and the per-unit statistics in individual. lags defaults to 1.	None

Returns:

Type

Description

dict

Keys: 'F_stat', 'p_value', 'df1', 'df2', 'chi2', 'chi2_p_value', 'caused', 'causing', 'reject'. chi2 is the Wald statistic W with its chi2(df1) p-value (Stata vargranger after var); F_stat = W / df1 is referred to F(df1, T - m) (Stata vargranger after var, small). The residual variance in W follows the fitted model's se_df: 'stata' (divisor T) or 'r' (divisor T - m, Stata var, small dfk). reject uses the F p-value at 5%.

Examples:

x Granger-causes y by construction, but not vice versa:

>>> import numpy as np
>>> import pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(5)
>>> T = 200
>>> x = np.zeros(T)
>>> y = np.zeros(T)
>>> for t in range(1, T):
...     x[t] = 0.5 * x[t - 1] + rng.normal()
...     y[t] = 0.3 * y[t - 1] + 0.4 * x[t - 1] + rng.normal()
>>> df = pd.DataFrame({"y": y, "x": x})
>>> vr = sp.var(df, variables=["y", "x"], lags=2)
>>> gc = sp.granger_causality(vr, caused="y", causing="x")
>>> print(gc["p_value"] < 0.05)  # x helps predict y
True

Or fit the VAR implicitly from data:

>>> gc2 = sp.granger_causality(data=df, caused="x", causing="y", lags=2)
>>> print(gc2["reject"])  # y does not Granger-cause x
False

References

[@granger1969investigating], [@dumitrescu2012testing]

irf

irf(var_result: VARResult, periods: int = 20, impulse: Optional[str] = None, response: Optional[str] = None, orthogonal: bool = True, cumulative: bool = False, sigma_df: Optional[str] = None, *, ci: Optional[str] = None, alpha: float = 0.05, reps: int = 1000, seed: Optional[int] = None, boot: str = 'kilian') -> Dict[str, Any]

Compute impulse response functions from VAR.

Simple responses are the MA coefficients Phi_s of the VAR; orthogonalised responses are Phi_s P with P the lower Cholesky factor of the residual covariance (variable order = Cholesky order).

Parameters:

Name	Type	Description	Default
var_result	VARResult	Estimated VAR model.	required
periods	int	Number of periods for IRF.	20
impulse	str	Impulse variable (if None, all).	None
response	str	Response variable (if None, all).	None
orthogonal	bool	Orthogonalized (Cholesky) IRF.	True
cumulative	bool	Return cumulative responses sum_{j<=s} (Stata cirf / coirf; vars::irf(cumulative = TRUE)).	False
sigma_df	(ml, unbiased)	Divisor of the residual covariance that the orthogonalisation uses: 'ml' = T (Stata irf create after var), 'unbiased' = T - m with m regressors per equation (vars::Psi, and Stata after var, dfk). Default follows the fitted model's se_df ('stata'/'ml' -> 'ml', 'r'/'unbiased' -> 'unbiased'). Irrelevant when orthogonal=False.	'ml'
ci	(None, asymptotic, bootstrap)	Add standard errors and a band. 'asymptotic': delta method on the joint normal limit of the coefficients and the residual covariance (Stata's irf create standard errors). 'bootstrap': residual bootstrap, the series regenerated recursively from its first lags observations and the VAR re-estimated (irf create, bs; vars::irf(boot = TRUE)).	None
alpha	float	The band has nominal pointwise coverage 1 - alpha.	0.05
reps	int	Bootstrap replications.	1000
seed	int	Seed of the bootstrap draws.	None
boot	(kilian, efron, hall)	Bootstrap band. 'kilian' is the bias-corrected bootstrap: it first estimates the small-sample bias of the lag coefficients from reps samples, then draws reps more from the bias-corrected VAR and corrects each of their estimates, with percentile bands. 'efron' is the plain percentile interval (the quantiles of the replicates, which is what Stata's irf create, bs and vars report). 'hall' reflects those quantiles about the estimate. The default costs twice the fits and departs from the two reference programs because their band undercovers: in a bivariate VAR(1) with a root of 0.92 and 60 observations, 90% bands for the own response at horizons 4 to 8 covered 38% ('efron'), 56% ('hall'), 62% (ci='asymptotic') and 86% ('kilian') of the time over 300 samples. With a largest root of 0.6 and 60 observations the figures were 77%, 70%, 78% and 89%; with a root of 0.2 the percentile and the bias-corrected band both covered 93%. The uncorrected two warn when the largest estimated root exceeds 0.9.	'kilian'

Returns:

Type	Description
dict	Keys: 'irf' (dict of arrays), 'periods'; with ci also 'se', 'lower', 'upper' (dicts with the same keys as 'irf') and 'ci' (the settings used). The bootstrap 'se' is the standard deviation of the replicates.

Notes

Both bands are pointwise and condition on the lag order. The asymptotic standard error of an orthogonalised response at horizon 0 comes from the residual covariance alone and is zero for the responses the Cholesky ordering sets to zero. Near a unit root the normal approximation is poor and the bootstrap replicates are biased towards zero persistence. boot='kilian' is built for that case; it removes the first-order bias and keeps every replicate stationary, and it is not a unit-root method either.

Examples:

>>> import statspai as sp
>>> import numpy as np, pandas as pd
>>> rng = np.random.default_rng(3)
>>> T = 200
>>> x = np.zeros(T); y = np.zeros(T)
>>> for t in range(1, T):
...     x[t] = 0.5 * x[t - 1] + rng.normal()
...     y[t] = 0.3 * y[t - 1] + 0.4 * x[t - 1] + rng.normal()
>>> df = pd.DataFrame({"y": y, "x": x})
>>> vr = sp.var(df, variables=["y", "x"], lags=2)
>>> out = sp.irf(vr, periods=10, impulse="x", response="y")
>>> sorted(out.keys())
['irf', 'periods']
>>> len(out["irf"]["x -> y"])  # s = 0, 1, ..., 10
11
>>> out["periods"][0]
0
>>> bool(np.isfinite(out["irf"]["x -> y"]).all())
True

References

[@lutkepohl2005new], [@runkle1987vector], [@kilian1998small], [@kilian2017structural]

varsoc

varsoc(data: DataFrame, variables: Optional[Sequence[str]] = None, *, maxlag: int = 4, trend: str = 'c', alpha: float = 0.05, exog: Optional[Sequence[str]] = None) -> DataFrame

Lag-order selection statistics for a VAR.

Fits VAR(0) .. VAR(maxlag) on the same observations (the first maxlag are set aside for every model, so the likelihoods are comparable) and reports, for each order, the log likelihood, the likelihood-ratio test against the order below, the final prediction error and the Akaike, Hannan-Quinn and Schwarz criteria.

Parameters:

Name	Type	Description	Default
data	DataFrame	The series, in time order.	required
variables	sequence of str	Endogenous variables. Default: every numeric column.	None
maxlag	int	Highest order considered.	4
trend	('c', 'ct', 'n')	Deterministic terms: constant, constant and trend, or none.	'c'
alpha	float	Level of the sequential likelihood-ratio test.	0.05

Returns:

Type	Description
DataFrame	Indexed by lag, with columns LL, LR, df, p, FPE, AIC, HQIC, SBIC. attrs['selected'] maps each criterion to the order it picks: the minimum of FPE / AIC / HQIC / SBIC, and for LR the highest order whose test against the order below rejects. attrs['n'] is the common number of observations.

Notes

The criteria are scaled by the number of observations, AIC = (-2 LL + 2 m) / T with m the total number of coefficients, as in Stata's varsoc and Lütkepohl's text; R's vars::VARselect reports log|Sigma| + 2 m / T, which ranks the orders the same way.

SBIC and HQIC are consistent for the true order; AIC and FPE overestimate it with positive probability in large samples but can forecast better in small ones. When they disagree, check the residuals of the smaller model (sp.estat(fit, 'varlmar')) before settling.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(1)
>>> T = 300
>>> y = np.zeros((T, 2))
>>> for t in range(1, T):
...     y[t] = 0.5 * y[t - 1] + rng.normal(size=2)
>>> table = sp.varsoc(pd.DataFrame(y, columns=["a", "b"]), maxlag=3)
>>> int(table.attrs["selected"]["SBIC"])
1

References

lutkepohl2005new

vec

vec(data: DataFrame, variables: Optional[Sequence[str]] = None, *, lags: int = 1, rank: int = 1, trend: str = 'c', alpha: float = 0.05) -> VECResult

Vector error-correction model (Johansen maximum likelihood).

Parameters:

Name	Type	Description	Default
data	DataFrame	The series in levels, in time order.	required
variables	sequence of str	The series of the system. Default: every numeric column. Their order matters for the normalisation: the first rank are the ones each cointegrating relation is solved for.	None
lags	int	Number of lagged differences, as in :func:statspai.johansen. Stata's vec, lags(p) counts lags of the levels: lags = p - 1.	1
rank	int	Number of cointegrating relations, 1 <= rank < K. Choose it with :func:statspai.johansen.	1
trend	('c', 'rc', 'n', 'rt', 'ct')	'c': unrestricted constant (the levels may drift); 'rc': constant inside the cointegrating relation only; 'n': no deterministic term; 'rt': a linear trend inside the relation and an unrestricted constant; 'ct': unrestricted constant and trend (a quadratic trend in the levels). Stata's names constant / rconstant / none / rtrend / trend are accepted.	'c'
alpha	float	One minus the confidence level (kept on the result for reporting).	0.05

Returns:

Type	Description
VECResult	alpha, beta (with beta_se), per-equation coefficient tables in coefs, the log likelihood and information criteria, and the methods equation_table(), stability(), lm_test() and summary().

Notes

With an unrestricted constant v the model writes v = alpha * mu + gamma with gamma orthogonal to alpha: mu is reported as the constant of the cointegrating relation (no standard error, it is not separately estimated) and gamma as the constant of each equation.

The cointegrating vector is super-consistent, but its normalisation is a choice: a relation that is insignificant in the first variable cannot be solved for it. If alpha of an equation is not different from zero, that series is weakly exogenous for the long-run parameters.

Examples:

>>> import numpy as np, pandas as pd
>>> import statspai as sp
>>> rng = np.random.default_rng(0)
>>> common = np.cumsum(rng.normal(size=400))
>>> df = pd.DataFrame({
...     "y": common + rng.normal(scale=0.5, size=400),
...     "x": 0.5 * common + rng.normal(scale=0.5, size=400),
... })
>>> fit = sp.vec(df, ["y", "x"], lags=1, rank=1)
>>> bool(abs(fit.beta.loc["x", "_ce1"] + 2.0) < 0.3)
True
>>> int(fit.stability().attrs["unit_moduli"])
1

References

johansen1991estimation; lutkepohl2005new