Mastering ARIMA: The Complete Guide to Time Series Forecasting

I took a common non-seasonal ARIMA recipe and ran it on a classic forecasting dataset: monthly airline passengers from 1949 to 1960. The model looked fine on a standard residual check. Then I gave the same data to a baseline that takes one line of code, and the ARIMA model lost on the 12-month holdout.

That result is the point of this guide. In this experiment the model did not fail because ARIMA mathematics is weak. It failed because the workflow around it was weak: nobody checked whether the data suited a non-seasonal model, whether the diagnostic was measuring what it appeared to measure, or whether the model beat something trivial.

Most tutorials teach you to fit the model. This one teaches you to choose it, test it, trust it, and walk away from it when needed. You will get a three-gate framework, working Python code, and an experiment you can reproduce yourself, with the data source, software versions and raw results listed.

What ARIMA Actually Does

ARIMA (Autoregressive Integrated Moving Average) is a statistical model that forecasts a single time series from its own past values and past forecast errors, after removing trend through differencing.

It answers one narrow question: does this series have memory, and can that memory be modeled? Box and Jenkins formalized the approach in 1970, and it is still taught and used today. The model has three parts.

• AR (autoregressive): today’s value depends on recent values.
• I (integrated): the series is differenced to remove trend.
• MA (moving average): today’s value depends on recent forecast errors.

Picture monthly sales for a coffee brand. Last month’s sales pull this month’s number (AR). A cold snap two months ago caused a surprise spike, and the model learns to correct for that kind of shock (MA). Sales also grow year over year, so you difference the series first to get a stable signal (I).

Notice the MA part has nothing to do with the moving average on a stock chart. It models how past mistakes echo into the present.

The AR Component in One Formula

Here is the AR(p) model. The parameter p says how many past values the model looks at.

A φ near 0.9 means yesterday carries strongly into today. A φ near 0.2 means the series forgets quickly. Try it below with an AR(1) model.

ARIMA combines momentum (AR), shock correction (MA) and trend removal (I). It is simple in concept, and the difficulty lies in the workflow.

Gate 1: The ARIMA Readiness Test

Before opening a notebook, ask whether this series is a candidate at all. Five questions settle most cases.

1. Is it one variable? ARIMA models a single series. If promotions, weather or policy drive your numbers, you probably want ARIMAX, which adds external regressors.
2. Is the pattern mostly linear? Heavy non-linearity and explosive volatility break the model.
3. Is the spacing regular? Irregular timestamps need preprocessing.
4. Is the horizon short to medium? Forecasts flatten as the horizon grows. More on that below.
5. Are there structural breaks? A merger, a regulation or a pandemic splits your data into two different processes.

The statsmodels ARIMA class handles seasonal terms and exogenous regressors in the same interface, so upgrading is cheap. See the statsmodels ARIMA documentation for seasonal_order and exog.

The ARIMA Stoplight

Use this editorial framework to sort your situation in two minutes. It is a heuristic assembled from the sources in this guide, not an empirically validated rule, so treat the thresholds as starting points.

Stationarity: Test It, Don’t Eyeball It

A stationary series has a constant mean, stable variance and consistent autocorrelation over time, with no trend and no growing swings.

Two tests are standard, and they ask opposite questions.

• ADF test: the null hypothesis is a unit root, meaning non-stationary. A small p-value rejects it. The statsmodels adfuller documentation explains the options and notes that in version 0.15.0 the function still returns a plain tuple by default and warns about a future change, which is why the code below passes result_object=True. The original method comes from Dickey and Fuller’s 1979 paper.
• KPSS test: the null hypothesis is stationarity, either around a constant or around a trend, as the statsmodels kpssdocumentation describes. If ADF and KPSS agree, you can trust the verdict. If they disagree, investigate before modeling. We cover this in our practical guide to the KPSS test in Python.

How Much Data Do You Need?

Nobody can give you a universal number. A common practitioner rule of thumb is about 50 observations as a floor, 100 or more for comfort, and several full cycles if seasonality is present. Treat it as a heuristic, not a law.

What is well established is the small-sample correction: Hurvich and Tsai’s paper on regression and time series model selection in small samples derives a bias-corrected AIC (AICc) for that situation, and it penalizes complexity more when data is scarce.

Choosing p, d and q Without Guesswork

Most people obsess over p and q. Get d right first, because everything else depends on it.

How to choose d, step by step:

1. Plot the series and look for trend or growing variance.
2. If variance grows, apply a log or Box-Cox transformation first.
3. Run ADF (and KPSS) on the series.
4. If it is non-stationary, difference once and test again.
5. Rarely go beyond d = 2 (a practitioner rule of thumb). If you seem to need more, suspect a structural break or non-linearity.

A common workflow is to transform first, then difference. This is my recommendation rather than a documented rule: a log transform stabilizes variance, and differencing the logged series then removes trend. The SciPy boxcox documentation covers the lambda parameter, where a value of zero corresponds to a log transform. Box-Cox also requires strictly positive data; in my test SciPy 1.17.1 raised ValueError: Data must be positive. for zeros and negative values.

For p and q, the classic pattern guide is: PACF cuts off after lag p for a pure AR process, and ACF cuts off after lag q for a pure MA process. Penn State’s STAT 510 lesson on identifying ARIMA models shows how these patterns guide a first guess. Real plots are messy. In practice, fit several small (p, q) combinations and compare them, but only under the rule you will meet in the AIC section: same data, same differencing.

One more point most guides skip: what the constant does. Rob Hyndman’s note on constants in ARIMA models lists the long-run behavior. With d = 1 and no constant, forecasts head to a non-zero constant set by the last few observations. With d = 1 and a constant, they follow a straight line whose slope is the mean of the differenced data. With d = 0 and a constant, they go to the series mean.

For seasonal models the same rules apply with d replaced by d + D. This is why ARIMA forecasts flatten, and it matters in Python: the statsmodels ARIMA documentation says the default is a constant for models without integration and no trend for models with integration. In my test, a non-seasonal ARIMA(1,1,1) fitted to the 132 training observations of the airline data gave a flat 60-step forecast by default, while trend='t' produced a rising line.

Fix d first, transform before differencing, and treat ACF/PACF as hints rather than answers.

The Box-Jenkins Workflow

The Box-Jenkins method follows four iterative stages:

1. Identification: plot, test stationarity, transform, difference, inspect ACF and PACF.
2. Estimation: fit candidate models by maximum likelihood.
3. Diagnostic checking: examine residuals for leftover structure.
4. Forecasting: produce point forecasts and prediction intervals.

Stage 3 is the one people rush, and it is where this guide’s biggest lesson lives. If diagnostics fail, go back to Stage 1. Tuning parameters on a structurally wrong model does not rescue it.

A Worked Example: Where the Standard Recipe Breaks

Let’s use the classic airline passengers dataset: 144 monthly values from January 1949 to December 1960. It has strong trend, growing seasonal swings and a clear annual cycle. I ran everything below on October 16, 2026, holding out the final 12 months. The table lists every condition so you can reproduce the result or challenge it.

If you are choosing tools for this kind of work, our guide to choosing Python data libraries covers where statsmodels fits next to pandas and others. If your data loading feels slow, see why pandas can be slow and how to fix it.

Step 1: Load, Split and Test Stationarity

Python: Load, split, test stationarity:

import numpy as np import pandas as pd from statsmodels.tsa.stattools import adfuller, kpss
df = pd.read_csv("airline.csv", parse_dates=["Month"], index_col="Month") df.index.freq = "MS" y = df["Passengers"] y_train, y_test = y.iloc[:-12], y.iloc[-12:] # 132 training months, 12 test months
adf = adfuller(y_train, result_object=True) kp = kpss(y_train, regression="c", nlags="auto", result_object=True) print("Training observations:", len(y_train)) print("ADF p-value:", round(adf.pvalue, 4)) print("KPSS p-value:", kp.pvalue)

Result: the ADF p-value is 0.9929, so a unit root cannot be rejected. The KPSS p-value is 0.01, which is a boundary value: the statsmodels kpss documentation says the p-value is interpolated from a table and a boundary point is returned when the statistic falls outside it, and statsmodels prints an InterpolationWarning to say so. Read it as “p is at most 0.01”. Both tests agree that the raw series is not stationary.

Step 2: Stabilize Variance, Then Difference

Python: Box-Cox on training data, then test the differenced series:

from scipy.stats import boxcox
bc_values, lam = boxcox(y_train) bc = pd.Series(bc_values, index=y_train.index) print("Box-Cox lambda:", round(lam, 3)) print("ADF after 1 difference :", round(adfuller(bc.diff().dropna(), result_object=True).pvalue, 4)) print("ADF after 1 + seasonal (12) :", round(adfuller(bc.diff().diff(12).dropna(), result_object=True).pvalue, 4))

Lambda comes out at 0.129. After one ordinary difference, the ADF p-value is 0.0255, so the test passes at the 5% level. After adding a seasonal difference at lag 12 it is 0.0012.

Notice what this does and does not tell you. A pass after one difference would let you proceed with a plain ARIMA. But ADF asks only about a unit root; it cannot tell you whether a non-seasonal model is adequate. In my interpretation, that is exactly the trap in this example: the stationarity test is satisfied while a seasonal pattern remains that a non-seasonal model cannot absorb.

Step 3: Inspect ACF and PACF

Python: ACF and PACF of the stationary series:

import matplotlib.pyplot as plt from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
stationary = bc.diff().diff(12).dropna() fig, ax = plt.subplots(1, 2, figsize=(12, 4)) plot_acf(stationary, lags=24, ax=ax[0]) plot_pacf(stationary, lags=24, ax=ax[1]) plt.tight_layout() plt.show()

Step 4: Fit Two Models and Compare Them Against Baselines

Here is the part tutorials skip. We fit the non-seasonal ARIMA(2,1,2), a seasonal model, and two baselines that need no fitting at all.

Python: Fit ARIMA and SARIMA, forecast 12 months, compare with baselines:

from statsmodels.tsa.arima.model import ARIMA from scipy.special import inv_boxcox
def mape(actual, forecast): return float(np.mean(np.abs((actual - forecast) / actual)) * 100)
arima = ARIMA(bc, order=(2, 1, 2)).fit() sarima = ARIMA(bc, order=(1, 1, 1), seasonal_order=(0, 1, 1, 12)).fit()
forecasts = { "ARIMA(2,1,2)": inv_boxcox(arima.forecast(12), lam).values, "Naive (last value)": np.repeat(y_train.iloc[-1], 12), "Seasonal naive": y_train.iloc[-12:].values, # same month last year "SARIMA(1,1,1)(0,1,1,12)": inv_boxcox(sarima.forecast(12), lam).values, } for name, f in forecasts.items(): print(f"{name:26s} MAPE: {mape(y_test.values, f):.2f}%")

MAPE is the mean absolute percentage error: the average size of the miss as a percentage of the actual value, lower is better. It is fine here because every actual value is far from zero. FPP3’s discussion of accuracy measures explains why percentage errors break down near zero and recommends scaled errors such as MASE when you compare across series.

Read that table twice. On this holdout, the non-seasonal ARIMA(2,1,2) was worse than repeating last year’s numbers, and the seasonal version cut the error by roughly three quarters. This is also the logic behind the baseline gate. FPP3 states the principle plainly: forecasting methods you develop should be compared with simple benchmarks, and a method that does not beat them is not worth considering.

Before vs. after:

• Before: a plausible-looking ARIMA with a tidy forecast line, no baseline, and one number nobody questioned.
• After: a baseline gate that exposes the weak specification, and a seasonal specification that earns its place.

Here is the output you should see if you run Steps 1 to 6 on the same data and versions (the InterpolationWarning from KPSS appears on stderr):

Expected output (stdout):

Training observations: 132 ADF p-value: 0.9929 KPSS p-value: 0.01 Box-Cox lambda: 0.129 ADF after 1 difference : 0.0255 ADF after 1 + seasonal (12) : 0.0012 ARIMA(2,1,2) MAPE: 12.73% Naive (last value) MAPE: 14.25% Seasonal naive MAPE: 9.99% SARIMA(1,1,1)(0,1,1,12) MAPE: 2.97% loglikelihood_burn ARIMA: 1 SARIMA: 13 First ARIMA residual: 6.494 | first training value: 6.494 ARIMA Ljung-Box p, all residuals : 0.9996 ARIMA Ljung-Box p, burn-in removed : 4.7467168640214974e-35 SARIMA Ljung-Box p, all residuals : 0.25 SARIMA Ljung-Box p, burn-in removed : 0.3941 windows: 9 mean MAPE ARIMA: 11.95 | SARIMA: 4.38 | seasonal naive: 9.42 SARIMA beats seasonal naive in 8 of 9 windows ARIMA beats seasonal naive in 2 of 9 windows

Gate 2: The Diagnostic That Lies

Here is the finding I did not expect. Run the residual check on the ARIMA(2,1,2) model in two ways and compare.

Python: Residual diagnostics, naive vs. honest:

from statsmodels.stats.diagnostic import acorr_ljungbox
def ljung_box_p(fit, n_arma_params, trim=True): resid = fit.resid.iloc[fit.loglikelihood_burn:] if trim else fit.resid lags = min(24, len(resid) // 5) # 2 x period for monthly data, capped at a fifth of the sample return float(acorr_ljungbox(resid, lags=[lags], model_df=n_arma_params).iloc[0, 1])
print("loglikelihood_burn ARIMA:", arima.loglikelihood_burn, " SARIMA:", sarima.loglikelihood_burn) print("First ARIMA residual:", round(arima.resid.iloc[0], 3), "| first training value:", round(bc.iloc[0], 3)) print("ARIMA Ljung-Box p, all residuals :", round(ljung_box_p(arima, 4, trim=False), 4)) print("ARIMA Ljung-Box p, burn-in removed :", ljung_box_p(arima, 4)) print("SARIMA Ljung-Box p, all residuals :", round(ljung_box_p(sarima, 3, trim=False), 4)) print("SARIMA Ljung-Box p, burn-in removed :", round(ljung_box_p(sarima, 3), 4)) sarima.plot_diagnostics(figsize=(12, 7)) plt.tight_layout() plt.show()

On my run the numbers were:

• ARIMA(2,1,2), all residuals: Ljung-Box p-value of 0.9996. A perfect pass.
• ARIMA(2,1,2), burn-in removed: p-value of about 5 × 10⁻³⁵. A decisive fail.
• SARIMA, all residuals: p-value of 0.25. SARIMA, burn-in removed: 0.3941. Both pass, so this trap mattered mostly for the non-seasonal model.

What happened? With differencing, the first residual is not an ordinary forecast error. In my run it was 6.494, identical to the first Box-Cox value in the training series, and loglikelihood_burn was 1 for the ARIMA model and 13 for the seasonal one (the regular difference plus the seasonal difference at lag 12).

My interpretation is that these are start-up residuals, and that one huge artificial value hides the autocorrelation in the rest. I verified the behavior in statsmodels 0.15.0 by inspection, but I did not find loglikelihood_burn described in the documentation pages I opened, so treat it as observed behavior and re-check it after upgrading.

Two more details in that code come from FPP3’s section on residual diagnostics. The suggested number of lags is 10 for non-seasonal data and twice the seasonal period for seasonal data, capped at a fifth of the sample, and for ARIMA models the degrees of freedom should be adjusted. That is what model_df does in the acorr_ljungbox function. The same section adds a caution that matters here: checking residuals is important, but it is not a good way to select a forecasting method.

Rule to remember: after differencing, remove the burn-in residuals, adjust the degrees of freedom, and never treat a clean residual check as proof of forecast quality.

When Things Go Wrong: Failure Modes I Hit or Verified

The AIC Trap

AIC (Akaike Information Criterion) scores a model by rewarding fit and penalizing the number of parameters. Lower is better, but only among comparable models.

There are two traps here, and the second is the one almost nobody mentions.

1. Trap 1: AIC measures in-sample fit. It cannot tell you how a model behaves on a different time period. Use it to shortlist candidates, then validate out of sample. The small-sample correction (AICc) comes from Hurvich and Tsai’s work on model selection in small samples.
2. Trap 2: AIC is invalid across different differencing orders. Hyndman and Athanasopoulos note that comparing information criteria is only valid for ARIMA models with the same orders of differencing. The reason is that differencing changes the data the likelihood is computed on.

My experiment makes the danger concrete. The ARIMA(2,1,2) model had an AIC of −40.5 and the seasonal model had −271.2 (AICc −40.0 and −270.9). The gap looks enormous, and the seasonal model is indeed better on the holdout, but the two numbers are not on the same scale because the differencing differs. The verdict came from out-of-sample error, not from the AIC gap.

A small grid search shows the first trap. Among 16 non-seasonal candidates (p and q from 0 to 3, d = 1), ARIMA(3,1,3) had the lowest AIC at −66.9, but its fit raised a ConvergenceWarning, and its holdout MAPE was 8.48%, still behind the seasonal model. Its burn-in-removed residuals failed the Ljung-Box test as well. The next best candidates, ARIMA(3,1,2) at −44.5 and ARIMA(2,1,1) at −43.3, fit without warnings. A low AIC from a fit that did not converge is not a result.

Use the calculator below. It computes AIC and AICc from the formulas above and refuses to give a verdict when the differencing differs. Its default input values are illustrative placeholders, not results from the experiment.

The thresholds in the calculator (2 and 6) are common rules of thumb rather than hard laws. Treat them as prompts to look closer.

AIC shortlists. Out-of-sample validation decides.

Gate 3: Rolling-Origin Validation Against a Baseline

One train/test split tells you how a model did in one window. Time series change, so you want to know how it does across many.

How rolling-origin validation works, step by step:

1. Train on the first T observations.
2. Forecast the next h steps and record the error.
3. Move the origin forward and retrain on a longer history.
4. Repeat until you run out of data.
5. Look at the whole error distribution, not just the mean.

FPP3’s section on time series cross-validation calls this evaluation on a rolling forecasting origin, notes that the earliest observations are not used as test sets, and describes how to adapt it to multi-step errors. It is also the habit that separated serious M5 entrants: according to the competition paper, the top four methods and most of the top 50 validated on at least the last four 28-day windows.

Python: Rolling-origin validation (expanding window):

def rolling_mape(order, seasonal=(0, 0, 0, 0), start=84, horizon=12, step=6): scores = [] for t in range(start, len(y) - horizon + 1, step): b, l = boxcox(y.iloc[:t]) # lambda from the training window only fit = ARIMA(pd.Series(b, index=y.index[:t]), order=order, seasonal_order=seasonal).fit() pred = inv_boxcox(fit.forecast(horizon), l) scores.append(mape(y.iloc[t:t + horizon].values, pred)) return scores
def rolling_seasonal_naive(start=84, horizon=12, step=6): return [mape(y.iloc[t:t + horizon].values, y.iloc[t - 12:t].values) for t in range(start, len(y) - horizon + 1, step)]
a = rolling_mape((2, 1, 2)) s = rolling_mape((1, 1, 1), (0, 1, 1, 12)) n = rolling_seasonal_naive() print("windows:", len(a)) print("mean MAPE ARIMA:", round(np.mean(a), 2), "| SARIMA:", round(np.mean(s), 2), "| seasonal naive:", round(np.mean(n), 2)) print("SARIMA beats seasonal naive in", sum(x < z for x, z in zip(s, n)), "of", len(s), "windows") print("ARIMA beats seasonal naive in", sum(x < z for x, z in zip(a, n)), "of", len(a), "windows")