ARIMA in Python: ARIMA, seasonal ARIMA and ARMAX models

Most economic time series remember their past. This quarter's output, this month's sales or today's price depend on recent values and on recent surprises. An ARIMA model describes that memory with a handful of parameters. The fitted model summarizes the dynamics of the series, tells you whether a shock fades or persists, and produces forecasts with error bands. The same machinery fixes a common problem in time-series regressions: autocorrelated errors. A regression with ARMA disturbances keeps the regression and models the error dynamics instead of ignoring them.

What an ARIMA model estimates

An ARIMA(p, d, q) model has three parts. The series is first differenced d times until it has a stable mean and variance. The differenced series then follows an ARMA process: an autoregressive (AR) part with p lags of the series itself and a moving-average (MA) part with q lags of past shocks. With regressors in the mean, OpenEconometrics writes the model as a regression whose disturbance follows the ARIMA process:

y_t = x_t'b + u_t
φ(L) (1 − L)^d u_t = θ(L) ε_t,        ε_t ~ N(0, σ²)
φ(L) = 1 − φ₁L − … − φ_pL^p          θ(L) = 1 + θ₁L + … + θ_qL^q

L is the lag operator. AR coefficients enter the lag polynomial with a minus sign and MA coefficients with a plus sign. The outcome and every regressor are differenced, so the constant is the mean of the differenced outcome: the level mean when d = 0 and the drift per period when d = 1.

The model needs two conditions on the polynomials. The AR part must be stationary, so that shocks to the differenced series die out. The MA part must be invertible, so that the shocks can be recovered from the data. A fitted model whose inverse AR root is close to one is a sign that the series needs another difference; an inverse MA root close to one is a sign that it was differenced once too often.

The main variants:

• ARIMA for a single series. Choose d with a unit-root test (see unit-root tests), then compare a few small (p, q) orders by the correlogram and information criteria.
• Seasonal ARIMA adds multiplicative seasonal AR, MA and differencing terms at lag s (12 for monthly data, 4 for quarterly). The (0,1,1)×(0,1,1) "airline" model is the classic case.
• ARMAX, a regression with ARMA errors, when you care about the effect of x on y but the residuals are autocorrelated. Harvey and Phillips (1979) give the exact maximum-likelihood treatment.

Two estimators are common. Exact maximum likelihood evaluates the Gaussian likelihood with the presample values drawn from the stationary distribution, the likelihood a Kalman filter produces. Conditional sum of squares sets the presample values to zero. The two agree in long samples and can differ in short ones. Either way, the residuals should look like white noise. The portmanteau statistics of Box and Pierce (1970) and Ljung and Box (1978) test the first few residual autocorrelations jointly.

In OpenEconometrics

oe.arima fits ARIMA, seasonal ARIMA and ARMAX models. The call shape is oe.arima(data=df, y="y", order=(p, d, q), time="t"). With the series declared by tsset t, the model in the example below corresponds to the Stata command arima y, arima(1,1,1).

• time= is optional. It names an integer or datetime column, and the rows are sorted by it. Integer periods must be consecutive and distinct. Datetime values are ranked and treated as consecutive periods, with a recorded warning, because a calendar gap cannot be told from the sampling frequency. Without time=, the row order is the time order.
• seasonal=(P, D, Q) with period=s adds the seasonal terms (Stata's sarima(P,D,Q,s)), and x=[...] adds regressors for an ARMAX model.
• method="ml" (the default) is exact maximum likelihood; method="css" is conditional sum of squares, the workflow of Stata's condition option. CSS does not need a stationary AR part.
• covariance="opg" (the default, outer product of the gradients), "nonrobust" (observed information, Stata's vce(oim)) or "robust" (the sandwich). Tests are z tests.
• constant=False drops the constant, like Stata's noconstant.

The likelihood search starts from Hannan–Rissanen estimates (Hannan and Rissanen 1982) and, for mixed models in moderate samples, screens several starting points, because ARMA likelihoods can have more than one local maximum. The result reports the log likelihood, AIC and BIC in fit.metrics, the inverse roots of each polynomial in fit.extra["roots"], and four residual tests in fit.tests: a Wald test of the model, Ljung–Box, Jarque–Bera and a one-lag ARCH-LM test. The Ljung–Box entry also carries p_value_adjusted, with the degrees of freedom reduced by the number of ARMA parameters.

oe.forecast(fit, steps=h) (also exported as oe.arima_forecast) returns dynamic forecasts of the outcome in levels with standard errors and normal intervals. Pass future regressor values with exog= for an ARMAX model, or a longer table with data= to forecast from a later origin with the parameters held fixed. For identification, oe.corrgram prints autocorrelations, partial autocorrelations and Ljung–Box statistics, and oe.wntestq runs the portmanteau test on any column. The ARIMA family guide gives the likelihood, covariance and forecast formulas.

Supported options

Options of oe.arima

Generated from oe.capabilities() of openecon 0.3.20a1. Stata names identify the corresponding workflow, not numerical parity.

Example

The script simulates 300 periods of a series whose first difference has a drift of 0.2 and ARMA(1,1) dynamics with AR coefficient 0.5, MA coefficient 0.3 and unit shock variance. The true model is therefore ARIMA(1,1,1) with a constant of 0.2 and σ = 1. The script fits that model and forecasts four periods ahead.

import numpy as np
import pandas as pd
import openecon as oe


def simulate(n=300, seed=11):
    """A series whose first difference is ARMA(1,1) with drift 0.2: ARIMA(1,1,1)."""
    rng = np.random.default_rng(seed)
    e = rng.normal(size=n + 1)
    u = np.zeros(n)
    prev = 0.0
    for t in range(n):
        prev = 0.5 * prev + e[t + 1] + 0.3 * e[t]  # (1 - 0.5L) u_t = (1 + 0.3L) e_t
        u[t] = prev
    y = 100 + np.cumsum(0.2 + u)                    # integrate the differences
    return pd.DataFrame({"t": np.arange(1, n + 1), "y": y})


if __name__ == "__main__":
    df = simulate()
    fit = oe.arima(data=df, y="y", order=(1, 1, 1), time="t")
    print(fit.summary())
    print(oe.forecast(fit, steps=4).round(3))

With openecon 0.3.20a1 from PyPI the script prints:

ARIMA(1,1,1) regression — y
Observations: 299  |  Covariance: opg  |  Confidence: 95%

Term        Estimate  Std. error        z      P>|stat|     CI lower  CI upper
----------  --------  ----------  -------  ------------  -----------  --------
[D.y]
Intercept   0.263243    0.136652  1.92638     0.0540575  -0.00458994  0.531075
[ARMA]
ARMA:L1.ar  0.485304   0.0726929  6.67609   2.45406e-11     0.342829   0.62778
ARMA:L1.ma  0.340949   0.0791062  4.31002    1.6324e-05     0.185904  0.495995
/sigma      0.915421   0.0412938  22.1685  3.45968e-109     0.834487  0.996356

Excluded observations: 1

log_likelihood: -398.189  |  aic: 804.378  |  bic: 819.18  |  sigma: 0.915421
n_differenced: 299  |  iterations: 5
Wald chi2 test that all coefficients except the constant are zero: chi2(2) = 239.234, p = 1.125e-52
Ljung-Box Q(40) test of the residuals: chi2(40) = 35.0311, p = 0.6932
Jarque-Bera normality test of the residuals: chi2(2) = 1.07971, p = 0.5828
ARCH-LM(1) test of the residuals: chi2(1) = 0.816763, p = 0.3661
Warning: Differencing uses the first 1 observation(s) as starting values; the estimation sample has 299 observations.
   period  forecast  std_error   ci_low  ci_high
0     301   177.651      0.915  175.857  179.445
1     302   177.786      1.906  174.050  181.521
2     303   177.987      2.791  172.516  183.457
3     304   178.220      3.564  171.234  185.206

Reading the output

• The first observation is used up. With d = 1 the first row only supplies the starting level, so 299 differenced observations enter the likelihood, as the warning and "Excluded observations: 1" say.
• The dynamics are recovered well. The AR coefficient is 0.485 (95% interval 0.343 to 0.628) and the MA coefficient 0.341 (0.186 to 0.496); both intervals contain the true 0.5 and 0.3.
• The drift is imprecise. The Intercept in the D.y equation is the drift. It is 0.263 with a standard error of 0.137, so the interval from −0.005 to 0.531 contains the true 0.2 but also zero (p = 0.054). The mean of a persistent series is hard to pin down, even with 300 periods.
• σ is a little low. The estimate is 0.915 and its interval, 0.834 to 0.996, narrowly misses the true value of 1. One in twenty 95% intervals misses by chance; this is such a draw, not a defect of the estimator.
• The residuals look like white noise. Ljung–Box Q(40) = 35.03 with p = 0.69 finds no remaining autocorrelation, and the Jarque–Bera (p = 0.58) and ARCH-LM (p = 0.37) tests find no departure from normality or volatility clustering.
• Forecast uncertainty grows quickly. The one-step standard error, 0.915, equals σ. Because the series is integrated, shocks are permanent and the standard error reaches 3.564 after four steps. These intervals do not include parameter uncertainty.

Limitations

• oe.arima accepts no weights.
• The common oe.predict and oe.margins adapters do not cover oe.arima; use its own result tables and family helpers.
• One series per fit; there is no panel ARIMA.
• The series must be regularly spaced without holes. missing="drop" removes incomplete rows only at the start or end; a missing value inside the series raises time_gaps, and integer time values must increase by one (there is no delta option).
• AR and MA lags are consecutive (1 to p, 1 to q); subset lag lists such as Stata's ar(1 4) are not available.
• method="ml" models the differenced data as a stationary process and requires a stationary AR part; use method="css" or more differencing for a series with a unit root. The exact diffuse likelihood of the undifferenced series is not implemented.
• Forecast standard errors ignore parameter uncertainty, and with a datetime time column the forecast period is the step number.
• No numerical parity with Stata is claimed: the inventory's stata_parity_validated flag is False.

References

• Box, G. E. P. and Pierce, D. A. (1970). Distribution of residual autocorrelations in autoregressive-integrated moving average time series models. Journal of the American Statistical Association, 65(332), 1509–1526. doi:10.1080/01621459.1970.10481180
• Hannan, E. J. and Rissanen, J. (1982). Recursive estimation of mixed autoregressive-moving average order. Biometrika, 69(1), 81–94. doi:10.1093/biomet/69.1.81
• Harvey, A. C. and Phillips, G. D. A. (1979). Maximum likelihood estimation of regression models with autoregressive-moving average disturbances. Biometrika, 66(1), 49–58. doi:10.1093/biomet/66.1.49
• Ljung, G. M. and Box, G. E. P. (1978). On a measure of lack of fit in time series models. Biometrika, 65(2), 297–303. doi:10.1093/biomet/65.2.297