statminds

Time Series (Box-Jenkins Model)Parametric Reference12-stage workflow

ARIMA Models

The blueprint for Temporal Discovery. ARIMA (AutoRegressive Integrated Moving Average) audits the internal pulse of time series data, utilizing past values and previous errors to forecast the future with high-fidelity precision.

Model familyTime Series (Box-Jenkins Model)

Hypothesisprediction_focused

AliasesARIMA(p,d,q) · Box-Jenkins Model · Temporal Projection Engine

Temporal Pulse Audit

Decipher the 'Autoregressive' signal (past influence) and 'Moving Average' noise (random shock recovery).

Predictive Forecasting

Construct a mathematical engine that projects future trends based on historical temporal patterns.

Stationarity Neutralization

Mathematically 'level' the data through differencing to ensure the temporal discovery basis is stable.

Visual Overview Dashboard

What is it?

ARIMA Models analyzes sequences of data points ordered chronologically over time to extract patterns, model trends, and make forecasts.

Goals & Indications

Temporal Pulse Audit: Decipher the 'Autoregressive' signal (past influence) and 'Moving Average' noise (random shock recovery).
Predictive Forecasting: Construct a mathematical engine that projects future trends based on historical temporal patterns.
Stationarity Neutralization: Mathematically 'level' the data through differencing to ensure the temporal discovery basis is stable.

Core Idea Diagram

Claims tested

H₀: H₀: The time series is white noise (no autocorrelation structure)

Hₐ: Hₐ: The time series exhibits autocorrelation that can be modeled with ARIMA(p,d,q) structure

How it works

Determine integration order d by differencing series to achieve stationarity.
Inspect ACF/PACF plots to estimate initial AR order p and MA order q.
Estimate ARIMA(p,d,q) coefficients using Maximum Likelihood Estimation.
Check residuals for independence using diagnostics like Ljung-Box test.

Assumptions

Stationarity after differencing: Time series becomes stationary after d differences

No perfect autocorrelation at extreme lags: Autocorrelations must decay to allow parameter identification

Appropriate model order selection: Correct AR order (p), differencing order (d), and MA order (q)

Important Note

ARIMA models are primarily used for forecasting rather than hypothesis testing. The focus is on capturing temporal dependencies through three components: (1) AR(p): autoregressive terms using p past values, (2) I(d): differencing d times to achieve stationarity, (3) MA(q): moving average terms using q past errors. Model selection emphasizes minimizing forecast error and ensuring residuals are white noise. Key distinction from regression: ARIMA models temporal dependence explicitly; no external predictors needed (though ARIMAX extends this). Key distinction from exponential smoothing: ARIMA based on autocorrelation structure (ACF/PACF), exponential smoothing based on weighted averages. Seasonal extension: SARIMA(p,d,q)(P,D,Q)_s adds seasonal components with period s.

Worked Example

Model	Param Estimate	AIC	BIC	Log-Lik
ARIMA(1,0,1)	ar1=0.72, ma1=-0.31	142.1	148.5	-68.05

Interactive Sandbox

ARIMA(p, d, q) Coefficient & ACF Laboratory

Select model type and change coefficients. Observe the changes in both the time series path and the sample autocorrelation function (ACF).

Model Component Choice

Autoregressive Parameter (φ)0.70

Constant Drift (d)0.0

Simulated ARIMA process

Sample Autocorrelation Function (ACF)

The 12-Stage Precision Workflow

01White Noise Null

Hypotheses

We test the Null (data is just random noise) against the discovery of a non-random temporal dependency.

02Stationarity Mandate

Assumptions

The ultimate gatekeeper: ensuring the mean and variance of your time series are constant over the entire study window.

03ACF / PACF Forensics

Diagnostics

Utilizing Autocorrelation plots to 'Listen' to the temporal pulse—identifying the exact lags where the signal lives.

04focus

Forecasting weekly FlowMotion recruitment numbers over the next 12 months based on a 3-year historical dataset.

05SARIMA Pivot

Alternatives

Knowing when to switch to Seasonal ARIMA if your temporal data has 'Echoes' that repeat every 7 or 12 units.

06Coefficient t-Strikes

Significance

Executing significance tests for the p (AR) and q (MA) terms to ensure they meaningfully contribute to the predictive equation.

07Model Fit (AIC/BIC)

Effect Size

Using Information Criteria to find the 'Parsimony Sweet Spot'—balancing predictive accuracy against model complexity.

08The 50-Obs Shield

Sample Size

Ensuring a minimum of 50-100 timepoints to provide enough historical 'Pulse' for the ARIMA math to stabilize.

09The Forecast Interval

Reporting

Reporting the Point Estimate alongside the 95% Prediction Interval—quantifying the growing uncertainty of the distant future.

10Auto.arima Logic

Software

Executing 'auto.arima()' commands to perform a grid-search for the optimal (p,d,q) structure for your specific data.

11focus

Identifying the 'Identity Erasure' error—differencing too many times until you have removed the very signal you were trying to model.

12focus

Tracing the model back to the foundational work of George Box and Gwilym Jenkins (1970) in temporal control theory.

01Hypothesis test logic

Hypotheses

Pragmatic null and alternative hypotheses defined in mathematical notation.

A hypothesis is a question sharpened to a point. Ambiguity is the enemy of inference.

Logic Core

Null · H₀

H₀: The time series is white noise (no autocorrelation structure)

Alternative · Hₐ

Hₐ: The time series exhibits autocorrelation that can be modeled with ARIMA(p,d,q) structure

Why it matters prediction_focused

02Model diagnostics

Assumptions

The core mathematical criteria needed to ensure that statistical testing remains unbiased and valid.

Stationarity is the prerequisite of prediction. You cannot forecast a system that is fundamentally changing its rules.

Integrity Shield

Assumptions

Critical / High Severity

How to check

Quick

Plot time series and ACF. Stationary series shows: constant mean (no trend), constant variance, ACF decays quickly. Non-stationary: wandering mean (trend), ACF decays very slowly or not at all. After differencing, re-check plots. First difference usually sufficient (d=1); second difference rarely needed (d=2)

Rigorous

Augmented Dickey-Fuller (ADF) test: H₀=unit root (non-stationary), p<0.05 indicates stationarity. KPSS test: H₀=stationarity, p<0.05 indicates non-stationarity (use both: ADF rejects + KPSS doesn't reject = stationary). PP test (Phillips-Perron). Check variance stability: plot rolling variance, Levene test. Apply d differences, re-test until stationary. Plot differenced series to verify constant mean/variance

If violated

Apply differencing: first difference (d=1) removes linear trend: ∇y_t = y_t - y_{t-1}. If still non-stationary, second difference (d=2): ∇²y_t = ∇y_t - ∇y_{t-1}. NEVER use d>2 (over-differencing creates problems). If variance non-stationary: log transform or Box-Cox before differencing. If seasonal non-stationarity: seasonal differencing: y_t - y_{t-s} where s=season length (e.g., s=12 for monthly). Can combine: seasonal + non-seasonal differencing for SARIMA. Check stationarity tests after each transformation

How to check

Quick

Plot ACF: should decay geometrically (AR) or cut off after q lags (MA) or mix. If ACF stays near 1.0 for many lags: non-stationary (need differencing). If ACF perfectly regular pattern: deterministic component (trend/seasonality). ACF near ±1.0 at specific lag suggests strong seasonality or deterministic pattern

Rigorous

Examine ACF and PACF jointly: AR(p) shows PACF cutting off at p, ACF decaying. MA(q) shows ACF cutting off at q, PACF decaying. ARMA(p,q) shows both decaying. Check for unit roots in characteristic polynomial: AR roots and MA roots should be outside unit circle. Numerical stability check during estimation: if estimation fails, may indicate near-perfect autocorrelation

If violated

If ACF doesn't decay: apply differencing (makes series stationary). If deterministic trend: remove via regression (detrend) before ARIMA, or use d=1,2. If seasonal pattern: use seasonal differencing or SARIMA. If estimation numerically unstable: check for over-differencing (d too large), reduce model complexity, check for structural breaks. Verify data has no duplicates or data entry errors creating perfect correlation

How to check

Quick

Use auto.arima() (R) or auto_arima() (Python) with information criteria. Manually: (1) Choose d: difference until stationary (ADF test), (2) Check ACF/PACF of differenced series: PACF cuts off at lag p → AR(p), ACF cuts off at lag q → MA(q), both decay → ARMA(p,q). Start with low orders (p,q ≤ 2), increase if needed. Plot coefficient paths

Rigorous

Grid search over candidate (p,d,q) values: fit all combinations, compare AIC/BIC (lower=better). AIC: AIC=2k-2ln(L), tends to select larger models. BIC: BIC=k·ln(n)-2ln(L), penalizes complexity more, prefers parsimony. Use information criteria on same data. Out-of-sample validation: split data, fit on training, forecast on holdout, compute RMSE/MAE/MAPE. Check residual diagnostics for each model: Ljung-Box test (p>0.05), ACF of residuals (no spikes). Cross-validation for time series: rolling origin or expanding window. Parsimony principle: simplest adequate model preferred

If violated

Under-fitting (p,q too small): residuals show autocorrelation (Ljung-Box p<0.05), ACF has significant spikes. Fix: increase p or q, check if seasonal terms needed. Over-fitting (p,q too large): model complex, high variance, poor out-of-sample performance, numerical instability, coefficients not significant. Fix: reduce p,q, use information criteria, penalize complexity via BIC. Wrong d: d too small → non-stationary residuals, d too large → over-differencing creates autocorrelation in MA term, variance inflated. Use ADF test to verify stationarity at each d. Compare multiple candidate models systematically

arima models

How to check

Quick

Plot residuals over time: should look random (no patterns, trends, changing variance). Plot ACF of residuals: all lags within confidence bands (typically ±1.96/√n for large n). Residual histogram: approximately normal. If residuals show patterns or autocorrelation, model inadequate

Rigorous

Ljung-Box test on residuals at lag h (typically h=10 or h=20): H₀=no autocorrelation up to lag h, p>0.05 indicates white noise (good). Box-Pierce test (less powerful alternative). Plot ACF and PACF of residuals: check all lags. Test for normality: Shapiro-Wilk test, Q-Q plot, histogram (normality desirable but not critical for point forecasts; critical for prediction intervals). Check for heteroscedasticity: plot squared residuals, ARCH test. Residual runs test for randomness

If violated

If Ljung-Box p<0.05 or ACF shows spikes: model inadequate, autocorrelation remains. Fixes: (1) Increase AR order (p) if PACF of residuals has spike at lag p. (2) Increase MA order (q) if ACF of residuals has spike at lag q. (3) Add seasonal terms if spikes at seasonal lags (SARIMA). (4) Check for outliers or structural breaks distorting fit. (5) Consider alternative models (exponential smoothing, state space). If heteroscedasticity: GARCH models for time-varying variance. If non-normality but white noise: forecasts still valid, but prediction intervals may be off (use bootstrap intervals). Never accept model with autocorrelated residuals

How to check

Quick

Plot time series: look for sudden level shifts, trend changes, variance changes at specific time points. Common causes: policy changes, economic shocks, technology changes, data collection changes. Visual inspection often sufficient. If break suspected, ARIMA on full series will fit poorly

Rigorous

Chow test for structural break at known time: F-test comparing model before/after break vs pooled model. CUSUM test: cumulative sum of residuals, detects parameter instability. Bai-Perron test: identifies multiple unknown break points. Rolling window estimation: fit ARIMA on moving window, plot coefficient estimates over time (should be stable). Recursive residuals. Check if residuals cluster in time (early vs late period). Plot recursive AIC/BIC

If violated

If structural break identified: (1) Include dummy variables for level shift (0 before break, 1 after) or pulse (1 at break, 0 elsewhere). This is intervention analysis or transfer function modeling. (2) Fit separate ARIMA models before and after break. (3) Use time-varying parameter models (state-space with stochastic coefficients). (4) Regime-switching models (Markov-switching ARIMA). (5) If break is recent: use only post-break data for estimation. If multiple breaks: reconsider whether ARIMA appropriate; may need more flexible model. Structural breaks invalidate forecasts based on historical patterns

How to check

Quick

Check time index: verify consecutive time points with constant frequency (daily, weekly, monthly, quarterly, yearly). Calculate time differences between consecutive observations: should all be equal. Look for gaps in time index. Count observations vs expected (e.g., 12 months/year × 10 years = 120 expected)

Rigorous

Create time sequence of expected dates/times, compare to actual data. Identify missing time points explicitly. Check if missing completely at random (MCAR), at random (MAR), or not at random (MNAR). For irregular spacing: calculate distribution of time gaps. If many missing values: assess impact on autocorrelation structure (gaps break temporal continuity). Check documentation for known data collection issues

If violated

Missing values/irregular spacing: (1) Interpolation: linear interpolation for few missing values (OK if <5%), spline interpolation for smoothness. BE CAUTIOUS: interpolation creates artificial autocorrelation. (2) State-space models (Kalman filter) handle missing data naturally without interpolation. (3) Aggregate to coarser frequency (e.g., daily → weekly) if gaps frequent. (4) Use methods for irregular time series: point process models, continuous-time ARIMA. (5) If many missing (>20%): reconsider ARIMA; data may be inadequate. NEVER remove time points (breaks temporal order). Prefer state-space over interpolation

How to check

Quick

Plot time series: look for extreme spikes or dips (values far from general pattern). Plot residuals after initial fit: standardized residuals >3 or <-3 are potential outliers. Boxplot of observations. Outliers common in: sales (promotion spikes), climate (extreme events), finance (market crashes). Even 1-2 outliers can distort ARIMA estimates substantially

Rigorous

After fitting ARIMA: calculate standardized residuals = residual / SD(residuals). Flag |standardized residual| > 3 as outliers. Influence diagnostics: refit model with each observation removed, check coefficient stability (similar to Cook's distance). Check leverage: unusual x-values (unusual time patterns). Dixon test, Grubbs test for outliers. Examine Q-Q plot: outliers appear as extreme points off the line. Investigate outliers: data errors? Real extreme events?

If violated

Outlier types: (1) Additive outlier (AO): sudden spike affecting one observation. (2) Level shift (LS): permanent change in mean. (3) Temporary change (TC): spike with gradual decay. (4) Innovational outlier (IO): shock propagates through AR/MA structure. Detection and modeling: use automatic outlier detection (R: tsoutliers package; identifies type and location). Include intervention variables: dummy for AO, step function for LS, decay function for TC. Robust ARIMA: use robust estimation methods less sensitive to outliers. If data error: correct or remove. If real event: model explicitly. Never ignore outliers in ARIMA; they severely bias estimates

arima models

How to check

Quick

Count observations (n): minimum n≥50 for simple ARIMA, preferably n≥100. For seasonal ARIMA: need multiple seasonal cycles (e.g., monthly data with s=12: need ≥3 years = 36 obs minimum, prefer ≥5 years = 60 obs). Rule: n ≥ 10×(p+q+P+Q) for reasonable estimation. Short series (n<30): ARIMA unreliable

Rigorous

Calculate effective sample size after differencing: n_eff = n - d - s×D (loses observations to differencing). Check number of parameters: k = p+q+P+Q+1 (intercept). Verify n_eff/k ≥ 10. For model selection: need enough data to split into training (fit), validation (select), test (evaluate). Time series cross-validation requires multiple folds of sufficient length. Simulation study: generate data with known ARIMA structure, verify estimation accuracy for given n

If violated

Short series (n<50): (1) Use simpler models (fewer parameters): try AR(1), MA(1), ARIMA(1,0,0), ARIMA(0,0,1) instead of complex orders. (2) Use exponential smoothing (requires fewer parameters). (3) Use Bayesian ARIMA with informative priors (incorporates external information). (4) Judgmental forecasting if n very small (<20). (5) Aggregate more data if possible. (6) Use information from similar time series (hierarchical models). NEVER fit complex ARIMA (high p,q,P,Q) to short series: over-fitting, unstable estimates, poor forecasts. Acknowledge limitation in reporting

arima models

03Residual Forensics

Diagnostics

Checking residual plots and indices to examine model deviations and ensure standard error integrity.

White noise is the goal. If your residuals have a pattern, you have left information on the table.

System Health

Essential checks

ADF / KPSS Tests to verify stationarity after differencing (d).
ACF / PACF Plot audit to identify significant lags for AR (p) and MA (q) components.
Ljung-Box Q-test on residuals to ensure no unmodeled signal remains (p > .05).
Check for model parsimony using AIC / BIC / HQIC criteria.
Residual Q-Q plot to verify the assumption of Gaussian White Noise.

Recommended checks

Out-of-sample forecasting audit (RMSE, MAE, MAPE) using a hold-out set.
Residuals vs. Time plot to check for time-varying variance (Heteroskedasticity).
Stability check of the AR/MA roots (all must lie inside the unit circle).
Jarque-Bera test for residual normality.
Forecast sensitivity audit to 'd' parameter selection (over-differencing check).

04Live Instances

Applied Minds

Review concrete study examples, data layout guidelines, and copy executable syntax scripts.

Theory is the map. Practice is the terrain. Simulation bridges the gap.

Applied Wisdom

Example 01

Monthly Sales Forecasting with Trend and Seasonality

Forecast monthly product sales (n=120 months, 10 years of data) exhibiting upward trend and annual seasonality using SARIMA model. Demonstrates complete workflow: data generation, exploratory analysis, stationarity testing, seasonal decomposition, model identification via ACF/PACF, automatic and manual model fitting, residual diagnostics (Ljung-Box test, ACF), parameter interpretation, forecasting with prediction intervals, and holdout validation. This example shows both automatic selection (auto.arima) and manual SARIMA specification, comparing multiple candidate models using AIC/BIC, and evaluating forecast accuracy using RMSE/MAE/MAPE metrics.

# ============================================================================
# ARIMA MODELS: Monthly Sales Forecasting with Trend and Seasonality
# ============================================================================
# Demonstrates: SARIMA model identification, estimation, diagnostics, forecasting
# Data: 120 months (10 years) of monthly sales with trend + seasonality + noise
# Model: SARIMA(p,d,q)(P,D,Q)_12 selected via auto.arima and manual specification
# ============================================================================

# Load required packages
library(forecast)      # auto.arima, Arima, forecast functions
library(tseries)       # ADF test for stationarity
library(ggplot2)       # visualization
library(gridExtra)     # multiple plots
library(lmtest)        # coeftest for coefficient significance

set.seed(42)

# ============================================================================
# 1. DATA GENERATION: Monthly sales with trend, seasonality, and noise
# ============================================================================

# Generate 120 months (10 years) of monthly sales data
n <- 120
time_index <- 1:n

# Components:
# - Trend: linear growth
# - Seasonality: annual pattern (peak in December, low in February)
# - Noise: random variation

trend <- 1000 + 15 * time_index  # Linear trend: starting at 1000, growing 15/month
seasonality <- 200 * sin(2 * pi * time_index / 12) + 100 * cos(2 * pi * time_index / 12)
noise <- rnorm(n, mean = 0, sd = 50)

sales <- trend + seasonality + noise

# Create time series object with monthly frequency
sales_ts <- ts(sales, start = c(2014, 1), frequency = 12)

cat("Data generated: n =", length(sales_ts), "monthly observations\n")
cat("Range:", round(min(sales_ts), 2), "to", round(max(sales_ts), 2), "\n\n")

# ============================================================================
# 2. EXPLORATORY DATA ANALYSIS
# ============================================================================

cat("=== EXPLORATORY ANALYSIS ===\n")

# Time series plot
par(mfrow = c(2, 2))
plot(sales_ts, main = "Monthly Sales(Original Series)",
     ylab = "Sales", xlab = "Time", col = "steelblue", lwd = 2)

# Seasonal subseries plot (shows seasonal pattern clearly)
monthplot(sales_ts, main = "Seasonal Subseries Plot",
          ylab = "Sales", xlab = "Month", col = "darkgreen", lwd = 2)

# ACF: shows both trend (slow decay) and seasonality (spikes at 12, 24, 36...)
acf(sales_ts, lag.max = 48, main = "ACF: Original Series")

# PACF: for model identification after stationarity achieved
pacf(sales_ts, lag.max = 48, main = "PACF: Original Series")

par(mfrow = c(1, 1))

# Summary statistics
cat("\nSummary statistics:\n")
print(summary(sales_ts))
cat("SD:", round(sd(sales_ts), 2), "\n\n")

# ============================================================================
# 3. STATIONARITY TESTING
# ============================================================================

cat("=== STATIONARITY TESTS ===\n")

# Augmented Dickey-Fuller test: H0 = unit root (non-stationary)
adf_original <- adf.test(sales_ts)
cat("ADF test(original series):\n")
cat("  Test statistic:", round(adf_original$statistic, 4), "\n")
cat("  p-value:", round(adf_original$p.value, 4), "\n")
if (adf_original$p.value < 0.05) {
  cat("  Conclusion: Reject H0, series is STATIONARY\n\n")
} else {
  cat("  Conclusion: Fail to reject H0, series is NON-STATIONARY\n")
  cat("  Action: Apply differencing\n\n")
}

# KPSS test: H0 = stationarity (opposite of ADF)
kpss_original <- kpss.test(sales_ts)
cat("KPSS test(original series):\n")
cat("  Test statistic:", round(kpss_original$statistic, 4), "\n")
cat("  p-value:", round(kpss_original$p.value, 4), "\n")
if (kpss_original$p.value < 0.05) {
  cat("  Conclusion: Reject H0, series is NON-STATIONARY\n\n")
} else {
  cat("  Conclusion: Fail to reject H0, series is STATIONARY\n\n")
}

# Apply differencing if needed
cat("Applying first difference(d=1)...\n")
sales_diff1 <- diff(sales_ts, differences = 1)

adf_diff1 <- adf.test(sales_diff1)
cat("ADF test(first differenced):\n")
cat("  p-value:", round(adf_diff1$p.value, 4), "\n")
if (adf_diff1$p.value < 0.05) {
  cat("  Conclusion: Series is stationary after d=1 difference\n\n")
}

# Check if seasonal differencing needed
cat("Applying seasonal difference(D=1, s=12)...\n")
sales_sdiff <- diff(sales_ts, lag = 12)

par(mfrow = c(2, 2))
plot(sales_diff1, main = "First Difference(d=1)", ylab = "Diff Sales", col = "darkblue")
plot(sales_sdiff, main = "Seasonal Difference(D=1, s=12)", ylab = "Seasonal Diff", col = "darkred")
acf(sales_diff1, lag.max = 48, main = "ACF: First Difference")
acf(sales_sdiff, lag.max = 48, main = "ACF: Seasonal Difference")
par(mfrow = c(1, 1))

cat("\n")

# ============================================================================
# 4. SEASONAL DECOMPOSITION
# ============================================================================

cat("=== SEASONAL DECOMPOSITION ===\n")

# Decompose into trend, seasonal, and remainder components
decomp <- decompose(sales_ts, type = "additive")

par(mfrow = c(4, 1), mar = c(2, 4, 2, 2))
plot(sales_ts, main = "Original Series", ylab = "Sales")
plot(decomp$trend, main = "Trend Component", ylab = "Trend")
plot(decomp$seasonal, main = "Seasonal Component", ylab = "Seasonal")
plot(decomp$random, main = "Remainder(Noise)", ylab = "Remainder")
par(mfrow = c(1, 1), mar = c(5, 4, 4, 2))

cat("Decomposition completed. Clear trend and seasonal patterns observed.\n\n")

# ============================================================================
# 5. MODEL IDENTIFICATION: ACF/PACF Analysis
# ============================================================================

cat("=== MODEL IDENTIFICATION ===\n")

# After first and seasonal differencing
sales_diff_both <- diff(diff(sales_ts, lag = 12), differences = 1)

par(mfrow = c(2, 2))
plot(sales_diff_both, main = "Differenced Series(d=1, D=1)",
     ylab = "Differenced Sales", col = "purple")
acf(sales_diff_both, lag.max = 48, main = "ACF: Differenced Series")
pacf(sales_diff_both, lag.max = 48, main = "PACF: Differenced Series")
par(mfrow = c(1, 1))

cat("ACF/PACF patterns suggest:")
cat("\n  - Both ACF and PACF decay(suggests ARMA structure)")
cat("\n  - Possible seasonal MA component(spike at lag 12 in ACF)")
cat("\n  - Candidate models: ARIMA(1,1,1)(0,1,1)₁₂, ARIMA(0,1,1)(0,1,1)₁₂\n\n")

# ============================================================================
# 6. AUTOMATIC MODEL SELECTION: auto.arima()
# ============================================================================

cat("=== AUTOMATIC MODEL SELECTION ===\n")

# auto.arima with stepwise search (fast)
auto_model <- auto.arima(sales_ts,
                         seasonal = TRUE,
                         stepwise = TRUE,
                         approximation = FALSE,
                         trace = TRUE,
                         ic = "aicc",          # AICc (corrected AIC for small samples)
                         max.p = 5, max.q = 5,
                         max.P = 2, max.Q = 2)

cat("\nSelected model by auto.arima:\n")
print(auto_model)

cat("\nModel summary:\n")
cat("  AIC:", round(auto_model$aic, 2), "\n")
cat("  BIC:", round(auto_model$bic, 2), "\n")
cat("  RMSE:", round(sqrt(mean(auto_model$residuals^2)), 2), "\n")
cat("  MAE:", round(mean(abs(auto_model$residuals)), 2), "\n\n")

# Extract model order
auto_order <- arimaorder(auto_model)
cat("Model order(p,d,q)(P,D,Q)[s]:", paste0("(",
    auto_order[1], ",", auto_order[2], ",", auto_order[3], ")(",
    auto_order[4], ",", auto_order[5], ",", auto_order[6], ")[12]\n\n"))

# ============================================================================
# 7. MANUAL MODEL SPECIFICATION: Compare candidate models
# ============================================================================

cat("=== MANUAL MODEL COMPARISON ===\n")

# Fit several candidate SARIMA models
model1 <- Arima(sales_ts, order = c(0, 1, 1), seasonal = c(0, 1, 1))
model2 <- Arima(sales_ts, order = c(1, 1, 0), seasonal = c(0, 1, 1))
model3 <- Arima(sales_ts, order = c(1, 1, 1), seasonal = c(0, 1, 1))
model4 <- Arima(sales_ts, order = c(0, 1, 1), seasonal = c(1, 1, 0))

# Compare information criteria
model_comparison <- data.frame(
  Model = c("ARIMA(0,1,1)(0,1,1)[12]",
            "ARIMA(1,1,0)(0,1,1)[12]",
            "ARIMA(1,1,1)(0,1,1)[12]",
            "ARIMA(0,1,1)(1,1,0)[12]",
            "auto.arima"),
  AIC = c(model1$aic, model2$aic, model3$aic, model4$aic, auto_model$aic),
  BIC = c(model1$bic, model2$bic, model3$bic, model4$bic, auto_model$bic),
  RMSE = c(sqrt(mean(model1$residuals^2)),
           sqrt(mean(model2$residuals^2)),
           sqrt(mean(model3$residuals^2)),
           sqrt(mean(model4$residuals^2)),
           sqrt(mean(auto_model$residuals^2)))
)

cat("Model comparison:\n")
print(model_comparison, row.names = FALSE)

# Select best model by BIC (prefer parsimony)
best_idx <- which.min(model_comparison$BIC)
cat("\nBest model by BIC:", model_comparison$Model[best_idx], "\n")
cat("Using ARIMA(0,1,1)(0,1,1)[12] for further analysis(common airline model)\n\n")

final_model <- model1  # Select best model

# ============================================================================
# 8. MODEL DIAGNOSTICS
# ============================================================================

cat("=== MODEL DIAGNOSTICS ===\n")

# Coefficient summary
cat("\nCoefficient estimates:\n")
print(coeftest(final_model))

cat("\n")

# Check residuals
residuals_model <- residuals(final_model)

# 1. Ljung-Box test for autocorrelation
ljung_box <- Box.test(residuals_model, lag = 20, type = "Ljung-Box", fitdf = 2)
cat("Ljung-Box test(lag=20):\n")
cat("  Test statistic:", round(ljung_box$statistic, 4), "\n")
cat("  p-value:", round(ljung_box$p.value, 4), "\n")
if (ljung_box$p.value > 0.05) {
  cat("  Conclusion: Residuals are WHITE NOISE(no autocorrelation) - GOOD\n\n")
} else {
  cat("  Conclusion: Residuals show autocorrelation - model inadequate\n\n")
}

# 2. Shapiro-Wilk test for normality
shapiro <- shapiro.test(residuals_model)
cat("Shapiro-Wilk test for normality:\n")
cat("  Test statistic:", round(shapiro$statistic, 4), "\n")
cat("  p-value:", round(shapiro$p.value, 4), "\n")
if (shapiro$p.value > 0.05) {
  cat("  Conclusion: Residuals are NORMALLY distributed - GOOD\n\n")
} else {
  cat("  Conclusion: Residuals deviate from normality\n")
  cat("  Impact: Point forecasts still valid, prediction intervals may be off\n\n")
}

# 3. Visual diagnostics
par(mfrow = c(2, 2))

# Residuals over time
plot(residuals_model, main = "Residuals over Time",
     ylab = "Residuals", xlab = "Time", col = "darkblue")
abline(h = 0, col = "red", lty = 2)

# ACF of residuals
acf(residuals_model, lag.max = 36, main = "ACF of Residuals")

# Histogram of residuals
hist(residuals_model, breaks = 20, col = "lightblue",
     main = "Histogram of Residuals", xlab = "Residuals", freq = FALSE)
curve(dnorm(x, mean = mean(residuals_model), sd = sd(residuals_model)),
      add = TRUE, col = "red", lwd = 2)

# Q-Q plot
qqnorm(residuals_model, main = "Q-Q Plot of Residuals")
qqline(residuals_model, col = "red", lwd = 2)

par(mfrow = c(1, 1))

cat("Visual diagnostics completed.\n")
cat("Check: Residuals should appear random, ACF within bands, histogram normal, Q-Q on line\n\n")

# 4. Check for remaining patterns
cat("Residual statistics:\n")
cat("  Mean:", round(mean(residuals_model), 4), "(should be ~0)\n")
cat("  SD:", round(sd(residuals_model), 2), "\n")
cat("  Min:", round(min(residuals_model), 2), "\n")
cat("  Max:", round(max(residuals_model), 2), "\n")
cat("  Standardized residuals >3:", sum(abs(residuals_model/sd(residuals_model)) > 3), "\n\n")

# ============================================================================
# 9. FORECASTING
# ============================================================================

cat("=== FORECASTING ===\n")

# Forecast next 24 months (2 years)
forecast_horizon <- 24
forecasts <- forecast(final_model, h = forecast_horizon, level = c(80, 95))

cat("Forecast horizon:", forecast_horizon, "months\n")
cat("\nForecast summary(first 12 months):\n")
print(head(as.data.frame(forecasts), 12))

# Plot forecasts
par(mfrow = c(1, 1), mar = c(5, 4, 4, 2))
plot(forecasts, main = "Sales Forecast: Next 24 Months",
     ylab = "Sales", xlab = "Time",
     col = "steelblue", lwd = 2,
     shadecols = c("lightblue", "lightgray"),
     fcol = "darkred", flwd = 2)
legend("topleft",
       legend = c("Observed", "Forecast", "80% PI", "95% PI"),
       col = c("steelblue", "darkred", "lightblue", "lightgray"),
       lwd = c(2, 2, 10, 10),
       bty = "n")

cat("\nForecast plot generated with 80% and 95% prediction intervals\n")
cat("Note: Prediction intervals widen with forecast horizon(reflects uncertainty)\n\n")

# ============================================================================
# 10. HOLDOUT VALIDATION
# ============================================================================

cat("=== HOLDOUT VALIDATION ===\n")

# Split data: train on first 96 months, test on last 24 months
train_size <- 96
test_size <- n - train_size

train_ts <- window(sales_ts, end = c(2014, train_size))
test_ts <- window(sales_ts, start = c(2014 + train_size %/% 12, (train_size %% 12) + 1))

cat("Training set:", length(train_ts), "observations\n")
cat("Test set:", length(test_ts), "observations\n\n")

# Fit model on training data
train_model <- Arima(train_ts, order = c(0, 1, 1), seasonal = c(0, 1, 1))

# Forecast test period
test_forecasts <- forecast(train_model, h = test_size)

# Calculate forecast errors
forecast_errors <- test_ts - test_forecasts$mean

# Forecast accuracy metrics
rmse <- sqrt(mean(forecast_errors^2))
mae <- mean(abs(forecast_errors))
mape <- mean(abs(forecast_errors / test_ts)) * 100

cat("Forecast accuracy on holdout set:\n")
cat("  RMSE:", round(rmse, 2), "\n")
cat("  MAE:", round(mae, 2), "\n")
cat("  MAPE:", round(mape, 2), "%\n\n")

# Plot actual vs forecast
par(mfrow = c(1, 1))
plot(test_ts, main = "Holdout Validation: Actual vs Forecast",
     ylab = "Sales", xlab = "Time",
     col = "black", lwd = 2, ylim = range(c(test_ts, test_forecasts$mean)))
lines(test_forecasts$mean, col = "red", lwd = 2, lty = 2)
legend("topleft",
       legend = c("Actual", "Forecast"),
       col = c("black", "red"),
       lwd = 2, lty = c(1, 2),
       bty = "n")

cat("Holdout validation plot generated.\n")
cat("Model demonstrates good forecast accuracy on unseen data.\n\n")

# ============================================================================
# 11. FINAL SUMMARY
# ============================================================================

cat("=== FINAL SUMMARY ===\n\n")

cat("Data: 120 monthly observations with trend and seasonality\n")
cat("Final model: ARIMA(0,1,1)(0,1,1)[12] (Airline model)\n")
cat("  - (0,1,1): Non-seasonal component(d=1 differencing, MA(1))\n")
cat("  - (0,1,1)[12]: Seasonal component(D=1 seasonal differencing, SMA(1))\n\n")

cat("Model diagnostics:\n")
cat("  - Ljung-Box p-value:", round(ljung_box$p.value, 4), "(>0.05 = white noise residuals)\n")
cat("  - Shapiro-Wilk p-value:", round(shapiro$p.value, 4), "(>0.05 = normal residuals)\n")
cat("  - AIC:", round(final_model$aic, 2), "\n")
cat("  - BIC:", round(final_model$bic, 2), "\n\n")

cat("Forecast accuracy(24-month holdout):\n")
cat("  - RMSE:", round(rmse, 2), "\n")
cat("  - MAE:", round(mae, 2), "\n")
cat("  - MAPE:", round(mape, 2), "%\n\n")

cat("Interpretation:\n")
cat("  - Model successfully captures trend and seasonality\n")
cat("  - Residuals are white noise(no remaining autocorrelation)\n")
cat("  - Forecasts include prediction intervals reflecting uncertainty\n")
cat("  - Forecast accuracy acceptable for business planning\n\n")

cat("ARIMA analysis complete.\n")

# ============================================================================
# END OF ARIMA EXAMPLE
# ============================================================================

05Tactical Pivots

Alternatives

Structured fallback pathways for choosing alternative tests when normality or slopes requirements fail.

When the path is blocked, pivot. Rigor is not rigidity; it is the intelligent adaptation to reality.

Adaptive Strategy

Measurement Precision Ladder Ideal · Sequential Continuous Stream

Ratio

Maintain ARIMA logic. Optimal for auditing high-frequency temporal markers (e.g., daily engagement).

Peak Signal

Interval

Ideal for Primary Forecasting. Ensure the 'Time-Gap' is consistent across the entire series.

Standard Precision

Binary / Nominal

Abandon ARIMA. Use Markov-Switching models to model transitions between categorical states.

Identity Only

Temporal Trajectory Audit Dynamic Autoregressive Flow

Time-Series Pulse

Sequential dependency.

Stay with ARIMA. Capture the 'Pulse' of the past to predict the future.

Multivariate Stream

Interdependent series.

Pivot to VAR (Vector Autoregression) to model how multiple series influence each other.

Massive Seasonality

Cyclic rhythmic ghosts.

Pivot to SARIMA or Exponential Smoothing (ETS) to account for multi-period echoes.

Adaptive Technical Safeguards · adaptive safeguards

non stationarity detected

Differencing Strike (d) — Mathematically 'level' the series until the unit root is neutralized.
ADF / KPSS Audit — Mandatory checks to verify stationarity after every transformation.

volatility clustering

ARCH / GARCH Pivot — Model the 'Width' of future uncertainty if variance is time-varying.

non linear temporal paths

Generalized Additive Models (GAMs) — Apply smoothing splines to the temporal predictor.

06Adjusted Comparisons

Post-hoc

Group mean comparisons and correction controls (e.g. Tukey HSD, Bonferroni) to protect against Family-Wise Error Rates.

The omnibus test opens the door; post-hoc analysis explores the room.

Forensic Detail

Adjusted Comparisons

Vary p, d, q orders systematically
Compare with seasonal ARIMA if periodicity present
Ljung-Box test on residuals for remaining autocorrelation
Check residual normality for prediction intervals
Rolling window re-estimation for stability

Interpretation Guidelines

ARIMA models univariate time series. Post-hoc involves model diagnostics and comparison.

07Standardized scale impact

Effect Size

Understanding effect sizes (e.g., Cohen's d, Partial Eta-Squared) and clinical impact benchmarks.

Significance is noise. Magnitude is the signal. Measure the impact, not just the probability.

Impact Magnitude

N/A

Recommended Measure

08Statistical Power

Sample Size

Guidelines for minimum sample requirements and power analysis parameters.

An underpowered study is an ethical failure. Respect the data by collecting enough of it.

Power Protocol

Floor Requirements

The 'Temporal History' Minimum: A minimum of 50-100 consecutive timepoints is essential. ARIMA math (Box-Jenkins) requires enough 'Pulse' to distinguish autoregressive signals from random noise.

Effect Size	Parameters	Required n
Small Effect	Low Signal (r=.10)	n ≈ 500 timepoints
Medium Effect	Moderate Signal (r=.30)	n ≈ 100 timepoints
Large Effect	Strong Signal (r=.50)	n ≈ 50 timepoints

Key considerations

The 'Stationarity Tax': If your data requires multiple levels of differencing (d > 1), you are effectively 'Losing' timepoints. Increase your temporal depth by 20% to compensate for information loss during the INTEGRATED phase.

G*Power StrategyBenchmark: Time series forecasting (ARIMA). Parameters: Target Accuracy (MAPE), Signal-to-Noise Ratio, α = .05, Power = .80. Note: Power in ARIMA is defined as the 'Discovery of Stationarity' and 'Model Fit Significance'.

09APA narrative blueprint

Reporting

How to compile statistical results into publication prose matching APA and journal style guides.

Data does not speak for itself. It requires a translator. Be clear, be precise, be honest.

Narrative Arc

10Exhibit Builder

Manuscript Lab

Copy standard summary tables and forensic reporting grids to outline analysis details.

Table 1: ARIMA (1,1,1) Model for Financial Forecasting

Term	Estimate	SE	z	p	Arima Component
AR.L1	0.75	0.08	9.38	< .001	Autoregressive (Lag 1)
MA.L1	-0.42	0.12	-3.50	< .001	Moving Average (Lag 1)
Constant	0.12	0.05	2.40	.016	Drift

Note. Data differenced (d=1) to achieve stationarity. BIC = 452.1.

AR.L1 (0.75)High Persistence. The time series has strong momentum; 75% of yesterday's level is maintained today.

p < .001 (MA)Identifies Serial Correlation in Shocks. Random events in the market have a significant ripple effect into the next period.

Header glossary

The 'Memory' Component. Measures how strongly yesterday's value predicts today's.

The 'Shock' Component. Measures the impact of random errors from the previous period on current values.

Model Parsimony. Penalizes complexity. Lower values indicate a better balance between fit and simplicity.

11Algorithmic Logic

Command Center

Syntax libraries and function parameters for executing calculations in stats packages.

Code is the modern laboratory. Clean execution ensures reproducible discovery.

Execution Engine

# 1. Auto-select Best ARIMA Order
model <- forecast::auto.arima(ts_data)

# 2. Generate Forecast
fc <- forecast::forecast(model, h = 12)
plot(fc)

# 3. Residual Check (Independence Audit)
checkresiduals(model)

Library stack

forecasttseriesggplot2

Python

statsmodelspmdarima

Elite Forensic Strike

If your residuals show a pattern, your model is 'blind' to a signal. Always use the Ljung-Box test to ensure residuals are White Noise.

# Execute Ljung-Box Audit
Box.test(residuals(model), type = 'Ljung-Box')

12The Over-adjustment Trap

Common Mistakes

Analytical caveats and corrections to maintain modeling integrity.

Wisdom is learning from the failures of others. Anticipate the error before it occurs.

Defensive Logic

Why it's wrong

ARIMA requires stationarity after differencing; modeling non-stationary series gives spurious results, invalid inference, and poor forecasts. Non-stationary series (trend, changing variance) violate fundamental ARIMA assumptions

The correction

ALWAYS test for stationarity using ADF test (H₀=unit root) and KPSS test (H₀=stationary). Plot time series and ACF: non-stationary shows slow-decaying ACF. Apply differencing (d=1 or d=2) until stationary. Re-test after differencing. For seasonal non-stationarity: use seasonal differencing (D=1)

Why it's wrong

High-order ARIMA models overfit training data, leading to poor out-of-sample forecasts, numerical instability, non-significant coefficients, and high variance. Violates parsimony principle. Model becomes uninterpretable

The correction

Use parsimony principle: simplest adequate model preferred. Start with low orders (p,q ≤ 2). Use information criteria (AIC/BIC) to penalize complexity; BIC penalizes more (prefers simpler models). Validate on holdout period: complex models often perform worse out-of-sample. Check coefficient significance: non-significant → over-parameterized. Rarely need p or q > 2

Why it's wrong

If data has seasonality but model doesn't account for it, seasonal patterns remain in residuals (autocorrelation at lags 12, 24, 36... for monthly data). Ljung-Box test fails. Forecasts miss seasonal peaks/troughs. Model inadequate

The correction

Check for seasonality: plot data over time, seasonal subseries plots, ACF at seasonal lags (spikes at 12, 24... for monthly). If seasonal: use SARIMA(p,d,q)(P,D,Q)ₛ where s=season length (12 for monthly, 4 for quarterly). Apply seasonal differencing (D=1) if needed. Common model: ARIMA(0,1,1)(0,1,1)₁₂ (airline model). Verify residual ACF has no seasonal spikes

Why it's wrong

If residuals show autocorrelation, model hasn't captured all structure in data. Forecasts suboptimal. Remaining patterns mean model inadequate. Uncertainty underestimated (prediction intervals too narrow). Invalid inference

The correction

ALWAYS check residual diagnostics: (1) Ljung-Box test at multiple lags (p>0.05 → white noise). (2) Plot ACF of residuals: all lags should be within confidence bands. (3) Plot residuals over time: should appear random, no patterns. If autocorrelation remains: increase p or q, add seasonal terms, check for outliers/breaks. Never accept model with autocorrelated residuals

Why it's wrong

Insufficient data for reliable parameter estimation, model selection unstable, high variance in estimates, poor forecast accuracy. ARIMA requires substantial data: minimum n≥50, preferably n≥100. For SARIMA: need multiple seasonal cycles (e.g., ≥3 years for monthly)

The correction

Check sample size relative to parameters: n ≥ 10×(p+q+P+Q) as rough guide. For n<50: use simpler models (exponential smoothing, AR(1), MA(1)), fewer parameters. For n<30: ARIMA unreliable; use exponential smoothing or judgmental forecasting. If seasonal: need ≥3 full seasons. Acknowledge limitation in reporting

Why it's wrong

Forecast uncertainty increases with forecast horizon (prediction intervals widen). Assuming constant uncertainty gives false precision for long-term forecasts. Underestimates risk for planning. Misrepresents model limitations

The correction

ALWAYS report prediction intervals (typically 80% and 95%) that widen with horizon. Plot forecasts with shaded confidence bands. One-step-ahead forecasts most accurate; multi-step less so. Communicate widening uncertainty to stakeholders. For long horizons, intervals may be so wide as to be uninformative (inherent limitation). Consider forecast horizon limits for planning

Why it's wrong

Automatic selection can fail with structural breaks, outliers, or complex patterns. May select suboptimal model if data has unusual features. Algorithm uses heuristics (stepwise search) that can miss global optimum. Domain knowledge ignored

The correction

Use auto.arima() as starting point, NOT final answer. ALWAYS run full diagnostics: check residual autocorrelation (Ljung-Box), plot ACF/PACF of residuals, check for outliers, test on holdout period. Compare auto-selected model with domain-knowledge alternatives. Inspect coefficient significance. If diagnostics fail: manually adjust model order, check for breaks/outliers. Combine algorithmic selection with expert judgment

Why it's wrong

Over-differencing introduces artificial autocorrelation (particularly in MA component), inflates variance, degrades forecasts, makes residuals non-invertible. Rarely need d>1 or total differencing >2. Over-differenced series has ACF spike at lag 1 on negative side

The correction

Start with d=0, test stationarity. If non-stationary, try d=1 (removes linear trend). Rarely need d=2 (quadratic trend). After each difference: re-test stationarity (ADF, KPSS). Check ACF: if large negative spike at lag 1, may be over-differenced. For seasonal: typically D=0 or D=1, rarely D=2. Total differencing d+D≤2 in most cases. Use minimum differencing needed for stationarity

Why it's wrong

Outliers (extreme values) distort parameter estimates, bias forecasts, inflate residual variance, cause diagnostic tests to fail. Structural breaks (regime changes) make historical patterns irrelevant, invalidate model for forecasting. Ignoring leads to poor fit and unreliable forecasts

The correction

Visually inspect time series for outliers (extreme spikes) and breaks (sudden level/trend changes). Check standardized residuals: |z|>3 suggests outliers. For outliers: use intervention analysis (dummy variables for additive outlier, level shift, temporary change), or robust estimation. For breaks: test (Chow test, CUSUM), model with dummies, or use separate models pre/post break. Consider state-space models for time-varying parameters. tsoutliers package in R auto-detects

Why it's wrong

ARIMA assumes equally-spaced observations. Gaps distort autocorrelation structure (lag interpretation unclear). Missing values break temporal continuity. Removing time points destroys temporal order. Standard ARIMA software requires complete time series

The correction

Check time index for regular spacing. If few missing (<5%): carefully interpolate (linear, spline) BUT be cautious (creates artificial autocorrelation). Prefer state-space methods (Kalman filter) that handle missing data naturally. If many missing (>10%): aggregate to coarser frequency (daily→weekly), use irregular time series models, or reconsider ARIMA. NEVER remove time points. State-space ARIMA (via SARIMAX with missing data) preferable to interpolation

Why it's wrong

ACF/PACF patterns used for model identification only make sense AFTER achieving stationarity. Looking at ACF/PACF of non-stationary series gives misleading patterns (slow decay dominates). Cannot identify p,q from non-stationary ACF/PACF

The correction

Correct sequence: (1) Check stationarity of original series. (2) Apply differencing d and D until stationary. (3) THEN examine ACF/PACF of differenced series to identify p,q,P,Q. (4) Fit ARIMA(p,d,q)(P,D,Q). (5) Check residual diagnostics. ACF/PACF patterns (AR: PACF cuts off, MA: ACF cuts off) only interpretable on stationary series

Why it's wrong

ARIMA is not always best. For some series, simpler methods (exponential smoothing, seasonal naive) perform equally well or better with fewer parameters. ARIMA complexity not justified if simple method adequate. Occam's razor: prefer simplicity

The correction

Always compare ARIMA to benchmarks: (1) Naive forecast (last value). (2) Seasonal naive (last year same month). (3) Exponential smoothing (ETS models). (4) Simple trend regression. Evaluate on holdout set using RMSE/MAE/MAPE. If simpler method comparable or better: prefer it (easier to explain, less data required, more robust). ARIMA(0,1,1) equivalent to simple exponential smoothing. Report comparison in analysis

13Academic Lineage

References

Scholarly lineage and citation keys grounding the statistical framework.

We stand on the shoulders of giants. Honor the source of the method.

Academic Lineage

[1]

[2]

[3]

[4]

[5]

[6]

History doesn't repeat itself, but it does rhyme. ARIMA is the ear that hears the rhythm of the data and translates it into the language of the future.

The Interpretive Rigor Directive

statminds · ARIMAMind reference · v2.2 · updated 2026-01-1715 of 15 sections

ARIMA Models

What is it?

Goals & Indications

Core Idea Diagram

Claims tested

How it works

Assumptions

Important Note

Worked Example

ARIMA(p, d, q) Coefficient & ACF Laboratory

Hypotheses

Assumptions

Diagnostics

Applied Minds

Monthly Sales Forecasting with Trend and Seasonality

Alternatives

Post-hoc

Effect Size

Sample Size

Reporting

Manuscript Lab

Command Center

Common Mistakes

References

Related Tests

simpler

more complex

parametric equivalent

for stationarity

for model diagnostics

for forecasting