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:

  1. Regress x on deterministic terms (constant, or constant+trend).

  2. Compute partial sums S_t = sum_{i=1}^{t} e_i of residuals.

  3. 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.

Variables:

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

  • lags (NDArray) – Lag indices corresponding to each acf value.

  • n_obs (int) – Number of observations in the original series.

  • conf_level (float) – Confidence level used for the confidence bands.

  • ci_upper (NDArray) – Upper confidence bound per lag.

  • ci_lower (NDArray) – Lower confidence bound per lag.

  • kind (str) – ‘correlation’ for ACF or ‘partial’ for PACF.

Parameters:

_result (Result[ACFParams])

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.

Variables:

  • statistic (float) – Test statistic value.

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

  • method (str) – Name of the test, e.g. ‘Augmented Dickey-Fuller’, ‘KPSS’.

  • alternative (str) – Alternative hypothesis description, e.g. ‘stationary’ for ADF, ‘unit root’ for KPSS.

  • n_lags (int) – Number of lags used in the test.

  • n_obs (int) – Number of observations used (after differencing/lag adjustments).

  • critical_values (dict[str, float] | None) – Critical values at standard significance levels, e.g. {‘1%’: -3.43, ‘5%’: -2.86, ‘10%’: -2.57}.

Parameters:

_result (Result[StationarityParams])

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.

Variables:

  • mean (NDArray) – Point forecasts of length n_ahead.

  • lower (dict[float, NDArray]) – Lower prediction-interval bounds keyed by confidence level as a fraction (e.g. {0.8: ..., 0.95: ...}).

  • upper (dict[float, NDArray]) – Upper prediction-interval bounds keyed by confidence level.

  • n_ahead (int) – Forecast horizon.

  • model (ETSSpec) – Model specification used.

  • fitted (ETSSolution) – The underlying fitted 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
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.

Variables:

  • error (str) – Error type: 'A' (additive) or 'M' (multiplicative).

  • trend (str) – Trend type: 'N' (none), 'A' (additive), or 'Ad' (damped).

  • season (str) – Season type: 'N' (none), 'A' (additive), or 'M' (multiplicative).

  • period (int) – Seasonal period (1 when no seasonal component).

  • damped (bool) – True when trend is 'Ad'.

Parameters:
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):

  1. Difference the series d times.

  2. Fit ARMA(p, q) to the differenced series.

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

  1. Seasonally difference D times (lag m).

  2. Difference d times.

  3. 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:

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.

Variables:

  • order (tuple[int, int, int]) – (p, d, q) — AR order, differencing order, MA order.

  • seasonal_order (tuple[int, int, int, int] or None) – (P, D, Q, m) — seasonal orders and period, or None.

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

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

  • seasonal_ar (NDArray) – Seasonal AR coefficients (length P). Empty if non-seasonal.

  • seasonal_ma (NDArray) – Seasonal MA coefficients (length Q). Empty if non-seasonal.

  • mean (float or None) – Estimated mean of the differenced series (None if include_mean=False).

  • sigma2 (float) – Estimated innovation variance.

  • vcov (NDArray) – Variance-covariance matrix of the estimated coefficients.

  • residuals (NDArray) – Innovation residuals (length of the differenced series).

  • fitted_values (NDArray) – One-step-ahead fitted values (length of the differenced series).

  • log_likelihood (float) – Maximized log-likelihood value.

  • aic (float) – Akaike information criterion.

  • aicc (float) – Corrected AIC.

  • bic (float) – Bayesian information criterion.

  • n_obs (int) – Length of the original (undifferenced) series.

  • n_used (int) – Number of observations used in estimation (after differencing).

  • method (str) – Estimation method: 'css', 'ml', or 'css-ml'.

  • converged (bool) – Whether the optimizer converged.

  • n_iter (int) – Number of optimizer iterations.

Parameters:

_result (Result[ARIMAParams])

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).

    On the GPU path the per-iteration likelihood is torch.compile``d by default (measured ~1.75-2.14x end-to-end at K=1000), which costs a ONE-TIME ~1 s compilation on the first call in a process — a strictly one-shot call is slower end-to-end; repeated calls win from the second onward. Falls back to the identical eager computation (with a ``RuntimeWarning) if the compile machinery is unavailable or fails; the resolved mode is reported in Solution.info['nll_compiled'].

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:
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.

Variables:

  • order (tuple[int, int, int]) – (p, d, q). The d that was applied before the ARMA fit.

  • ar (NDArray) – AR coefficients, shape (K, p).

  • ma (NDArray) – MA coefficients, shape (K, q).

  • sigma2 (NDArray) – Innovation variance per series, shape (K,).

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

  • n_iter (int) – Maximum iteration count across all series.

  • converged (NDArray) – 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.

  • n_series (int) – Number of series K.

  • n_used (int) – Length of each (post-differencing) series.

  • method (str) – Which fitter ran: 'whittle-batch-gpu', 'whittle-loop-cpu'.

Parameters:

_result (Result[ARMABatchParams])

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.

Variables:

  • mean (NDArray) – Point forecasts on the original (un-differenced) scale, length n_ahead.

  • se (NDArray) – Standard errors of forecasts, length n_ahead.

  • lower (dict[float, NDArray]) – Lower prediction-interval bounds keyed by confidence level as a fraction (e.g. 0.80, 0.95).

  • upper (dict[float, NDArray]) – Upper prediction-interval bounds keyed by confidence level.

  • n_ahead (int) – Forecast horizon.

  • order (tuple[int, int, int]) – The (p, d, q) order of the 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]
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:
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.

Variables:

  • best_model (ARIMASolution) – The fitted model with the best information criterion.

  • best_order (tuple[int, int, int]) – The (p, d, q) order of the best model.

  • best_seasonal (tuple[int, int, int, int] | None) – The (P, D, Q, m) seasonal order, or None.

  • best_aic (float) – Value of the chosen information criterion for the best model.

  • models_fitted (int) – Total number of models successfully evaluated.

  • search_results (list[tuple[tuple, float]]) – (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)).

Parameters:

_result (Result[AutoARIMAParams])

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

  1. Trend via centered moving average of length period.

  2. De-trend: subtract (additive) or divide (multiplicative).

  3. Seasonal: average de-trended values at each seasonal position, then center so they sum to zero (additive) or average to 1 (multiplicative).

  4. 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:
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.

Variables:

  • observed (NDArray) – Original time series.

  • trend (NDArray) – Trend component. May contain NaN at edges for classical decomposition.

  • seasonal (NDArray) – Seasonal component (repeats with the given period).

  • residual (NDArray) – Remainder after removing trend and seasonal.

  • period (int) – Seasonal period used.

  • kind (str) – 'additive' or 'multiplicative'.

  • method (str) – 'classical' or 'stl'.

Parameters:

_result (Result[DecompositionParams])

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