Univariate linear models#

This page covers the linear Gaussian models for the conditional mean of a single series: autoregressions, the Box-Jenkins ARMA family with differencing, seasonal blocks and regressors, fractionally integrated ARMA, and Harvey’s structural (unobserved-components) models. They answer “what is the persistence, memory and deterministic structure of this series, and which specification does the likelihood prefer?”. They are not for a time-varying variance (see Volatility models), regime or threshold dynamics (see Nonlinear and regime-switching univariate models), or several series at once (see Vector autoregressions).

The family at a glance#

Class

Module

Model in one line

Use when

AR

cultivars.univariate.autoregression

\(y_t = c + \delta t + \sum_i \phi_i y_{t-i} + \varepsilon_t\), by OLS (CSS) or exact ML

a pure autoregression, a fast OLS fit, or a linear trend

ARMA, ARIMA, SARIMA, ARIMAX, SARIMAX

cultivars.univariate.box_jenkins

multiplicative seasonal ARIMA, optionally a regression with ARIMA errors, by exact ML through the Kalman filter

moving-average terms, integer differencing, seasonality, exogenous regressors

ARFIMA

cultivars.univariate.fractional_integration

\(\phi(L)(1-L)^d(y_t-\mu) = \theta(L)\varepsilon_t\), \(|d| < 0.5\)

hyperbolically decaying autocorrelations (long memory)

UnobservedComponents

cultivars.univariate.unobserved_components

level/trend + damped cycle + trigonometric seasonal + irregular, by exact ML and the Kalman smoother

a decomposition with bands, model-based detrending, and the only built-in forecast in this family

Every model follows the package’s three-object discipline (see Architecture): the constructor validates and freezes a specification, fit() estimates it, and the returned result is an immutable record. Printing a result prints its summary table.

Autoregressions: AR#

\[y_t = c + \delta t + \sum_{i=1}^{p} \phi_i\, y_{t-i} + \varepsilon_t, \qquad \varepsilon_t \sim N(0, \sigma^2).\]

Two estimators sit behind method. "css" (the default) conditions on the first \(p\) observations and solves the lagged regression by OLS, so the likelihood is the conditional one on \(n - p\) observations. "exact" maximizes the exact Gaussian likelihood through the companion state-space form, keeping the search inside the stationary region with the partial-autocorrelation map of Monahan (1984); its effective sample is all \(n\) observations (Hamilton, 1994, ch. 5).

import numpy as np
from cultivars.univariate.autoregression import AR

rng = np.random.default_rng(0)
e = rng.standard_normal(600)
z = np.zeros(600)
for t in range(2, 600):
    z[t] = 0.5 + 0.6 * z[t - 1] - 0.2 * z[t - 2] + e[t]
y = z[100:]                      # drop the start-up transient; n = 500

css = AR(y, order=2).fit()
exact = AR(y, order=2, method="exact").fit()
print(css.nobs, exact.nobs)
print({k: round(v, 3) for k, v in css.params.items()})
print(css.is_stationary, round(css.stability.max_modulus, 3))

498 500
{'const': 0.49, 'ar.L1': 0.556, 'ar.L2': -0.205, 'sigma2': 1.003}
True 0.453

Specification choices:

order

Non-negative integer \(p\). order=0 fits the deterministic terms alone.

trend

"n" (none), "c" (constant, default) or "ct" (constant and linear trend). "ct" is available under CSS only; combining it with method="exact" raises SpecificationError at construction.

method

"css" (OLS, \(n - p\) observations) or "exact" (exact ML, \(n\) observations).

What the result (ARResult) gives you:

  • params: a dict in table order. Names are const, trend, ar.L1 .. ar.Lp and sigma2; absent terms are omitted.

  • The same numbers as fields: const (None when trend="n"), trend_coeff, ar_params (an array) and sigma2.

  • llf, nobs, n_params (counts sigma2), information_criteria, compare, likelihood_ratio_test.

  • stability (companion eigenvalues, max_modulus, n_unit_roots, n_explosive) and is_stationary.

  • fittedvalues and resid on the effective sample, endog in full, to_pandas / to_polars / to_dict.

  • simulate(n, seed=..., burn=...).

const is the intercept inside the recursion. The implied mean is \(c / (1 - \sum_i \phi_i)\):

print(round(css.const / (1 - css.ar_params.sum()), 3), round(y.mean(), 3))

0.756 0.759

Box-Jenkins models: ARMA to SARIMAX#

\[\phi(L)\,\Phi(L^s)\,\bigl(w_t - x_t'\beta\bigr) = \theta(L)\,\Theta(L^s)\,\varepsilon_t, \qquad w_t = (1 - L)^d (1 - L^s)^D y_t,\]

with \(\phi(L) = 1 - \phi_1 L - \cdots\), the moving average in the plus convention \(\theta(L) = 1 + \theta_1 L + \cdots\), and \(x_t\) the deterministic terms followed by any exogenous regressors (Box, Jenkins, Reinsel and Ljung, 2015). Differencing is applied before estimation; the ARMA part is cast in Harvey’s state-space form and the exact likelihood is evaluated by the Kalman filter (Harvey, 1989; Durbin and Koopman, 2012). Stationarity and invertibility are imposed through the partial-autocorrelation reparameterization.

SARIMAX is the one real specification; the other four classes are narrowing constructors that fix or forbid arguments. All five return ARMAResult.

Class

Signature

Restriction

ARMA

ARMA(endog, *, order=(p, q), trend="c")

d = 0, no seasonal block, no regressors

ARIMA

ARIMA(endog, *, order=(p, d, q), trend="c")

no seasonal block, no regressors

SARIMA

SARIMA(endog, *, order, seasonal_order=(P, D, Q, s), trend="c")

seasonal block required, no regressors

ARIMAX

ARIMAX(endog, *, order, exog, trend="c")

regressors required (None or zero columns raise)

SARIMAX

SARIMAX(endog, *, order, seasonal_order=(0, 0, 0, 0), trend="c", exog=None)

none

An ARMA(1, 1) on simulated data, with the moving-average sign convention visible:

from cultivars.univariate.box_jenkins import ARIMA, ARIMAX, ARMA, SARIMA

rng = np.random.default_rng(2)
e = rng.standard_normal(500)
x = np.zeros(500)
for t in range(1, 500):
    x[t] = 0.8 * x[t - 1] + e[t] - 0.4 * e[t - 1]

arma = ARMA(x, order=(1, 1)).fit()
print(arma.order, arma.nobs, arma.n_params)
print({k: round(v, 3) for k, v in arma.params.items()})
print(arma.is_stationary, arma.is_invertible)

(1, 0, 1) 500 4
{'const': -0.158, 'ar.L1': 0.739, 'ma.L1': -0.345, 'sigma2': 1.022}
True True

A simulated \(e_t - 0.4 e_{t-1}\) comes back as ma.L1 near \(-0.4\).

Specification choices:

order

(p, d, q) ((p, q) for ARMA), non-negative integers.

seasonal_order

(P, D, Q, s); s must be at least 2 whenever P, D or Q is positive.

trend

"n", "c" or "ct", applied on the differenced scale. With d = 1 and trend="c" the constant is the drift of the levels.

exog

(n,) or (n, k) regressors, differenced alongside endog. This is a regression with ARIMA errors, not a transfer function.

Parameter names extend the AR scheme: const, trend, x1 .. xk (in the order of the beta array), ar.L1 .., ma.L1 .., seasonal terms at their lag multiples (ar.S.L12, ma.S.L12, ma.S.L24 ..), then sigma2. The fields are ar_params, ma_params, seasonal_ar_params, seasonal_ma_params, beta, sigma2, plus order, seasonal_order, seasonal_period, trend and k_exog.

A seasonal model and an integrated one:

s = np.zeros(500)
for t in range(13, 500):
    s[t] = 0.5 * s[t - 1] + 0.4 * s[t - 12] - 0.2 * s[t - 13] + e[t]
seasonal = SARIMA(s, order=(1, 0, 0), seasonal_order=(1, 0, 0, 12)).fit()
print(list(seasonal.params))
print(seasonal.stability.eigenvalues.shape)

levels = np.cumsum(0.3 + e)            # random walk with drift 0.3
rw = ARIMA(levels, order=(0, 1, 0)).fit()
print(rw.nobs, round(rw.beta[0], 3))

['const', 'ar.L1', 'ar.S.L12', 'sigma2']
(13,)
499 0.246

Stationarity and invertibility are judged on the multiplied polynomials \(\phi(L)\Phi(L^s)\) and \(\theta(L)\Theta(L^s)\), so the seasonal AR(1) at \(s = 12\) over an AR(1) has 13 companion eigenvalues. stability and invertibility are StabilityAssessment records.

A regression with ARIMA errors recovers a levels relationship even when the error is integrated:

X = rng.standard_normal((500, 2))
ylev = X @ [2.0, -1.0] + np.cumsum(rng.standard_normal(500))
reg = ARIMAX(ylev, order=(0, 1, 0), exog=X, trend="n").fit()
print(list(reg.params), reg.beta.round(2))

['x1', 'x2', 'sigma2'] [ 2.01 -0.91]

Long memory: ARFIMA#

\[\phi(L)\,(1 - L)^d\,(y_t - \mu) = \theta(L)\,\varepsilon_t, \qquad (1 - L)^d = \sum_{k \ge 0} \frac{\Gamma(k - d)}{\Gamma(-d)\,\Gamma(k + 1)} L^k.\]

For \(0 < d < 0.5\) the process is stationary with autocorrelations decaying as \(k^{2d-1}\) (Granger and Joyeux, 1980; Hosking, 1981). ARFIMA estimates \(\mu\), \(d\), the short-memory block and \(\sigma^2\) jointly: the series is fractionally differenced with the filter cut at truncation terms and the ARMA block’s exact Kalman likelihood is evaluated on the result. The search starts at a local Whittle estimate of \(d\) and confines \(d\) to \((-0.499, 0.499)\) through a scaled tanh.

from scipy.special import gammaln
from cultivars.univariate.fractional_integration import ARFIMA

rng = np.random.default_rng(3)
n, d = 500, 0.3
k = np.arange(n + 300)
psi = np.exp(gammaln(k + d) - gammaln(d) - gammaln(k + 1))  # (1-L)^(-d)
lm = np.convolve(rng.standard_normal(n + 300), psi)[300:300 + n]

arfima = ARFIMA(lm, order=(0, 0)).fit()
print({k: round(v, 3) for k, v in arfima.params.items()})
print(arfima.truncation, arfima.is_stationary, arfima.has_long_memory)

short = ARMA(lm, order=(1, 0)).fit()
table = arfima.compare(short, criterion="bic")
print([row[0] for row in table.rows])

{'mean': 0.181, 'd': 0.319, 'sigma2': 0.983}
500 True True
['ARFIMA(0, d, 0)', 'ARIMA(1, 0, 0)']

Specification choices: order=(p, q) for the short-memory block; trend is "c" (estimate \(\mu\), reported as mean) or "n"; truncation defaults to the sample size. The result (ARFIMAResult) carries d, mean, ar_params, ma_params, sigma2, truncation, the properties is_stationary (requires \(|d| < 0.5\) and a stable AR block), is_mean_reverting (\(d < 1\)), has_long_memory (\(|d| > 10^{-3}\), a threshold not a test), plus stability, invertibility and the comparison methods. fittedvalues and resid are on the fractionally differenced scale. There is no simulate.

Structural decomposition: UnobservedComponents#

\[y_t = \mu_t + \psi_t + \gamma_t + \varepsilon_t, \qquad \varepsilon_t \sim N(0, \sigma^2_\varepsilon),\]

with the trend \(\mu_t\) a local level, a local linear trend or a smooth trend, \(\psi_t\) a damped stochastic cycle with estimated period, and \(\gamma_t\) a trigonometric seasonal whose harmonics share one variance (Harvey, 1989; Harvey and Jaeger, 1993). The variances (and the cycle’s damping and frequency) are estimated by exact ML through the Kalman filter under an approximate-diffuse prior (variance \(10^6\)) on the nonstationary states; every component is reported as the full-sample smoothed estimate with its smoothed standard deviation.

Argument

Meaning

Default

trend="level"

local level \(\mu_{t+1} = \mu_t + \eta_t\)

default

trend="lltrend"

local linear trend: level and slope both random walks

trend="smooth"

integrated random walk (slope moves, level has no innovation); the smoothed level is the HP trend at \(\lambda = \sigma^2_\varepsilon / \sigma^2_\zeta\)

cycle=True

damped cycle, \(\rho\) and \(\lambda_c\) estimated; the period is confined to roughly 2 to 126 observations

False

seasonal=s

trigonometric seasonal of period s (at least 2)

None

 
from cultivars.univariate.unobserved_components import UnobservedComponents

rng = np.random.default_rng(5)
n = 200
level = np.cumsum(0.1 * rng.standard_normal(n))
pattern = np.tile([1.0, -0.5, -1.0, 0.5], n // 4)
q = level + pattern + 0.3 * rng.standard_normal(n)

uc = UnobservedComponents(q, seasonal=4).fit()
print(uc.n_params, round(uc.sigma2_irregular, 3), round(uc.sigma2_level, 3))
print(list(uc.to_pandas().columns))
print(bool(np.allclose(uc.signal + uc.irregular, q)))

mean, std = uc.forecast(4)
print(mean.round(2), std.round(2))

3.0 0.067 0.018
['observed', 'level', 'seasonal', 'irregular']
True
[-0.49 -1.92 -2.37 -1.04] [0.35 0.37 0.39 0.41]

The result (UnobservedComponentsResult) exposes:

  • The estimated variances as fields: sigma2_irregular, sigma2_level, sigma2_slope, sigma2_cycle, sigma2_seasonal, and cycle_rho, cycle_freq (each None when the component is absent), plus cycle_period.

  • Smoothed paths level, level_std, slope, cycle, seasonal, signal (everything but the irregular) and irregular. Asking for a component the specification lacks raises SpecificationError.

  • forecast(steps) returning a (mean, std) pair of arrays; the standard errors include state uncertainty and the irregular, not parameter uncertainty.

  • state_space: the fitted system as a LinearGaussianSSM, re-applicable to new data (see State-space models and filtering).

  • simulate(n, seed=..., burn=0): unit-root components start at zero, not from the diffuse prior.

  • llf, nobs, n_params, information_criteria, compare, likelihood_ratio_test, summary.

Comparing, testing and exporting#

Every result on this page shares one comparison surface. information_criteria is a named tuple (aic, bic, hqic) computed from llf, nobs and n_params. compare(*others, criterion="aic") returns a ranked SummaryTable (columns model, llf, k, the criterion and its gap to the best); there are deliberately no < / > operators. likelihood_ratio_test(unrestricted) treats the caller as the restricted model and returns a record with statistic, df and pvalue. Both refuse results with different nobs.

Order selection the way the package intends it, on a common sample:

fits = [AR(y, order=p, method="exact").fit() for p in (1, 2, 3)]
table = fits[0].compare(*fits[1:], criterion="bic")
print(table.columns)
for row in table.rows:
    print(row[0], row[3], row[4])

lr = fits[1].likelihood_ratio_test(fits[2])
print(lr.df, round(lr.statistic, 3), round(lr.pvalue, 3))

('model', 'llf', 'k', 'BIC', 'dBIC')
AR(2) 1445.784 0.000
AR(3) 1451.739 5.956
AR(1) 1460.832 15.048
1 0.259 0.611

to_pandas(index=None), to_polars() and to_dict() return the aligned per-observation series. For AR, ARMA-family and ARFIMA results the columns are observed, fitted and resid, trimmed to the effective sample; for UnobservedComponents they are the decomposition. The default index is 0 .. nobs - 1 of the effective sample, so pass index= to align with dates:

frame = css.to_pandas(index=np.arange(2, 500))
print(frame.shape, frame.index[0])

(498, 3) 2

Simulating from a fit#

AR, the Box-Jenkins results and UnobservedComponents provide simulate(n=200, *, seed=None, burn=...), which returns an (n,) array at the fitted parameters with Gaussian innovations. AR adds the mean path implied by the deterministic terms and refuses a fit whose AR coefficients sum to one. The Box-Jenkins simulator generates the stationary ARMA part, adds the deterministic block on the differenced scale and integrates from zero, so the first \(d + sD\) entries are exactly zero; it refuses fits with regressors.

print(css.simulate(3, seed=0).round(3))
path = ARIMA(levels, order=(1, 1, 0)).fit().simulate(4, seed=0)
print(path[0], path.shape)

[0.34  1.983 1.36 ]
0.0 (4,)

How the pieces connect#

Nuances and gotchas#

  • No forecast() on AR, ARMA-family or ARFIMA results. ARResult, ARMAResult and ARFIMAResult have no forecast or predict method. Compute point forecasts from the estimates (for an AR with trend="c", iterate \(c + \sum_i \phi_i y_{T+h-i}\)), or fit UnobservedComponents when you need forecasts with standard errors out of the box.

    print(hasattr(css, "forecast"), hasattr(arma, "forecast"),
      hasattr(arfima, "forecast"), hasattr(uc, "forecast"))

    False False False True

  • No standard errors. None of the four results reports standard errors, t-statistics or confidence intervals; each summary says “Standard errors are not yet available for this estimator.” Use likelihood_ratio_test for nested hypotheses.

  • Effective sample sizes differ. nobs is \(n - p\) for AR(method="css"), \(n\) for AR(method="exact"), \(n - d - sD\) for the Box-Jenkins models (MA and AR lags do not shorten it; the exact likelihood uses every differenced observation), \(n\) for ARFIMA and UnobservedComponents. compare and likelihood_ratio_test raise DimensionError on differing nobs, so AR orders fitted by CSS cannot be ranked against each other: use method="exact".

    from cultivars.exceptions import DimensionError
try:
    css.compare(AR(y, order=3).fit())
except DimensionError as err:
    print(err)

    information criteria are only comparable across a common sample; got differing observation counts [497, 498].

  • Mixed scales in to_pandas after differencing. For an ARIMA with d > 0 (or D > 0) the observed column is the level series trimmed to the last nobs rows, while fitted and resid are on the differenced scale, so observed - fitted is not resid. Difference endog yourself to compare like with like.

    df = rw.to_pandas()
print(round(df["observed"].iloc[-1], 2), round(df["fitted"].iloc[-1], 2))
print(bool(np.allclose(np.diff(levels) - rw.fittedvalues, rw.resid)))

    123.29 0.25
True

  • Two intercept conventions. ARResult.const is the recursion intercept \(c\); ARMAResult.beta[0] (params["const"]) is the mean of the differenced series. On the same AR(1) data the two differ by the factor \(1 - \phi\):

    ar1 = AR(y, order=1).fit()
arma10 = ARMA(y, order=(1, 0)).fit()
print(round(ar1.const / (1 - ar1.ar_params[0]), 2), round(arma10.beta[0], 2))

    0.76 0.76

  • Stability verdicts on constrained fits verify the transform. The Box-Jenkins, exact-AR and ARFIMA estimators cannot leave the stationary and invertible region, so is_stationary is True even on a random walk; read stability.max_modulus (close to one) instead. Only AR(method="css") can report a non-stationary fit.

    walk = ARMA(np.cumsum(rng.standard_normal(400)), order=(1, 0)).fit()
print(walk.is_stationary, round(walk.stability.max_modulus, 3))

    True 0.984

  • Near-cancelling ARMA roots are weakly identified, and the optimizer has one start. When \(\phi_1 + \theta_1\) is small the ARMA(1, 1) likelihood is nearly flat along \(\phi_1 = -\theta_1\), and the fit can stop at a common-factor point. On the series below (true \(\phi_1 = 0.6\), \(\theta_1 = -0.4\)) the estimate is a pair of almost cancelling roots; an independent exact-likelihood maximization of the same data reaches a slightly higher llf (by about 0.03) at \((0.71, -0.62)\). Compare against the AR(1) and white-noise fits before reading an ARMA(1, 1) coefficient.

    rng1 = np.random.default_rng(1)
e1 = rng1.standard_normal(500)
x1 = np.zeros(500)
for t in range(1, 500):
    x1[t] = 0.6 * x1[t - 1] + e1[t] - 0.4 * e1[t - 1]
weak = ARMA(x1, order=(1, 1)).fit()
print(weak.ar_params.round(3), weak.ma_params.round(3))

    [0.986] [-0.965]

  • ARFIMA’s \(d\) is boundary-capped. \(d\) is confined to \((-0.499, 0.499)\), so a unit-root series returns d near 0.499 with is_stationary=True. Treat \(d\) at the boundary as a request for integer differencing (fit ARIMA with d=1, or ARFIMA on the differences). \(d\) and a near-unit AR root trade off in short samples; keep order small.

  • Exogenous regressors are differenced too. In ARIMAX / SARIMAX with d > 0 the regressors are differenced with endog; supply them in levels when the relationship is in levels. A result with regressors refuses simulate.

  • Linear trend under exact ML is refused. AR(y, order=p, trend="ct", method="exact") raises SpecificationError.

  • Unobserved-components likelihoods are approximate-diffuse. llf is comparable across specifications on one sample, not with software that uses an exact-diffuse initialization. The result has no params dict, fittedvalues or resid; read the named variance fields. Fits with cycle=True are noticeably slower (two optimizer starts).

  • Likelihood-ratio tests on variances are non-standard. A test of UnobservedComponents(y) against UnobservedComponents(y, seasonal=4) puts a variance on the boundary under the null; the reported \(\chi^2\) p-value is conservative.

See also#