ARIMA + Prophet

ARIMA Model Overview

ARIMA (AutoRegressive Integrated Moving Average) combines three components:

DfARIMA

AutoRegressive Integrated Moving Average — a class of models that describes a time series as a linear function of its own past values (AR), past forecast errors (MA), and differencing to achieve stationarity (I). The model is parameterized by orders (p, d, q) for the AR, differencing, and MA components respectively.

\text{ARIMA}(p, d, q)

p = AR order (AutoRegressive)
d = degree of differencing (Integrated)
q = MA order (Moving Average)

Architecture Diagram

ARIMA Components:
┌─────────────────────────────────────────────────────────────┐
│                                                             │
│  AR(p): y_t = c + phi_1*y_{t-1} + ... + phi_p*y_{t-p} + epsilon_t │
│         -> Uses past values to predict current               │
│                                                             │
│  I(d):  Δ^d y_t  -> Differencing to achieve stationarity   │
│                                                             │
│  MA(q): y_t = c + epsilon_t + theta_1*epsilon_{t-1} + ... + theta_q*epsilon_{t-q} │
│         -> Uses past errors to predict current               │
│                                                             │
└─────────────────────────────────────────────────────────────┘

Mathematical Formulation

\phi(B)(1-B)^d y_t = \theta(B)\epsilon_t

Where:

$B$ is the backshift operator: $By_t = y_{t-1}$
$\phi(B) = 1 - \phi_1 B - \phi_2 B^2 - \ldots - \phi_p B^p$
$\theta(B) = 1 + \theta_1 B + \theta_2 B^2 + \ldots + \theta_q B^q$
$\epsilon_t \sim \text{WN}(0, \sigma^2)$ (white noise)

The backshift operator $B$ simplifies the ARIMA notation. For example, first-order differencing $(1-B)y_t = y_t - y_{t-1}$ removes a linear trend. The AR polynomial $\phi(B)$ captures autocorrelation, and the MA polynomial $\theta(B)$ captures the structure of forecast errors.

ACF and PACF for Model Selection

import numpy as np
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX
import pandas as pd
import warnings
warnings.filterwarnings('ignore')

# Generate AR(2) process
np.random.seed(42)
n = 500
ar_params = [0.6, -0.3]
y_ar = np.zeros(n)
for t in range(2, n):
    y_ar[t] = ar_params[0]*y_ar[t-1] + ar_params[1]*y_ar[t-2] + np.random.randn()

fig, axes = plt.subplots(1, 3, figsize=(15, 4))
axes[0].plot(y_ar, linewidth=0.8)
axes[0].set_title('AR(2) Process')

plot_acf(y_ar, lags=20, ax=axes[1], title='ACF - Tails off')
plot_pacf(y_ar, lags=20, ax=axes[2], title='PACF - Cuts off at lag 2')

plt.tight_layout()
plt.savefig('acf_pacf_ar.png', dpi=150)
plt.show()

Selection Rules

Architecture Diagram

ACF/PACF Patterns:
┌────────────┬─────────────────┬─────────────────┐
│  Process   │      ACF        │      PACF       │
├────────────┼─────────────────┼─────────────────┤
│  AR(p)     │  Tails off      │  Cuts off lag p │
│  MA(q)     │  Cuts off lag q │  Tails off      │
│  ARMA(p,q) │  Tails off      │  Tails off      │
└────────────┴─────────────────┴─────────────────┘

Fitting ARIMA Models

# Create real-world-like time series
dates = pd.date_range('2020-01-01', periods=365*3, freq='D')
t = np.arange(len(dates))

y_series = (
    0.01 * t +
    10 * np.sin(2 * np.pi * t / 365) +
    5 * np.sin(2 * np.pi * t / 30) +
    np.random.randn(len(t)) * 3
)
ts = pd.Series(y_series, index=dates)

# Split train/test
train_size = int(len(ts) * 0.8)
train, test = ts[:train_size], ts[train_size:]

print(f"Train: {len(train)} observations")
print(f"Test:  {len(test)} observations")

Auto-ARIMA (Grid Search)

DfAIC (Akaike Information Criterion)

A measure of model quality that balances goodness of fit against model complexity. Lower AIC indicates a better model. It penalizes models with more parameters to prevent overfitting: $\text{AIC} = 2k - 2\ln(\hat{L})$ where k is the number of parameters and $\hat{L}$ is the maximized likelihood.

# Manual ARIMA selection with AIC
import itertools

p_range = range(0, 4)
d_range = range(0, 3)
q_range = range(0, 4)

best_aic = np.inf
best_order = None
results_list = []

for p, d, q in itertools.product(p_range, d_range, q_range):
    try:
        model = ARIMA(train, order=(p, d, q))
        fitted = model.fit()
        results_list.append({
            'order': (p, d, q),
            'aic': fitted.aic,
            'bic': fitted.bic
        })
        if fitted.aic < best_aic:
            best_aic = fitted.aic
            best_order = (p, d, q)
    except:
        continue

results_df = pd.DataFrame(results_list).sort_values('aic')
print("Top 5 Models by AIC:")
print(results_df.head())
print(f"\nBest order: {best_order}")

Fit and Forecast with ARIMA

📝ARIMA Forecasting Workflow

# Fit best model
model = ARIMA(train, order=best_order)
fitted_model = model.fit()

print(fitted_model.summary())

# Forecast
forecast_result = fitted_model.get_forecast(steps=len(test))
forecast_mean = forecast_result.predicted_mean
forecast_ci = forecast_result.conf_int(alpha=0.05)

# Plot
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

axes[0].plot(train.index, train, label='Train', color='blue')
axes[0].plot(test.index, test, label='Test', color='green')
axes[0].plot(test.index, forecast_mean, label='Forecast', color='red', linestyle='--')
axes[0].fill_between(
    test.index,
    forecast_ci.iloc[:, 0],
    forecast_ci.iloc[:, 1],
    alpha=0.2, color='red'
)
axes[0].set_title('ARIMA Forecast')
axes[0].legend()

axes[1].plot(test.index[:60], test.values[:60], label='Actual', marker='o', markersize=3)
axes[1].plot(test.index[:60], forecast_mean.values[:60], label='Forecast', marker='x', markersize=3)
axes[1].fill_between(
    test.index[:60],
    forecast_ci.iloc[:60, 0],
    forecast_ci.iloc[:60, 1],
    alpha=0.2
)
axes[1].set_title('Forecast (First 60 Days)')
axes[1].legend()

plt.tight_layout()
plt.savefig('arima_forecast.png', dpi=150)
plt.show()

SARIMA: Seasonal ARIMA

For data with seasonal patterns:

\text{SARIMA}(p, d, q)(P, D, Q, s)

DfSARIMA

Seasonal ARIMA extends ARIMA by adding seasonal AR, differencing, and MA components. The notation (p, d, q)(P, D, Q, s) specifies non-seasonal orders (p, d, q) and seasonal orders (P, D, Q) with seasonal period s. Seasonal differencing $(1-B^s)^D$ removes seasonal patterns.

\Phi_P(B^s)\phi(B)(1-B)^d(1-B^s)^D y_t = \Theta_Q(B^s)\theta(B)\epsilon_t

Where $s$ is the seasonal period.

# Seasonal ARIMA for monthly data (s=12)
sarima_model = SARIMAX(
    train,
    order=(1, 1, 1),
    seasonal_order=(1, 1, 1, 12),
    enforce_stationarity=False,
    enforce_invertibility=False
)

sarima_fitted = sarima_model.fit(disp=False)
print(sarima_fitted.summary())

# Forecast
sarima_forecast = sarima_fitted.get_forecast(steps=len(test))
sarima_mean = sarima_forecast.predicted_mean
sarima_ci = sarima_forecast.conf_int()

# Plot
plt.figure(figsize=(12, 5))
plt.plot(train.index[-120:], train.values[-120:], label='Train')
plt.plot(test.index, test.values, label='Test', color='green')
plt.plot(test.index, sarima_mean, label='SARIMA Forecast', color='red', linestyle='--')
plt.fill_between(test.index, sarima_ci.iloc[:, 0], sarima_ci.iloc[:, 1], alpha=0.2, color='red')
plt.title('SARIMA(1,1,1)(1,1,1,12) Forecast')
plt.legend()
plt.tight_layout()
plt.savefig('sarima_forecast.png', dpi=150)
plt.show()

When selecting SARIMA parameters, first determine the non-seasonal (p, d, q) using ACF/PACF on the seasonally differenced series. Then determine the seasonal (P, D, Q) by examining ACF/PACF at lags that are multiples of the seasonal period s. The auto_arima function automates this process.

Facebook Prophet

Prophet is designed for business time series with strong seasonal patterns, missing data, historical trend changes, and holiday effects.

Mathematical Model

y(t) = g(t) + s(t) + h(t) + \epsilon_t

Where:

$g(t)$ = trend function (logistic or linear)
$s(t)$ = periodic seasonal component
$h(t)$ = holiday/event effects
$\epsilon_t$ = error term

DfProphet

A forecasting library by Facebook (Meta) designed for business time series with strong seasonal effects, trend changes, and holiday effects. Prophet fits an additive regression model with trend, seasonal, and holiday components. The trend uses either a linear function or a logistic growth model with automatic changepoint detection.

from prophet import Prophet

# Prepare data for Prophet (requires 'ds' and 'y' columns)
df_prophet = pd.DataFrame({
    'ds': ts.index,
    'y': ts.values
})

# Split
train_p = df_prophet[:train_size]
test_p = df_prophet[train_size:]

# Initialize and fit Prophet
prophet_model = Prophet(
    growth='linear',
    yearly_seasonality=True,
    weekly_seasonality=False,
    daily_seasonality=False,
    changepoint_prior_scale=0.05,
    seasonality_prior_scale=10.0
)

prophet_model.fit(train_p)

# Make future dataframe
future = prophet_model.make_future_dataframe(periods=len(test_p))
forecast = prophet_model.predict(future)

# Plot
fig1 = prophet_model.plot(forecast)
plt.title('Prophet Forecast')
plt.savefig('prophet_forecast.png', dpi=150)
plt.show()

# Components
fig2 = prophet_model.plot_components(forecast)
plt.savefig('prophet_components.png', dpi=150)
plt.show()

Adding Holidays

# Define custom holidays
holidays = pd.DataFrame({
    'holiday': 'special_event',
    'ds': pd.to_datetime(['2021-06-15', '2022-06-15', '2023-06-15']),
    'lower_window': -1,
    'upper_window': 1,
})

# Fit with holidays
prophet_holiday = Prophet(
    holidays=holidays,
    holidays_prior_scale=10.0
)
prophet_holiday.fit(train_p)

forecast_holiday = prophet_holiday.predict(future)

Adding Custom Seasonality

# Add monthly seasonality
prophet_custom = Prophet(
    yearly_seasonality=False,
    weekly_seasonality=False,
    daily_seasonality=False
)

# Add custom Fourier seasonality
prophet_custom.add_seasonality(
    name='monthly',
    period=30.5,  # ~30.5 days
    fourier_order=5
)

prophet_custom.fit(train_p)
forecast_custom = prophet_custom.predict(future)

Prophet models seasonality using Fourier series. The fourier_order parameter controls the flexibility of the seasonal pattern: higher values capture more complex patterns but risk overfitting. The default yearly seasonality uses Fourier order 10, which can capture up to 5 harmonics per year.

Forecasting Evaluation Metrics

📝Comprehensive Forecast Evaluation

from sklearn.metrics import (
    mean_absolute_error,
    mean_squared_error,
    mean_absolute_percentage_error
)

def evaluate_forecast(actual, predicted, model_name=''):
    """Calculate comprehensive forecast metrics."""
    mae = mean_absolute_error(actual, predicted)
    rmse = np.sqrt(mean_squared_error(actual, predicted))
    mape = mean_absolute_percentage_error(actual, predicted) * 100

    # Mean Absolute Scaled Error (MASE)
    naive_forecast = actual.shift(1).dropna()
    actual_trimmed = actual[1:]
    mae_naive = mean_absolute_error(actual_trimmed, naive_forecast)
    mase = mae / mae_naive if mae_naive > 0 else np.inf

    print(f"\n{model_name} Evaluation:")
    print(f"  MAE:  {mae:.4f}")
    print(f"  RMSE: {rmse:.4f}")
    print(f"  MAPE: {mape:.2f}%")
    print(f"  MASE: {mase:.4f}")

    return {'mae': mae, 'rmse': rmse, 'mape': mape, 'mase': mase}

# Evaluate ARIMA
arima_metrics = evaluate_forecast(test, forecast_mean, 'ARIMA')

# Evaluate SARIMA
sarima_metrics = evaluate_forecast(test, sarima_mean, 'SARIMA')

# Evaluate Prophet
prophet_pred = forecast['yhat'].iloc[-len(test):]
prophet_metrics = evaluate_forecast(test.values, prophet_pred.values, 'Prophet')

# Comparison
print("\n" + "="*50)
print("Model Comparison:")
print(f"{'Model':<12} {'MAE':>8} {'RMSE':>8} {'MAPE':>8}")
print("-"*40)
for name, metrics in [('ARIMA', arima_metrics), ('SARIMA', sarima_metrics), ('Prophet', prophet_metrics)]:
    print(f"{name:<12} {metrics['mae']:>8.4f} {metrics['rmse']:>8.4f} {metrics['mape']:>7.2f}%")

DfMASE (Mean Absolute Scaled Error)

A scale-free forecast accuracy metric that compares the model's MAE to the MAE of a naive (random walk) forecast. MASE < 1 means the model outperforms the naive forecast. Unlike MAPE, MASE is well-defined when actual values are zero or near zero.

Residual Diagnostics

from statsmodels.stats.diagnostic import acorr_ljungbox

def check_residuals(model, lags=20):
    """Check if residuals are white noise."""
    residuals = model.resid

    fig, axes = plt.subplots(2, 2, figsize=(12, 8))

    axes[0, 0].plot(residuals, linewidth=0.5)
    axes[0, 0].set_title('Residuals')
    axes[0, 0].axhline(y=0, color='r', linestyle='--')

    axes[0, 1].hist(residuals, bins=30, edgecolor='black', alpha=0.7)
    axes[0, 1].set_title('Residual Distribution')

    plot_acf(residuals, lags=lags, ax=axes[1, 0], title='ACF of Residuals')

    from scipy import stats
    stats.probplot(residuals, dist="norm", plot=axes[1, 1])
    axes[1, 1].set_title('Q-Q Plot')

    plt.tight_layout()
    plt.savefig('residual_diagnostics.png', dpi=150)
    plt.show()

    # Ljung-Box test
    lb_test = acorr_ljungbox(residuals, lags=[10, 20], return_df=True)
    print("\nLjung-Box Test (H0: residuals are white noise):")
    print(lb_test)
    print("\nIf p-values > 0.05, residuals are likely white noise ✓")

check_residuals(fitted_model)

Residual diagnostics are critical for validating an ARIMA model. If residuals show significant autocorrelation (Ljung-Box test p < 0.05), the model has not captured all the temporal structure. Q-Q plots check normality of residuals — deviations from the diagonal suggest non-normal error distribution, which affects prediction intervals but not point forecasts.

Model Selection Guide

Architecture Diagram

┌─────────────────────────────────────────────────────────────────┐
│                    Which model to choose?                       │
├─────────────────────────────────────────────────────────────────┤
│                                                                 │
│  Data has strong seasonality?                                   │
│  ├─ Yes → SARIMA or Prophet                                     │
│  └─ No  → ARIMA                                                 │
│                                                                 │
│  Need to add holidays/events?                                   │
│  ├─ Yes → Prophet (easiest)                                     │
│  └─ No  → ARIMA/SARIMA                                          │
│                                                                 │
│  Data has changepoints?                                         │
│  ├─ Yes → Prophet (best) or piecewise ARIMA                     │
│  └─ No  → ARIMA                                                 │
│                                                                 │
│  Quick baseline needed?                                         │
│  ├─ Yes → auto_arima                                            │
│  └─ No  → Prophet (tunable)                                     │
│                                                                 │
│  Interpretability important?                                    │
│  ├─ Yes → ARIMA (coefficients interpretable)                    │
│  └─ No  → Prophet or ensemble                                   │
│                                                                 │
└─────────────────────────────────────────────────────────────────┘

Key Takeaways

ARIMA is the foundation: AR captures autocorrelation, I makes series stationary, MA captures shocks
ACF/PACF guide order selection: AR cuts off in PACF, MA cuts off in ACF
SARIMA extends ARIMA for seasonal data with $(P,D,Q,s)$ parameters
Prophet excels at business forecasting with holidays and trend changes
Always check residuals: They should be white noise (no autocorrelation)
Use AIC/BIC for model comparison: lower is better
Multiple metrics (MAE, RMSE, MAPE, MASE) give a complete picture of forecast quality

📋Summary: ARIMA + Prophet

ARIMA(p,d,q) combines autoregression (AR), differencing (I), and moving average (MA) via the backshift operator: $\phi(B)(1-B)^d y_t = \theta(B)\epsilon_t$ .
ACF/PACF patterns guide order selection: AR(p) → PACF cuts off at lag p; MA(q) → ACF cuts off at lag q; ARMA → both tail off → use AIC/BIC.
SARIMA(p,d,q)(P,D,Q,s) adds seasonal components with period s. Seasonal differencing $(1-B^s)^D$ removes seasonal patterns.
Prophet fits $y(t) = g(t) + s(t) + h(t) + \epsilon_t$ with automatic changepoint detection and Fourier-based seasonality. Ideal for business data with holidays.
AIC $= 2k - 2\ln(\hat{L})$ balances fit and complexity for model selection — lower is better.
Residual diagnostics: Ljung-Box test (white noise check), Q-Q plot (normality), ACF of residuals (no autocorrelation).
MASE compares forecast accuracy to a naive benchmark: MASE < 1 means the model beats the random walk.
Forecast evaluation requires multiple metrics: MAE (scale-dependent), RMSE (penalizes large errors), MAPE (percentage), MASE (scale-free).
Always use log transforms or Box-Cox before ARIMA for multiplicative seasonality or variance stabilization.
Choose ARIMA for interpretable linear models, SARIMA for seasonal data, and Prophet for business forecasting with trend changes and holidays.

Practice Exercises

ARIMA Tuning: Fit ARIMA(p,d,q) for all combinations of p,d,q ∈ {0,1,2}. Which has lowest AIC?
Seasonal Data: Apply SARIMA to monthly temperature data. Try different seasonal periods (12, 365).
Prophet vs ARIMA: Compare both on a dataset with trend changes. Which adapts better?
Residual Analysis: After fitting, test residuals with Ljung-Box. What do significant p-values indicate?

ARIMA + Prophet

ARIMA + Prophet

ARIMA Model Overview

DfARIMA

Mathematical Formulation

ACF and PACF for Model Selection

Selection Rules

Fitting ARIMA Models

Auto-ARIMA (Grid Search)

DfAIC (Akaike Information Criterion)

Fit and Forecast with ARIMA

📝ARIMA Forecasting Workflow

SARIMA: Seasonal ARIMA

DfSARIMA

Facebook Prophet

Mathematical Model

DfProphet

Adding Holidays

Adding Custom Seasonality

Forecasting Evaluation Metrics

📝Comprehensive Forecast Evaluation

DfMASE (Mean Absolute Scaled Error)

Residual Diagnostics

Model Selection Guide

Key Takeaways

📋Summary: ARIMA + Prophet

Practice Exercises

Need Expert Data Science Help?