Time Series¶

Complete time series analysis: ACF/PACF, stationarity tests, ETS, ARIMA/SARIMA, automatic model selection, and decomposition.

Time series analysis module.

Phase 7A: ACF/PACF, differencing, and stationarity tests. Phase 7B: Exponential smoothing (ETS) models. Phase 7C: ARIMA models, forecasting, and automatic order selection. Phase 7D: Time series decomposition (classical and STL).

Public API:

acf(x) - Autocorrelation function (matches R stats::acf) pacf(x) - Partial autocorrelation function (matches R stats::pacf) diff(x) - Difference a time series (matches R base::diff) ndiffs(x) - Estimate differences for stationarity (matches R forecast::ndiffs) adf_test(x) - Augmented Dickey-Fuller unit root test (matches R tseries::adf.test) kpss_test(x) - KPSS stationarity test (matches R tseries::kpss.test) ets(y) - Fit an ETS state space model (matches R forecast::ets,

including “Z”-wildcard automatic selection; default model=”ZZZ”. Reported log-likelihood/AIC use the full- Gaussian convention — R’s concentrated numbers differ by a constant in n; model rankings/selection are identical. See timeseries._ets_fit and ._ets_select docstrings.)

forecast_ets(f) - Forecast from a fitted ETS model arima(y) - Fit an ARIMA model (R stats::arima numerics; interface subset documented on arima()) forecast_arima(f, y) - Forecast from a fitted ARIMA model auto_arima(y) - Automatic ARIMA order selection (matches R forecast::auto.arima) decompose(x) - Classical time series decomposition (matches R stats::decompose) stl(x) - STL decomposition (matches R stats::stl; identical parameters

reproduce R’s components to floating-point noise. Interface divergences, both deliberate: seasonal_window defaults to “periodic” where R requires it explicitly, and even/short loess spans raise instead of being silently rounded up.)

pystatistics.timeseries.acf(x, *, max_lag=None, conf_level=0.95, demean=True)[source]¶

Compute the autocorrelation function.

Matches R’s stats::acf(type=”correlation”).

Algorithm:: c(k) = (1/n) * sum_{t=1}^{n-k} (x_t - x_bar)(x_{t+k} - x_bar) acf(k) = c(k) / c(0)

Confidence intervals use Bartlett’s approximation under the white noise null hypothesis: +/- z_{alpha/2} / sqrt(n).

Parameters:

x (ArrayLike) – Time series (1-D array).
max_lag (int or None) – Maximum lag to compute. Default: min(10*log10(n), n-1), matching R.
conf_level (float) – Confidence level for CI bands. Must be in (0, 1).
demean (bool) – Whether to subtract the mean before computing. Default True, matches R.

Returns:

Result containing acf values, lags, and confidence bands.

Return type:

ACFSolution

Raises:

ValidationError – If inputs are invalid (NaN, non-1D, bad max_lag, etc.).

Notes

Validated against R stats::acf().

pystatistics.timeseries.pacf(x, *, max_lag=None, conf_level=0.95)[source]¶

Compute the partial autocorrelation function.

Matches R’s stats::pacf(). Uses the Durbin-Levinson recursion.

Algorithm:

phi_{1,1} = acf(1) For k = 2, 3, …:

phi_{k,k} = (acf(k) - sum_{j=1}^{k-1} phi_{k-1,j} * acf(k-j)) /
(1 - sum_{j=1}^{k-1} phi_{k-1,j} * acf(j))

phi_{k,j} = phi_{k-1,j} - phi_{k,k} * phi_{k-1,k-j} for j < k

Note: R’s pacf() does NOT include lag 0. Lags start from 1.

Confidence intervals: +/- z_{alpha/2} / sqrt(n) (same as ACF under white noise null hypothesis).

Parameters:

x (ArrayLike) – Time series (1-D array).
max_lag (int or None) – Maximum lag to compute. Default: min(10*log10(n), n-1), matching R.
conf_level (float) – Confidence level for CI bands. Must be in (0, 1).

Returns:

Result with kind=’partial’. Lags start at 1 (no lag 0).

Return type:

ACFSolution

Raises:

ValidationError – If inputs are invalid.

Notes

Validated against R stats::pacf().

pystatistics.timeseries.diff(x, differences=1, lag=1)[source]¶

Difference a time series.

Matches R’s base::diff().

For differences=1, lag=1: y[t] = x[t] - x[t-1] For differences=1, lag=12: y[t] = x[t] - x[t-12] (seasonal) For differences=2: apply differencing twice.

Parameters:

x (ArrayLike) – Time series (1-D array).
differences (int) – Number of times to difference. Default 1. Must be >= 1.
lag (int) – Lag for differencing. Default 1. Must be >= 1.

Returns:

Differenced series of length n - differences * lag.

Return type:

NDArray

Raises:

ValidationError – If inputs are invalid or the series is too short after differencing.

pystatistics.timeseries.ndiffs(x, *, test='kpss', significance_level=0.05, max_d=2)[source]¶

Estimate the number of differences needed for stationarity.

Matches R’s forecast::ndiffs(), including its default test. Repeatedly differences and tests for stationarity until either the test indicates stationarity or max_d is reached.

Parameters:

x (ArrayLike) – Time series.
test (str) – Stationarity test to use: ‘kpss’ (default, matching R forecast::ndiffs) or ‘adf’. The two tests have opposite null hypotheses (KPSS: stationary; ADF: unit root) and can recommend different d on borderline series.
significance_level (float) – Significance level for the test.
max_d (int) – Maximum number of differences to try.

Returns:

Recommended number of differences (0, 1, …, max_d).

Return type:

int

Raises:

ValidationError – If inputs are invalid.

pystatistics.timeseries.adf_test(x, *, n_lags=None, regression='ct')[source]¶

Augmented Dickey-Fuller test for unit root.

Matches R’s tseries::adf.test() (statistic and lag convention) with MacKinnon (1994) p-values (the surface used by statsmodels adfuller and by urca::ur.df’s critical values).

Tests H0: x has a unit root (non-stationary) vs H1: x is stationary

The test regression is:

Delta_x_t = alpha + beta*t + gamma*x_{t-1}

sum_{i=1}^{p} delta_i * Delta_x_{t-i} + eps_t

The test statistic is the t-statistic for gamma. The regression parameter controls which deterministic terms are included: - ‘nc’: no constant, no trend - ‘c’: constant only - ‘ct’: constant + linear trend (default). This is what R’s

tseries::adf.test always uses; with matching n_lags the statistic reproduces it exactly.

Parameters:

x (ArrayLike) – Time series (1-D array). Must have at least 3 observations.
n_lags (int or None) – Number of lagged difference terms. Default: floor((n-1)^(1/3)), matching R’s tseries::adf.test().
regression (str) – ‘nc’ (none), ‘c’ (constant), ‘ct’ (constant + trend, default).

Returns:

Test result with statistic, p-value, and critical values.

Return type:

StationaritySolution

Raises:

ValidationError – If inputs are invalid.

Notes

Statistic validated against both R tseries::adf.test() and statsmodels adfuller; p-values validated against statsmodels (MacKinnon 1994 surface, full range — both tails and the middle). tseries::adf.test itself interpolates a small table and caps its p-values at [0.01, 0.99]; inside that range the two agree, outside it this implementation keeps resolving while tseries saturates. Critical values use the MacKinnon (2010) finite-sample surface.

pystatistics.timeseries.kpss_test(x, *, regression='c', n_lags=None, lshort=True)[source]¶

KPSS test for stationarity.

Matches R’s tseries::kpss.test(), including its default bandwidth.

Tests H0: x is (level or trend) stationary vs H1: x has a unit root

NOTE: KPSS has the OPPOSITE null hypothesis to ADF.

Algorithm:

Regress x on deterministic terms (constant, or constant+trend).
Compute partial sums S_t = sum_{i=1}^{t} e_i of residuals.
Statistic: eta = (1/n^2) * sum S_t^2 / sigma^2_LR where sigma^2_LR is the long-run variance estimator using a Bartlett kernel: sigma^2_LR = gamma(0) + 2 * sum_{j=1}^{l} (1 - j/(l+1)) * gamma(j)

Parameters:

x (ArrayLike) – Time series (1-D array). Must have at least 3 observations.
regression (str) – ‘c’ (level stationarity, tseries null="Level") or ‘ct’ (trend stationarity, tseries null="Trend").
n_lags (int or None) – Number of lags for the Bartlett kernel. Default None uses the tseries::kpss.test rule selected by lshort: trunc(4*(n/100)^(1/4)) when lshort=True (tseries default) or trunc(12*(n/100)^(1/4)) when lshort=False. An explicit n_lags overrides lshort.
lshort (bool) – Bandwidth rule used when n_lags is None — the equivalent of tseries’s lshort argument. Default True.

Returns:

Test result with statistic, p-value, and critical values.

Return type:

StationaritySolution

Raises:

ValidationError – If inputs are invalid.

Notes

Validated against R tseries::kpss.test() for both null="Level" and null="Trend": at a matched bandwidth the statistic is reproduced exactly, and the p-value uses the same linear interpolation of the Kwiatkowski et al. (1992) table, clamped to tseries’s reporting range [0.01, 0.10].

class pystatistics.timeseries.ACFSolution(_result)[source]¶

Bases: SolutionReprMixin

Result from autocorrelation or partial autocorrelation computation.

Wraps a Result [ACFParams] envelope; every datum is exposed via a read-only @property so the public attribute surface is unchanged from the previous flat dataclass.

Parameters:: _result (Result[ACFParams])

acf¶

Autocorrelation values at each lag. For ACF, includes lag 0 (= 1.0). For PACF, starts at lag 1 (matching R’s pacf()).

Type:: NDArray

lags¶

Lag indices corresponding to each acf value.

Type:: NDArray

n_obs¶

Number of observations in the original series.

Type:: int

conf_level¶

Confidence level used for the confidence bands.

Type:: float

ci_upper¶

Upper confidence bound per lag.

Type:: NDArray

ci_lower¶

Lower confidence bound per lag.

Type:: NDArray

kind¶

‘correlation’ for ACF or ‘partial’ for PACF.

Type:: str

property acf: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property lags: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property n_obs: int¶

property conf_level: float¶

property ci_upper: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property ci_lower: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property kind: str¶

property info: dict¶

property timing: dict[str, float] | None¶

property backend_name: str¶

property warnings: tuple[str, ...]¶

summary()[source]¶

Return a human-readable summary of the ACF/PACF result.

Returns:: Multi-line summary string.
Return type:: str

class pystatistics.timeseries.StationaritySolution(_result)[source]¶

Bases: SolutionReprMixin

Result from a stationarity test (ADF, KPSS, PP).

Wraps a Result [StationarityParams] envelope; every datum is exposed via a read-only @property so the public attribute surface is unchanged from the previous flat dataclass.

Parameters:: _result (Result[StationarityParams])

statistic¶

Test statistic value.

Type:: float

p_value¶

p-value of the test. For KPSS, may be truncated to the table range (0.01 to 0.10).

Type:: float

method¶

Name of the test, e.g. ‘Augmented Dickey-Fuller’, ‘KPSS’.

Type:: str

alternative¶

Alternative hypothesis description, e.g. ‘stationary’ for ADF, ‘unit root’ for KPSS.

Type:: str

n_lags¶

Number of lags used in the test.

Type:: int

n_obs¶

Number of observations used (after differencing/lag adjustments).

Type:: int

critical_values¶

Critical values at standard significance levels, e.g. {‘1%’: -3.43, ‘5%’: -2.86, ‘10%’: -2.57}.

Type:: dict[str, float] | None

property statistic: float¶

property p_value: float¶

property method: str¶

property alternative: str¶

property n_lags: int¶

property n_obs: int¶

property critical_values: dict[str, float] | None¶

property info: dict¶

property timing: dict[str, float] | None¶

property backend_name: str¶

property warnings: tuple[str, ...]¶

summary()[source]¶

Return a human-readable summary of the stationarity test result.

Returns:: Multi-line summary string matching R-style output.
Return type:: str

pystatistics.timeseries.ets(y, *, model='ZZZ', period=1, damped=None, alpha=None, beta=None, gamma=None, phi=None, ic='aicc', tol=1e-10, max_iter=1000)[source]¶

Fit an ETS (ExponenTial Smoothing) state space model.

Matches R’s forecast::ets(): each component of the model string may be a concrete letter or a "Z" wildcard, and wildcards are resolved by fitting every admissible candidate and selecting the one minimising ic. The default model="ZZZ" performs full automatic selection, as in R. See the module docstring for the exact candidate table and the documented divergences from R.

Parameters:

y (ArrayLike) – Time series (1-D). Multiplicative components require strictly positive values.
model (str) – ETS model string — error, trend, season — e.g. 'ZZZ' (default: select everything), 'ANN', 'AAdN', 'MZZ', 'ZZN'.
period (int) – Seasonal period (e.g. 12 for monthly, 4 for quarterly). Seasonal models require 2 <= period <= 24.
damped (bool or None) – Damped-trend control. With a concrete model string it forces or forbids damping; with a 'Z' trend, None means both damped and undamped candidates are tried, True/False restricts the candidate set accordingly. damped=True with a trend fixed to 'N' raises (R: “Forbidden model combination”).
alpha (float or None) – Fix specific smoothing parameters (applies to every candidate). Fixed values must lie in R’s “usual” region (see _ets_fit.py); out-of-range values raise instead of being coerced into bounds.
beta (float or None) – Fix specific smoothing parameters (applies to every candidate). Fixed values must lie in R’s “usual” region (see _ets_fit.py); out-of-range values raise instead of being coerced into bounds.
gamma (float or None) – Fix specific smoothing parameters (applies to every candidate). Fixed values must lie in R’s “usual” region (see _ets_fit.py); out-of-range values raise instead of being coerced into bounds.
phi (float or None) – Fix specific smoothing parameters (applies to every candidate). Fixed values must lie in R’s “usual” region (see _ets_fit.py); out-of-range values raise instead of being coerced into bounds.
ic (str) – Selection criterion for wildcard models: 'aicc' (default, matching forecast::ets), 'aic', or 'bic'. AICc-based selection needs n >= 5 (see module docstring).
tol (float) – Convergence tolerance for the optimiser.
max_iter (int) – Maximum optimiser iterations.

Returns:

The fitted (for wildcards: selected) model. For wildcard requests, solution.info["selection"] records the requested string, the criterion, every candidate’s IC value, and any skipped candidates with the reason.

Return type:

ETSSolution

Raises:

ValidationError – On invalid inputs; on an explicitly requested component that cannot be honoured (e.g. multiplicative error on non-positive data, a seasonal letter with period of 1 or > 24, damped=True with trend 'N'); when no candidate is admissible; or when no candidate has a finite ic (series too short for AICc).
ConvergenceError – If candidates were attempted but none could be fitted.

class pystatistics.timeseries.ETSSolution(_result)[source]¶

Bases: SolutionReprMixin

Result from fitting an ETS model.

Wraps a Result [ETSParams] envelope; every datum is exposed via a read-only @property so the public attribute surface is unchanged from the previous flat dataclass.

Parameters:: _result (Result[ETSParams])

property spec: ETSSpec¶

property alpha: float¶

property beta: float | None¶

property gamma: float | None¶

property phi: float | None¶

property init_level: float¶

property init_trend: float | None¶

property init_season: ndarray[tuple[Any, ...], dtype[_ScalarT]] | None¶

property fitted_values: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property residuals: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property states: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property log_likelihood: float¶

Full Gaussian log-likelihood.

R’s forecast::ets reports the concentrated pseudo- log-likelihood -0.5*n*log(SSE); this value equals R’s plus the deterministic constant 0.5*n*[log(n/(2*pi)) - 1]. Model comparisons on the same data are identical either way.

property aic: float¶

AIC under the full-Gaussian log-likelihood.

Differs from forecast::ets’s printed AIC by the constant -n*[log(n/(2*pi)) - 1]; AIC differences and model rankings match R exactly (same parameter count k).

property aicc: float¶: AICc; same convention note as aic.

property bic: float¶: BIC; same convention note as aic.

property mse: float¶

property mae: float¶

property n_obs: int¶

property n_params: int¶

property converged: bool¶

property info: dict¶

property timing: dict[str, float] | None¶

property backend_name: str¶

property warnings: tuple[str, ...]¶

summary()[source]¶

Return a human-readable summary matching R’s forecast::ets() style.

Returns:: Multi-line summary.
Return type:: str

pystatistics.timeseries.forecast_ets(fitted, n_ahead=10, *, conf_level=(0.8, 0.95))[source]¶

Generate forecasts from a fitted ETS model.

Parameters:

fitted (ETSSolution) – A fitted ETS model (from ets()).
n_ahead (int) – Forecast horizon (number of steps ahead).
conf_level (float or sequence of float) – Prediction-interval confidence level(s) as fractions in (0, 1) (default (0.80, 0.95)). A single float requests one interval; a sequence requests several. Whole-percent values (e.g. 95) are rejected.

Returns:

Point forecasts and prediction intervals.

Return type:

ETSForecast

Raises:

ValidationError – If n_ahead < 1 or conf_level is invalid.

class pystatistics.timeseries.ETSForecast(mean, lower, upper, n_ahead, model, fitted)[source]¶

Bases: object

Forecast from a fitted ETS model.

Parameters:

mean (ndarray[tuple[Any, ...], dtype[_ScalarT]])
lower (dict[float, ndarray[tuple[Any, ...], dtype[_ScalarT]]])
upper (dict[float, ndarray[tuple[Any, ...], dtype[_ScalarT]]])
n_ahead (int)
model (ETSSpec)
fitted (ETSSolution)

mean¶

Point forecasts of length n_ahead.

Type:: NDArray

lower¶

Lower prediction-interval bounds keyed by confidence level as a fraction (e.g. {0.8: ..., 0.95: ...}).

Type:: dict[float, NDArray]

upper¶

Upper prediction-interval bounds keyed by confidence level.

Type:: dict[float, NDArray]

n_ahead¶

Forecast horizon.

Type:: int

model¶

Model specification used.

Type:: ETSSpec

fitted¶

The underlying fitted model.

Type:: ETSSolution

mean: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

lower: dict[float, ndarray[tuple[Any, ...], dtype[_ScalarT]]]¶

upper: dict[float, ndarray[tuple[Any, ...], dtype[_ScalarT]]]¶

n_ahead: int¶

model: ETSSpec¶

fitted: ETSSolution¶

summary()[source]¶

Return a human-readable forecast summary.

Returns:: Multi-line table of forecasts and intervals.
Return type:: str

class pystatistics.timeseries.ETSSpec(error, trend, season, period, damped)[source]¶

Bases: object

Specification of an ETS model type.

Parameters:

error (str)
trend (str)
season (str)
period (int)
damped (bool)

error¶

Error type: 'A' (additive) or 'M' (multiplicative).

Type:: str

trend¶

Trend type: 'N' (none), 'A' (additive), or 'Ad' (damped).

Type:: str

season¶

Season type: 'N' (none), 'A' (additive), or 'M' (multiplicative).

Type:: str

period¶

Seasonal period (1 when no seasonal component).

Type:: int

damped¶

True when trend is 'Ad'.

Type:: bool

error: str¶

trend: str¶

season: str¶

period: int¶

damped: bool¶

property name: str¶: Human-readable model name, e.g. 'ETS(A,Ad,N)'.

property n_states: int¶

1 (level) + trend? + season?.

Type:: Number of state variables

property n_params: int¶: Number of smoothing parameters (alpha, beta, gamma, phi).

pystatistics.timeseries.arima(y, *, order=(0, 0, 0), seasonal=None, include_mean=True, xreg=None, include_drift=False, fixed=None, method='css-ml', init=None, tol=1e-08, max_iter=1000, backend=None)[source]¶

Fit an ARIMA(p, d, q) or seasonal ARIMA(p, d, q)(P, D, Q)[m] model.

Matches the numerical behaviour of R’s stats::arima() — exact maximum likelihood via the same Kalman-filter approach, with results (log-likelihood, coefficients, information criteria, forecasts, standard errors) verified against R — and supports a documented subset of its interface (see R interface coverage below).

For non-seasonal ARIMA(p, d, q):

Difference the series d times.
Fit ARMA(p, q) to the differenced series.

For seasonal ARIMA(p, d, q)(P, D, Q)[m]:

Seasonally difference D times (lag m).
Difference d times.
Multiply out the seasonal and non-seasonal AR/MA polynomials to form effective ARMA coefficients.

- ``'css'``: Conditional sum of squares (fast, approximate).

- ``'ml'``: Exact maximum likelihood via the Kalman filter.

- ``'css-ml'``: CSS for initialization, then ML refinement (default).

Starting values:

AR: Yule-Walker estimates from autocorrelations.
MA: zeros.
Mean: sample mean of the differenced series.

Parameters:

y (ArrayLike) – Time series (1-D array).
order (tuple[int, int, int]) – (p, d, q) — AR order, differencing order, MA order.
seasonal (tuple[int, int, int, int] or None) – (P, D, Q, m) — seasonal AR order, seasonal differencing order, seasonal MA order, and period. None for non-seasonal models.
include_mean (bool) – Whether to include a mean term. Default True. Ignored (no mean is estimated) when the model has any differencing (d + D > 0), matching R stats::arima’s include.mean. With xreg / include_drift the mean is carried as the 'intercept' regression coefficient.
xreg (ArrayLike or None) – External regressors for regression with ARIMA errors: fit y = X @ beta + eta where eta follows the ARIMA process (R’s stats::arima(xreg=) / forecast::Arima(xreg=)). Shape (n,) or (n, k) with one row per observation of y. The regression coefficients appear on the solution as xreg_coef (named xreg1..xregk), with joint standard errors in vcov. Forecasting a model with xreg requires future regressor values (forecast_arima(..., new_xreg=)). Default None.
include_drift (bool) – Include a linear time-trend (drift) regressor — the models R reports “with drift” (forecast::Arima(include.drift=TRUE)). Reported as the 'drift' regression coefficient. Requires total differencing d + D <= 1 (the trend vanishes under higher-order differencing). Default False.
fixed (dict or ArrayLike or None) – Hold coefficients fixed during estimation (R’s stats::arima(fixed=) parameter masking). Primary form is a {name: value} mapping over the coefficient names (ar1, ma1, …, intercept, drift, xreg1, …) — e.g. fixed={'ma1': 0} holds ma1 at 0 and estimates the rest. A positional array aligned to the coefficient order with nan for free parameters (R’s convention) is also accepted. Fixed coefficients carry zero variance in vcov and do not count toward the information criteria. Default None.
method (str) – 'css', 'ml', or 'css-ml'. Default 'css-ml'.
init (ArrayLike or None) – Initial parameter values for the optimizer, in R coef() order: [ar_1..ar_p, ma_1..ma_q, sar_1..sar_P, sma_1..sma_Q, mean?] (the mean slot exists only when a mean is estimated, i.e. include_mean=True and d + D == 0). numpy.nan entries use the defaults (zero for coefficients, the sample mean of the differenced series for the mean). AR parts must be stationary (as in R); non-invertible MA parts are normalized to the invertible representative before optimization (R’s documented maInvert intent — R’s own implementation errors on such inits). Not supported with method='whittle'. Default None (internal starting values).
tol (float) – Convergence tolerance for the optimizer. Default 1e-8.
max_iter (int) – Maximum optimizer iterations. Default 1000.
backend (str | None)

Returns:

Fitted model with coefficients, residuals, and diagnostics.

Return type:

ARIMASolution

Raises:

ValidationError – If inputs are invalid.
ConvergenceError – If the optimizer fails to converge.

Notes

R interface coverage. Supported R parameters: order, seasonal, include.mean (include_mean), xreg (regression with ARIMA errors), include.drift (include_drift), fixed (parameter masking), method (‘css-ml’/’ml’/’css’), and init. Not exposed, by design: transform.pars, SSinit, kappa, n.cond, and optim.method/optim.control are knobs over R’s optimizer and state-space internals; pystatistics guarantees parity of RESULTS, not of internal knobs — AR stationarity and MA invertibility are handled internally (the latter with R’s own maInvert convention), the stationary state initialization is solved exactly, and the ML stage uses a better-of-two-starts strategy verified to reach equal-or-better optima than R on every reference model.

CSS convention. method='css' uses a zero-initialized conditional recursion over ALL observations, whereas R conditions on (and excludes) the first n.cond. Pure-CSS estimates can therefore differ slightly from R’s (coefficients typically ~1e-3, sigma2 ~1% on reference fits; weakly identified fits may reach different local optima) and are NOT covered by the parity guarantee. 'css-ml' and 'ml' results are — CSS supplies starting values only.

class pystatistics.timeseries.ARIMASolution(_result)[source]¶

Bases: SolutionReprMixin

Result from fitting an ARIMA model.

Wraps a Result [ARIMAParams] envelope; every datum is exposed via a read-only @property so the public attribute surface is unchanged from the previous flat dataclass.

Parameters:: _result (Result[ARIMAParams])

order¶

(p, d, q) — AR order, differencing order, MA order.

Type:: tuple[int, int, int]

seasonal_order¶

(P, D, Q, m) — seasonal orders and period, or None.

Type:: tuple[int, int, int, int] or None

ar¶

AR coefficients (length p). For seasonal models, these are the non-seasonal AR coefficients only.

Type:: NDArray

ma¶

MA coefficients (length q). For seasonal models, these are the non-seasonal MA coefficients only.

Type:: NDArray

seasonal_ar¶

Seasonal AR coefficients (length P). Empty if non-seasonal.

Type:: NDArray

seasonal_ma¶

Seasonal MA coefficients (length Q). Empty if non-seasonal.

Type:: NDArray

mean¶

Estimated mean of the differenced series (None if include_mean=False).

Type:: float or None

sigma2¶

Estimated innovation variance.

Type:: float

vcov¶

Variance-covariance matrix of the estimated coefficients.

Type:: NDArray

residuals¶

Innovation residuals (length of the differenced series).

Type:: NDArray

fitted_values¶

One-step-ahead fitted values (length of the differenced series).

Type:: NDArray

log_likelihood¶

Maximized log-likelihood value.

Type:: float

aic¶

Akaike information criterion.

Type:: float

aicc¶

Corrected AIC.

Type:: float

bic¶

Bayesian information criterion.

Type:: float

n_obs¶

Length of the original (undifferenced) series.

Type:: int

n_used¶

Number of observations used in estimation (after differencing).

Type:: int

method¶

Estimation method: 'css', 'ml', or 'css-ml'.

Type:: str

converged¶

Whether the optimizer converged.

Type:: bool

n_iter¶

Number of optimizer iterations.

Type:: int

property order: tuple[int, int, int]¶

property seasonal_order: tuple[int, int, int, int] | None¶

property ar: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property ma: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property seasonal_ar: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property seasonal_ma: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property mean: float | None¶

property sigma2: float¶

property vcov: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property residuals: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property fitted_values: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property log_likelihood: float¶

property aic: float¶

property aicc: float¶

property bic: float¶

property n_obs: int¶

property n_used: int¶

property method: str¶

property converged: bool¶

property n_iter: int¶

property xreg_coef: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property xreg_names: tuple[str, ...]¶

property include_drift: bool¶

property xreg: ndarray[tuple[Any, ...], dtype[_ScalarT]] | None¶

property info: dict¶

property timing: dict[str, float] | None¶

property backend_name: str¶

property warnings: tuple[str, ...]¶

property n_params: int¶

Total number of parameters in the reported coefficient vector (AR + MA + seasonal + mean/regression + sigma2).

For a regression-with-ARIMA-errors model the mean is carried as the 'intercept' regression coefficient, so it is counted in xreg_coef rather than via mean.

summary()[source]¶

R-style summary matching stats::arima() output.

Returns:: Multi-line summary string.
Return type:: str

pystatistics.timeseries.arima_batch(Y, *, order=(0, 0, 0), include_mean=True, method='whittle', tol=1e-05, max_iter=300, lr=0.05, backend=None)[source]¶

Fit K independent ARMA(p, d, q) models on the rows of Y.

Parameters:

Y (ArrayLike or torch.Tensor) – Shape (K, n). Each row is an independent time series.
order (tuple[int, int, int]) – (p, d, q). Differencing d is applied per-series before the ARMA fit.
include_mean (bool) – Whether to report the per-series sample mean of the differenced series. Default True. Whittle is centred internally regardless.
method (str) – Only 'whittle' is supported for batch fitting in 1.9.0.
tol (float) – Per-series gradient-norm (L∞) convergence tolerance.
max_iter (int) – Maximum Adam iterations (batched) / maximum per-series L-BFGS-B iterations (CPU loop).
lr (float) – Adam learning rate (GPU path only). 0.05 is the default tuned for typical Whittle NLL curvature — smaller if you see per-series non-convergence, larger if you want faster wall time on easy problems.
backend (str or None) – Compute backend = (device, precision). 'cpu' (default — loop over arima() with method='whittle', float64), 'gpu' (float32, require GPU), 'gpu_fp64' (float64, CUDA only — raises on MPS), 'auto' (GPU-float32 if CUDA present, else CPU loop). A torch.Tensor input routes to its own device automatically (see CONVENTIONS.md).

Returns:

Failed series — a non-stationary Whittle optimum or a per-series optimizer failure, on any backend — have their ar/ma/sigma2/mean rows set to NaN and converged=False, with a UserWarning naming the count (also recorded in .warnings). Identify them with np.isnan(result.sigma2). The contract is identical across backends; see _arima_batch_contract.

Return type:

ARMABatchSolution

Raises:

ConvergenceError – If every series in the batch fails.
ValidationError – On invalid method / order / shape / backend.

class pystatistics.timeseries.ARMABatchSolution(_result)[source]¶

Bases: SolutionReprMixin

Result from a batched ARMA fit.

Wraps a Result [ARMABatchParams] envelope; every datum is exposed via a read-only @property so the public attribute surface is unchanged from the previous flat dataclass.

Parameters:: _result (Result[ARMABatchParams])

order¶

(p, d, q). The d that was applied before the ARMA fit.

Type:: tuple[int, int, int]

ar¶

AR coefficients, shape (K, p).

Type:: NDArray

ma¶

MA coefficients, shape (K, q).

Type:: NDArray

sigma2¶

Innovation variance per series, shape (K,).

Type:: NDArray

mean¶

Per-series sample mean of the differenced series (None if include_mean=False).

Type:: NDArray | None

n_iter¶

Maximum iteration count across all series.

Type:: int

converged¶

Per-series boolean convergence flag, shape (K,). Always False for a failed series (whose estimates are NaN — see arima_batch()); on the float32 GPU path it can also be False on a valid fit whose Adam gradient stayed above tol.

Type:: NDArray

n_series¶

Number of series K.

Type:: int

n_used¶

Length of each (post-differencing) series.

Type:: int

method¶

Which fitter ran: 'whittle-batch-gpu', 'whittle-loop-cpu'.

Type:: str

property order: tuple[int, int, int]¶

property ar: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property ma: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property sigma2: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property mean: ndarray[tuple[Any, ...], dtype[_ScalarT]] | None¶

property n_iter: int¶

property converged: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property n_series: int¶

property n_used: int¶

property method: str¶

property info: dict¶

property timing: dict[str, float] | None¶

property backend_name: str¶

property warnings: tuple[str, ...]¶

summary()[source]¶

Compact summary of a batched ARMA fit across the K series.

Return type:: str

pystatistics.timeseries.forecast_arima(fitted, y_original, *, n_ahead=10, conf_level=(0.8, 0.95), new_xreg=None)[source]¶

Generate forecasts from a fitted ARIMA model.

Matches R’s predict.Arima() / forecast::forecast.Arima().

Parameters:

fitted (ARIMASolution) – A fitted ARIMA model (from arima()).
y_original (ArrayLike) – The original (un-differenced) time series that was passed to arima(). Needed to reverse the differencing.
n_ahead (int) – Forecast horizon (number of steps ahead). Default 10.
conf_level (float or sequence of float) – Prediction-interval confidence level(s) as fractions in (0, 1) (default (0.80, 0.95)). A single float requests one interval; a sequence requests several. Whole-percent values (e.g. 95) are rejected.
new_xreg (ArrayLike, optional) – Future values of the external regressors, shape (n_ahead, k) (or (n_ahead,) for a single regressor), required when the model was fit with xreg. Not needed for drift-only / intercept-only models (those future columns are synthesized). Ignored for models with no regressors.

Returns:

Point forecasts and prediction intervals on the original scale.

Return type:

ARIMAForecast

Raises:

ValidationError – If n_ahead < 1, conf_level is invalid, or new_xreg is missing/misshaped.

class pystatistics.timeseries.ARIMAForecast(mean, se, lower, upper, n_ahead, order)[source]¶

Bases: object

Forecast from a fitted ARIMA model.

Parameters:

mean (ndarray[tuple[Any, ...], dtype[_ScalarT]])
se (ndarray[tuple[Any, ...], dtype[_ScalarT]])
lower (dict[float, ndarray[tuple[Any, ...], dtype[_ScalarT]]])
upper (dict[float, ndarray[tuple[Any, ...], dtype[_ScalarT]]])
n_ahead (int)
order (tuple[int, int, int])

mean¶

Point forecasts on the original (un-differenced) scale, length n_ahead.

Type:: NDArray

se¶

Standard errors of forecasts, length n_ahead.

Type:: NDArray

lower¶

Lower prediction-interval bounds keyed by confidence level as a fraction (e.g. 0.80, 0.95).

Type:: dict[float, NDArray]

upper¶

Upper prediction-interval bounds keyed by confidence level.

Type:: dict[float, NDArray]

n_ahead¶

Forecast horizon.

Type:: int

order¶

The (p, d, q) order of the model.

Type:: tuple[int, int, int]

mean: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

se: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

lower: dict[float, ndarray[tuple[Any, ...], dtype[_ScalarT]]]¶

upper: dict[float, ndarray[tuple[Any, ...], dtype[_ScalarT]]]¶

n_ahead: int¶

order: tuple[int, int, int]¶

summary()[source]¶

Return a human-readable forecast summary.

Returns:: Multi-line table of point forecasts and intervals.
Return type:: str

pystatistics.timeseries.auto_arima(y, *, max_p=5, max_q=5, max_d=2, max_P=2, max_Q=2, max_D=1, period=1, ic='aicc', stepwise=True, allow_drift=True, tol=1e-08, max_iter=1000, method='css-ml', backend=None)[source]¶

Automatic ARIMA model selection.

Simplified version of R’s forecast::auto.arima().

For stepwise=True the Hyndman–Khandakar (2008) algorithm is used: start from a set of initial candidates and greedily explore neighbouring orders. For seasonal models (period > 1) the seasonal AR/MA orders (P, Q) are searched alongside (p, q).

For stepwise=False an exhaustive grid search over all p = 0 .. max_p, q = 0 .. max_q — and, for seasonal models, P = 0 .. max_P, Q = 0 .. max_Q — combinations is performed (much slower but thorough).

Parameters:

y (ArrayLike) – Time series.
max_p (int) – Maximum non-seasonal AR / MA orders.
max_q (int) – Maximum non-seasonal AR / MA orders.
max_d (int) – Maximum non-seasonal differencing order.
max_P (int) – Maximum seasonal orders.
max_Q (int) – Maximum seasonal orders.
max_D (int) – Maximum seasonal orders.
period (int) – Seasonal period (1 = non-seasonal).
ic (str) – Information criterion: 'aic', 'aicc', or 'bic'.
stepwise (bool) – Use stepwise search (default) or grid search.
allow_drift (bool) – Allow a drift (linear-trend) term to be selected when the total differencing order d + D == 1 — the models R reports “with drift” (the drift-allowance option of forecast::auto.arima). Each visited order is fit with and without drift and the better information criterion wins; the chosen model exposes include_drift and a 'drift' entry in best_model.xreg_coef. Default True.
tol (float) – Convergence tolerance passed to arima().
max_iter (int) – Maximum iterations passed to arima().
method (str) – Estimation method forwarded to every candidate fit. Default 'css-ml' matches R. Use 'whittle' with backend='gpu' to route each candidate through the frequency-domain GPU path.
backend (str or None) – Backend forwarded to every candidate fit. Default None → CPU (R-reference path). Pass 'gpu' or 'auto' to opt into the GPU path; only meaningful when the candidate fits actually support it (e.g. method='whittle').

Returns:

Best model and search history.

Return type:

AutoARIMASolution

Raises:

ValidationError – If inputs are invalid.
ConvergenceError – If no model converges.

class pystatistics.timeseries.AutoARIMASolution(_result)[source]¶

Bases: SolutionReprMixin

Result from automatic ARIMA order selection.

Wraps a Result [AutoARIMAParams] envelope; every datum is exposed via a read-only @property so the public attribute surface is unchanged from the previous flat dataclass.

Parameters:: _result (Result[AutoARIMAParams])

best_model¶

The fitted model with the best information criterion.

Type:: ARIMASolution

best_order¶

The (p, d, q) order of the best model.

Type:: tuple[int, int, int]

best_seasonal¶

The (P, D, Q, m) seasonal order, or None.

Type:: tuple[int, int, int, int] | None

best_aic¶

Value of the chosen information criterion for the best model.

Type:: float

models_fitted¶

Total number of models successfully evaluated.

Type:: int

search_results¶

(order, ic_value) pairs for every model tried (including those that failed, recorded with inf). For seasonal searches order is the pair ((p, d, q), (P, D, Q, m)).

Type:: list[tuple[tuple, float]]

summary()[source]¶

Human-readable summary of the search.

Return type:: str

property best_model: object¶

property best_order: tuple[int, int, int]¶

property best_seasonal: tuple[int, int, int, int] | None¶

property best_aic: float¶

property models_fitted: int¶

property search_results: list[tuple[tuple, float]]¶

property info: dict¶

property timing: dict[str, float] | None¶

property backend_name: str¶

property warnings: tuple[str, ...]¶

summary()[source]¶

Return a human-readable summary of the search.

Returns:: Multi-line summary.
Return type:: str

pystatistics.timeseries.decompose(x, period, *, kind='additive')[source]¶

Classical time series decomposition.

Matches R’s stats::decompose().

Algorithm¶

Trend via centered moving average of length period.
De-trend: subtract (additive) or divide (multiplicative).
Seasonal: average de-trended values at each seasonal position, then center so they sum to zero (additive) or average to 1 (multiplicative).
Residual: x - trend - seasonal (additive) or x / (trend * seasonal) (multiplicative).

param x:: Time series (1-D). Length must be >= 2 * period.
type x:: ArrayLike
param period:: Seasonal period (>= 2).
type period:: int
param kind:: 'additive' or 'multiplicative'.
type kind:: str
rtype:: DecompositionSolution
raises ValidationError:: On invalid inputs.

Parameters:

x (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str])
period (int)
kind (str)

Return type:

DecompositionSolution

pystatistics.timeseries.stl(x, period, *, seasonal_window='periodic', seasonal_degree=0, trend_window=None, trend_degree=1, lowpass_window=None, lowpass_degree=None, seasonal_jump=None, trend_jump=None, lowpass_jump=None, robust=False, n_inner=None, n_outer=None)[source]¶

Seasonal-Trend decomposition using Loess (STL).

Matches R’s stats::stl (Cleveland et al., 1990): identical parameters produce identical seasonal/trend/remainder components up to floating-point noise. Two interface divergences, both deliberate: seasonal_window defaults to "periodic" (R requires it explicitly), and invalid spans raise instead of being silently rounded up as R does.

Parameters:

x (ArrayLike) – Time series (1-D, finite). Length must exceed 2 * period.
period (int) – Seasonal period (>= 2), e.g. 12 for monthly data.
seasonal_window (int, "periodic", or None) – Loess span for cycle-subseries smoothing (odd, >= 3), or "periodic" for an exactly periodic seasonal (equivalent to a span of 10*n + 1 with degree 0, followed by cycle-position averaging). Default "periodic"; None also selects the default, as for every other window parameter. R: s.window.
seasonal_degree (int) – Loess degree (0 or 1) for the seasonal smoother. Default 0, matching R. Must be 0 when seasonal_window="periodic".
trend_window (int or None) – Loess span for the trend (odd, >= 3). Default nextodd(ceil(1.5*period / (1 - 1.5/seasonal_window))), as in R.
trend_degree (int) – Loess degree (0 or 1) for the trend. Default 1, matching R.
lowpass_window (int or None) – Loess span of the low-pass filter. Default nextodd(period).
lowpass_degree (int or None) – Loess degree of the low-pass filter. Default: trend_degree.
seasonal_jump (int or None) – Evaluation strides: each loess is evaluated every jump-th point and linearly interpolated between, exactly as in R. Defaults ceil(window/10) (R’s defaults). Pass 1 to evaluate the loess at every point with no interpolation (R’s own interpolation is a speed shortcut, Cleveland et al. 1990, sec. 3.4; with the default strides the output matches R’s bit-for-bit up to float noise).
trend_jump (int or None) – Evaluation strides: each loess is evaluated every jump-th point and linearly interpolated between, exactly as in R. Defaults ceil(window/10) (R’s defaults). Pass 1 to evaluate the loess at every point with no interpolation (R’s own interpolation is a speed shortcut, Cleveland et al. 1990, sec. 3.4; with the default strides the output matches R’s bit-for-bit up to float noise).
lowpass_jump (int or None) – Evaluation strides: each loess is evaluated every jump-th point and linearly interpolated between, exactly as in R. Defaults ceil(window/10) (R’s defaults). Pass 1 to evaluate the loess at every point with no interpolation (R’s own interpolation is a speed shortcut, Cleveland et al. 1990, sec. 3.4; with the default strides the output matches R’s bit-for-bit up to float noise).
robust (bool) – Enable outer-loop robustness iterations (bisquare weights on the remainder). Changes the n_inner/n_outer defaults to R’s (1, 15) from (2, 0).
n_inner (int or None) – Inner-loop passes (>= 1). Default: 1 if robust else 2. R: inner.
n_outer (int or None) – Robustness iterations (>= 0). Default: 15 if robust else 0. R: outer.

Returns:

With seasonal + trend + residual == observed exactly. info records the resolved windows/degrees/jumps, the iteration counts, and the final robustness weights.

Return type:

DecompositionSolution

Raises:

ValidationError – On invalid input, spans that are even or < 3, degrees outside {0, 1}, or a series with fewer than two full periods.

class pystatistics.timeseries.DecompositionSolution(_result)[source]¶

Bases: SolutionReprMixin

Result from time series decomposition.

Wraps a Result [DecompositionParams] envelope; every datum is exposed via a read-only @property so the public attribute surface is unchanged from the previous flat dataclass.

Parameters:: _result (Result[DecompositionParams])

observed¶

Original time series.

Type:: NDArray

trend¶

Trend component. May contain NaN at edges for classical decomposition.

Type:: NDArray

seasonal¶

Seasonal component (repeats with the given period).

Type:: NDArray

residual¶

Remainder after removing trend and seasonal.

Type:: NDArray

period¶

Seasonal period used.

Type:: int

kind¶

'additive' or 'multiplicative'.

Type:: str

method¶

'classical' or 'stl'.

Type:: str

property observed: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property trend: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property seasonal: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property residual: ndarray[tuple[Any, ...], dtype[_ScalarT]]¶

property period: int¶

property kind: str¶

property method: str¶

property info: dict¶

property timing: dict[str, float] | None¶

property backend_name: str¶

property warnings: tuple[str, ...]¶

summary()[source]¶

Return a human-readable summary of the decomposition.

Returns:: Multi-line summary string.
Return type:: str