Sina Bolouri
Back to Science articles

Introduction to Time Series

ARIMA, SARIMA, and Multi-Step Forecasting

How differencing, seasonal structure, and innovation recursions turn persistent time series into forecasts with explicit uncertainty.

Article 6Time-Series Forecastingintermediate23-29 minutespartially-verified
time seriesforecastingdata sciencearimasarimaunit rootsadf testseasonal differencingmulti-step forecastingforecast intervals

Introduction: The Problem Before the Method

Many important series are not stationary in their observed levels. Revenue grows, water levels drift, traffic follows recurring weekly cycles, and equipment features can change as wear accumulates. A stationary ARMA model applied directly to such data may confuse persistent level movement with short-memory dependence.

ARIMA models address this problem by separating two operations:

  1. Differencing transforms the observed level into a series whose dependence may be stationary.
  2. ARMA modeling explains the remaining stationary dynamics.

Seasonal ARIMA extends the same reasoning to repeated cycles. Forecasting then reverses the differencing operation so predictions return to the original scale.

This framework is powerful, but it creates several responsibilities. The analyst must justify the differencing order, distinguish stochastic trend from deterministic trend, avoid over-differencing, model seasonality at the correct period, and carry increasing uncertainty through the inverse transformation.

The Core Idea

An ARIMA((p,d,q)) model is

ϕ(B)(1B)dXt=θ(B)Wt.\phi(B)(1-B)^dX_t=\theta(B)W_t.

The differenced series

Yt=(1B)dXtY_t=(1-B)^dX_t

is modeled as a stationary ARMA((p,q)) process.

The notation separates:

  • (p): nonseasonal autoregressive order,
  • (d): number of ordinary differences,
  • (q): nonseasonal moving-average order.

For seasonal period (s), a multiplicative SARIMA model is commonly written

Φ(Bs)ϕ(B)(1B)d(1Bs)DXt=Θ(Bs)θ(B)Wt,\Phi(B^s)\phi(B) (1-B)^d(1-B^s)^D X_t = \Theta(B^s)\theta(B)W_t,

and denoted

ARIMA(p,d,q)(P,D,Q)s.\text{ARIMA}(p,d,q)(P,D,Q)_s.

The seasonal polynomial terms operate at lags (s,2s,\ldots), while ordinary terms operate at consecutive lags.

Why the Concept Exists

A random walk demonstrates why an ARMA model on levels can fail:

Xt=Xt1+Wt.X_t=X_{t-1}+W_t.

The process has a unit root because its AR polynomial is (1-B). Shocks accumulate permanently:

Xt=X0+j=1tWj.X_t=X_0+\sum_{j=1}^{t}W_j.

Its variance grows with (t), so the level is nonstationary. First differencing produces

(1B)Xt=Wt,(1-B)X_t=W_t,

which is stationary.

ARIMA treats this accumulated innovation mechanism explicitly. Instead of forcing a stationary model onto the level, it models the changes and then integrates forecasts back into levels.

Seasonal differencing solves an analogous problem when seasonal effects persist from one cycle to the next:

(1Bs)Xt=XtXts.(1-B^s)X_t=X_t-X_{t-s}.

Background and Prerequisites

Deterministic trend versus stochastic trend

A deterministic linear trend model is

Xt=a+bt+Yt,X_t=a+bt+Y_t,

where (Y_t) is stationary. Removing (a+bt) yields a stationary residual if the trend form is correct.

A random walk with drift is

Xt=c+Xt1+Wt.X_t=c+X_{t-1}+W_t.

Its first difference is

ΔXt=c+Wt.\Delta X_t=c+W_t.

The deterministic-trend model returns toward its fitted trend after a shock. The random walk with drift permanently incorporates the shock into its future level. This distinction changes long-horizon uncertainty and interpretation.

Unit roots

For an AR((p)) polynomial (\phi(z)), a root at (z=1) corresponds to a factor (1-B). A seasonal unit root may occur at roots associated with (1-B^s). Near-unit roots can resemble exact unit roots in finite samples, making diagnosis uncertain.

Assumptions

Differenced dynamics are stable

After applying the selected ordinary and seasonal differences, the remaining process should be reasonably stationary.

The seasonal period is meaningful

A monthly series may have (s=12), but this should reflect an actual yearly cycle. Multiple or noninteger seasonalities may require Fourier terms, dynamic harmonic regression, or more flexible models.

Innovations are sufficiently well behaved

Classical intervals often assume uncorrelated, constant-variance, approximately Gaussian innovations. Heavy tails or changing variance affect calibration.

Future structure resembles the fitted regime

ARIMA extrapolates differenced dynamics and deterministic terms. Structural breaks, interventions, and changing seasonality can invalidate forecasts.

How It Works

Choosing (d): ordinary integration

Plot the series and inspect whether shocks appear persistent and whether the ACF remains very high across many lags. Compare the original and differenced series. A first difference is common, but (d=2) is appropriate only when the first difference still has a changing level consistent with second-order integration.

Over-differencing often produces strong negative lag-one autocorrelation and unnecessarily large forecast variance.

Augmented Dickey-Fuller testing

The Dickey-Fuller framework tests a null hypothesis of a unit root. In an augmented form, a common regression is

ΔXt=α+βt+γXt1+i=1kδiΔXti+εt.\Delta X_t = \alpha+\beta t+\gamma X_{t-1} + \sum_{i=1}^{k}\delta_i\Delta X_{t-i} +\varepsilon_t.

The unit-root null is expressed through (\gamma=0), with critical values that are not the ordinary Student-(t) values. The augmentation lags help account for serial correlation in the regression errors (Dickey & Fuller, 1979; Said & Dickey, 1984).

The deterministic terms—none, intercept, or intercept plus trend—and lag order materially affect the test. Failure to reject may reflect a unit root, low power, short data, or a near-unit stationary process. The test should be combined with plots, domain knowledge, and model diagnostics.

Choosing seasonal structure

A seasonal spike in the ACF at (s), (2s), and later multiples can indicate unresolved seasonality. Seasonal differencing addresses evolving seasonal levels. Seasonal AR or MA terms address dependence at seasonal lags after differencing.

A useful order-identification sequence is:

  1. stabilize variance if needed;
  2. choose (d) and (D);
  3. inspect the ACF and PACF of the differenced series;
  4. propose low-order nonseasonal and seasonal terms;
  5. fit candidates;
  6. compare diagnostics and AICc;
  7. evaluate forecasts.

Multi-step forecasting

For a causal model with innovation representation

Xt=j=0ψjWtj,X_t=\sum_{j=0}^{\infty}\psi_jW_{t-j},

the (h)-step forecast error contains future innovations:

en(h)=Xn+hX^n+hn=j=0h1ψjWn+hj.e_n(h) = X_{n+h}-\widehat{X}_{n+h|n} = \sum_{j=0}^{h-1}\psi_jW_{n+h-j}.

Therefore,

Var(en(h))=σ2j=0h1ψj2.\operatorname{Var}(e_n(h)) = \sigma^2 \sum_{j=0}^{h-1}\psi_j^2.

The sum gains nonnegative terms as (h) increases, explaining why forecast intervals usually widen.

The Practical Procedure

1. Define the forecasting target and horizon

A model for monthly levels may be inappropriate if the decision concerns daily changes. The horizon influences which seasonal cycles and error patterns matter.

2. Audit the time index and interventions

Correct missing timestamps, duplicate periods, calendar anomalies, and known changes before interpreting unit-root behavior.

3. Transform variance before differencing

For positive data whose variability rises with level, a logarithm or Box-Cox transformation may improve additivity and interval behavior.

4. Select the minimum plausible (d) and (D)

Use plots, differenced plots, ACF behavior, and unit-root evidence. Avoid mechanical repeated differencing.

5. Identify low-order AR and MA terms

Inspect ACF/PACF after all selected differences. Seasonal and nonseasonal terms should be considered jointly.

6. Fit, diagnose, and compare

Check convergence, roots, residuals, Ljung-Box results, AICc, and coefficient uncertainty.

7. Produce horizon-specific forecasts

Use the fitted recursion and model-implied forecast-error variance. Report the transformation and interval assumptions.

8. Reconstruct the original scale

Undo differencing using known historical values. For nonlinear transformations, distinguish the median forecast obtained by simple back-transformation from a bias-adjusted mean forecast.

9. Backtest across multiple origins

Evaluate point errors and interval coverage at the horizons the system will actually use.

Mathematical or Technical Foundation

AR(1) multi-step forecast

For

Xt=ϕXt1+Wt,ϕ<1,X_t=\phi X_{t-1}+W_t, \qquad |\phi|<1,

the (h)-step forecast is

X^n+hn=ϕhXn.\widehat{X}_{n+h|n}=\phi^hX_n.

The forecast-error variance is

σh2=σ2j=0h1ϕ2j=σ21ϕ2h1ϕ2.\sigma_h^2 = \sigma^2 \sum_{j=0}^{h-1}\phi^{2j} = \sigma^2 \frac{1-\phi^{2h}}{1-\phi^2}.

As (h\to\infty), the forecast approaches the stationary mean and the error variance approaches the unconditional variance.

Random walk forecast

For

Xt=Xt1+Wt,X_t=X_{t-1}+W_t, X^n+hn=Xn,\widehat{X}_{n+h|n}=X_n,

and

Var(en(h))=hσ2.\operatorname{Var}(e_n(h))=h\sigma^2.

The point forecast remains at the latest level, but uncertainty grows without bound.

Random walk with drift

For

Xt=c+Xt1+Wt,X_t=c+X_{t-1}+W_t, X^n+hn=Xn+hc,\widehat{X}_{n+h|n}=X_n+hc,

with the same (h\sigma^2) innovation contribution when parameters are treated as known. Estimated-drift uncertainty adds further variance in practice.

AR(2) recursion

For

Xt=ϕ1Xt1+ϕ2Xt2+Wt,X_t=\phi_1X_{t-1}+\phi_2X_{t-2}+W_t, X^n+1n=ϕ1Xn+ϕ2Xn1,\widehat{X}_{n+1|n} = \phi_1X_n+\phi_2X_{n-1}, X^n+2n=ϕ1X^n+1n+ϕ2Xn,\widehat{X}_{n+2|n} = \phi_1\widehat{X}_{n+1|n}+\phi_2X_n,

and later horizons recursively replace unavailable future observations with their forecasts.

Worked Example

Source-derived example: AirPassengers

The monthly AirPassengers series has a rising level and seasonal amplitude. The uploaded material uses a logarithm, ordinary differencing, seasonal differencing, and a seasonal ARIMA model.

x <- log(AirPassengers)
y <- diff(diff(x, lag = 12), differences = 1)

acf(y)
pacf(y)

fit <- arima(
  x,
  order = c(0, 1, 1),
  seasonal = list(order = c(0, 1, 1), period = 12)
)

tsdiag(fit)
forecast <- predict(fit, n.ahead = 24)

This commonly used “airline” structure is a candidate, not a universal template for monthly data. Residual diagnostics, parameter estimates, and forecast evaluation determine whether it is suitable.

Source-derived example: BJsales

The course material also uses BJsales to illustrate differencing and forecast reconstruction. A compact workflow is:

x <- BJsales
dx <- diff(x)

fit <- arima(dx, order = c(1, 0, 1))
pred_dx <- predict(fit, n.ahead = 12)

pred_level <- tail(x, 1) + cumsum(pred_dx$pred)

The cumulative sum reverses first differencing. Forecast uncertainty must also be accumulated; simply adding transformed-scale standard errors is not correct.

Original explanatory example: condition-feature forecasting

Suppose an hourly vibration feature has:

  • a gradual persistent level,
  • a 24-hour operational cycle,
  • and short-run residual autocorrelation.

A candidate model might use a log transform, seasonal difference (1-B^{24}), ordinary difference (1-B), and low-order ARMA terms. Before deployment:

  1. verify that the absolute feature level remains available for safety thresholds;
  2. backtest at 1-, 6-, and 24-hour horizons;
  3. inspect interval coverage during both normal and high-load periods;
  4. retrain or include interventions after maintenance events.

The statistical forecast should complement, not replace, engineering limits and fault logic.

Interpreting the Results

The differencing orders describe how many persistent ordinary and seasonal components were removed. They should not be interpreted as physical laws.

AR and MA coefficients describe dependence in the transformed series, not directly in the original level. A negative MA coefficient after differencing can partly represent short-run correction introduced by the differencing operation.

Forecast intervals describe uncertainty conditional on the fitted model, estimated parameters, assumed innovation distribution, and future regime. They do not include every source of operational uncertainty unless the model explicitly represents it.

Real-World Applications

ARIMA is useful for demand, workload, inventory, environmental measurements, and transformed sensor features when differencing yields stable linear dynamics. SARIMA is useful when a dominant fixed seasonal period remains important.

In software capacity planning, seasonal ARIMA can model periodic request volume after logarithmic transformation. In predictive maintenance, it can forecast stationary changes or deviations, but absolute thresholds and maintenance interventions must remain explicit.

Common Misunderstandings

“The ADF test proves a series is stationary or nonstationary”

It tests a specific null under a selected regression and lag structure. Its power and conclusions depend on those choices.

“A deterministic trend and random walk with drift are interchangeable”

A deterministic-trend process returns toward a fixed path after shocks; a random walk with drift permanently accumulates them.

“Seasonal differencing and a seasonal AR term do the same thing”

Seasonal differencing removes a seasonal unit-root-like component. A seasonal AR term models stationary dependence at seasonal lags.

“Forecast intervals widen only because the model is uncertain”

Even with known parameters, future innovations accumulate with horizon. Parameter uncertainty adds another source.

“Exponentiating a log forecast always gives the mean forecast”

Direct exponentiation generally yields a median under a lognormal predictive distribution. Mean forecasts require a variance-based bias adjustment under the relevant assumptions.

Limitations and Failure Modes

ARIMA can over-difference near-stationary data, miss changing seasonality, and respond poorly to structural breaks. Parameter estimates near unit or invertibility boundaries can be unstable. Long-horizon forecasts may become dominated by drift or revert to a simple long-run pattern.

A fixed seasonal period is unsuitable for multiple, evolving, or calendar-driven seasonalities without additional regressors or harmonic terms. Missing data and irregular spacing require state-space handling or preprocessing rather than blind application of standard formulas.

Alternatives and Trade-Offs

Regression with deterministic trend and seasonal features is preferable when those structures are stable and interpretable. Dynamic regression combines predictors with ARIMA errors. Exponential smoothing models evolving level, trend, and seasonality directly. Structural state-space models represent components probabilistically. Prophet-like decompositions, TBATS, and dynamic harmonic regression may handle multiple or complex seasonalities, though each introduces its own assumptions and complexity.

The choice depends on whether nonstationarity is better represented by differencing a level or by evolving latent components.

Connection to Broader Topics

ARIMA is a special case of linear state-space modeling, enabling Kalman-filter treatment of missing data and recursive likelihood. Unit-root analysis connects time-series forecasting with econometrics. Seasonal differencing relates to seasonal filters and frequency-domain zeros.

Forecast-error variance leads naturally to probabilistic scoring, interval calibration, and decision-making under uncertainty.

Connection to Portfolio or Learning

A portfolio forecasting system could expose the full reconstruction chain:

raw observations
→ transformation
→ ordinary/seasonal differences
→ fitted ARMA dynamics
→ transformed forecasts
→ reconstructed levels
→ interval and backtest report

For KineticNode or condition monitoring, maintenance events should be stored as interventions rather than left for a differencing operation to absorb. This preserves operational meaning and prevents artificial forecast errors at known change points.

Key Takeaways

  1. ARIMA models stationary dynamics in a differenced representation and then reconstructs forecasts on the original scale.
  2. Ordinary and seasonal differencing address different persistent structures and should be used at the minimum defensible orders.
  3. The augmented Dickey-Fuller test supplies evidence about a unit-root null but cannot replace plots, context, and diagnostic comparison.
  4. SARIMA combines nonseasonal and seasonal AR, differencing, and MA polynomials at a specified period.
  5. Multi-step forecast uncertainty grows because each future horizon introduces additional unknown innovations.
  6. A random walk forecast remains at the latest level while its variance grows linearly with horizon.
  7. Inverse differencing and nonlinear back-transformation require explicit handling; point and interval forecasts cannot be reconstructed casually.
  8. Forecast validity remains conditional on a stable future regime and must be checked with rolling-origin evaluation.

Review Questions

  1. How does a stochastic trend differ from a deterministic trend after a shock?
  2. What null hypothesis is tested in an augmented Dickey-Fuller regression?
  3. Why can over-differencing create negative lag-one autocorrelation?
  4. What roles do (D) and (P) play differently in a seasonal ARIMA model?
  5. Why does the AR(1) forecast converge to the mean while a random-walk forecast does not?
  6. What must be retained to reconstruct forecasts after first and seasonal differencing?
  7. Why can a log-scale median forecast differ from the original-scale mean forecast?

Further-Learning Path

Stationarity and ARMA models
→ provide the stationary dependence foundation.

Unit roots, ARIMA, and SARIMA
→ extend the model to persistent level and seasonal structure.

State-space representation and Kalman filtering
→ support missing data, recursive estimation, and component interpretation.

Rolling-origin and probabilistic forecast evaluation
→ assess horizon-specific point and interval performance.

Intervention and dynamic-regression models
→ incorporate known events and external predictors.

  1. Diagnosing and Selecting Time-Series Modelsprerequisite. Provides adequacy and model-comparison tools.
  2. Exponential Smoothing as a State-Space Modelcomparison. Represents evolving components without explicit ARIMA differencing choices.
  3. Unit Roots and the Augmented Dickey-Fuller Testdeeper theory. Develops deterministic terms, lag selection, power, and alternatives.
  4. Rolling-Origin Evaluation for Seasonal Forecastsimplementation. Tests models at operational horizons.
  5. Intervention Analysis for Maintenance and Deployment Eventspractical application. Prevents known changes from being misread as ordinary stochastic dynamics.

References

Author not identified. (n.d.). Introduction to time series [Course notes, Modules 7-8]. Full citation details could not be confirmed.

Box, G. E. P., Jenkins, G. M., Reinsel, G. C., & Ljung, G. M. (2015). Time series analysis: Forecasting and control (5th ed.). Wiley.

Brockwell, P. J., & Davis, R. A. (2016). Introduction to time series and forecasting (3rd ed.). Springer. https://doi.org/10.1007/978-3-319-29854-2

Dickey, D. A., & Fuller, W. A. (1979). Distribution of the estimators for autoregressive time series with a unit root. Journal of the American Statistical Association, 74(366a), 427–431. https://doi.org/10.1080/01621459.1979.10482531

Hyndman, R. J., & Athanasopoulos, G. (2021). Forecasting: Principles and practice (3rd ed.). OTexts. https://otexts.com/fpp3/

R Core Team. (n.d.). R documentation. R Foundation for Statistical Computing. https://stat.ethz.ch/R-manual/R-devel/library/stats/html/00Index.html

Said, S. E., & Dickey, D. A. (1984). Testing for unit roots in autoregressive-moving average models of unknown order. Biometrika, 71(3), 599–607. https://doi.org/10.1093/biomet/71.3.599

Shumway, R. H., & Stoffer, D. S. (2025). Time series analysis and its applications: With R examples (5th ed.). Springer. https://doi.org/10.1007/978-3-031-70584-7

Research and Verification Notes

  • Uploaded material used: M7L1 for ARIMA, unit roots, ADF, fitting, and the SES connection; M7L2 for SARIMA and AirPassengers; M8L1 for multi-step AR forecasts, random walk with drift, innovation-based forecast variance, AirPassengers, and BJsales.
  • Primary-source verification: Dickey and Fuller (1979) and Said and Dickey (1984) support the unit-root testing discussion.
  • Authoritative verification: ARIMA/SARIMA notation and forecast formulas were checked against Brockwell and Davis (2016), Shumway and Stoffer (2025), Box et al. (2015), and Hyndman and Athanasopoulos (2021).
  • Official documentation: Current R arima and predict.Arima documentation was checked for implementation context.
  • Clarifications added: Deterministic versus stochastic trend, ADF limitations, over-differencing, interval widening, and retransformation bias.
  • Original material: The condition-feature workflow and portfolio reconstruction chain are original explanatory additions.
  • Code status: R code was not executed in this environment.
  • Review warning: Confirm transformation order, coefficient signs, drift conventions, and forecast-object behavior in the exact R package used for publication.