Time Series Forecasting

What is Time Series Forecasting?

A time series is a sequence of observations indexed by time: {y₁, y₂, …, yₙ} where the ordering is fundamental. Forecasting predicts future values ŷₙ₊ₕ given past history.

Properties that distinguish time series from i.i.d. data: - Temporal dependency: yₜ depends on yₜ₋₁, yₜ₋₂, … - Non-stationarity: mean, variance, and autocorrelation can change over time - Trend: long-run increase or decrease - Seasonality: periodic patterns (hourly, weekly, yearly) - Cycles: irregular multi-year oscillations - Noise: random fluctuations

Decomposition

Separate a time series into structural components:

Additive: yₜ = Trend + Seasonal + Residual (when amplitudes don't scale with level) Multiplicative: yₜ = Trend × Seasonal × Residual (when seasonality scales with trend)

from statsmodels.tsa.seasonal import seasonal_decompose
import pandas as pd

ts = pd.read_csv("data.csv", index_col="date", parse_dates=True)["value"]
result = seasonal_decompose(ts, model="additive", period=12)
result.plot()

Stationarity

Most classical methods require stationarity: statistical properties (mean, variance, autocorrelation) don't change over time.

Augmented Dickey-Fuller (ADF) test: - H₀: unit root (non-stationary) - p < 0.05: reject H₀ → evidence of stationarity

from statsmodels.tsa.stattools import adfuller

result = adfuller(ts)
print(f"ADF Statistic: {result[0]:.4f}")
print(f"p-value: {result[1]:.4f}")
print(f"Critical values: {result[4]}")

Making series stationary: - Differencing: ∇yₜ = yₜ − yₜ₋₁ (removes trend) - Seasonal differencing: ∇_s yₜ = yₜ − yₜ₋ₛ (removes seasonality) - Log transform: stabilizes variance for exponentially growing series - Box-Cox transform: generalization of log

Classical Statistical Models

Moving Average (MA)

Simple moving average smooths noise: MA(q): ŷₜ = (1/q) Σₖ₌₀^{q-1} yₜ₋ₖ

Exponential Moving Average (EMA): down-weights older observations: ŷₜ = α yₜ + (1−α) ŷₜ₋₁

Exponential Smoothing (Holt-Winters)

Handles level, trend, and seasonality with separate smoothing equations:

Level: lₜ = α(yₜ − sₜ₋ₘ) + (1−α)(lₜ₋₁ + bₜ₋₁) Trend: bₜ = β(lₜ − lₜ₋₁) + (1−β)bₜ₋₁ Season: sₜ = γ(yₜ − lₜ) + (1−γ)sₜ₋ₘ Forecast: ŷₜ₊ₕ = lₜ + h·bₜ + sₜ₊ₕ₋ₘ

from statsmodels.tsa.holtwinters import ExponentialSmoothing

model = ExponentialSmoothing(ts, trend="add", seasonal="add", seasonal_periods=12)
fit   = model.fit(optimized=True)
forecast = fit.forecast(12)   # 12 steps ahead

ARIMA (AutoRegressive Integrated Moving Average)

Combines three components: - AR(p): regression on p lagged values: yₜ = c + φ₁yₜ₋₁ + … + φₚyₜ₋ₚ + εₜ - I(d): d-th order differencing to achieve stationarity - MA(q): regression on q lagged errors: yₜ = c + εₜ + θ₁εₜ₋₁ + … + θqεₜ₋q

Full model: ARIMA(p, d, q)

Choosing p, d, q: - d: number of differences needed for stationarity (ADF test) - p: from PACF plot — look for cutoff - q: from ACF plot — look for cutoff - Or use auto_arima to search automatically

from statsmodels.tsa.arima.model import ARIMA
from pmdarima import auto_arima

# Manual
model = ARIMA(ts, order=(2, 1, 2))
fit   = model.fit()
print(fit.summary())
forecast = fit.forecast(steps=12)

# Automatic order selection
auto_model = auto_arima(ts, seasonal=True, m=12,
                          stepwise=True, trace=False)
print(auto_model.summary())

SARIMA (Seasonal ARIMA)

Extends ARIMA with seasonal terms: SARIMA(p,d,q)(P,D,Q)_m

where m is the seasonal period (12 for monthly, 7 for daily with weekly seasonality).

Vector Autoregression (VAR)

Multivariate ARIMA: model multiple correlated time series jointly. Each variable is a linear function of lagged values of all variables.

from statsmodels.tsa.vector_ar.var_model import VAR

model = VAR(df_multivariate)   # df with multiple columns
results = model.fit(maxlags=12, ic="aic")   # AIC selects lag order
forecast = results.forecast(df_multivariate.values[-12:], steps=6)

Facebook Prophet

Additive model designed for business time series with: - Piecewise linear or logistic trend with changepoints - Fourier series seasonality at multiple frequencies - Holiday effects

from prophet import Prophet

df = pd.DataFrame({"ds": dates, "y": values})
model = Prophet(
    yearly_seasonality=True,
    weekly_seasonality=True,
    daily_seasonality=False,
    changepoint_prior_scale=0.05,   # flexibility of trend
    seasonality_prior_scale=10.0
)
model.add_country_holidays(country_name="US")
model.fit(df)

future  = model.make_future_dataframe(periods=90)
forecast = model.predict(future)
model.plot(forecast)
model.plot_components(forecast)

Strengths: robust to missing data, automatic changepoint detection, handles multiple seasonalities, accessible API. Weaknesses: univariate only, not best for noisy high-frequency data.

Machine Learning for Time Series

Feature Engineering

Transform time series into tabular form by creating features from the past:

def make_lag_features(df, target_col, lags, windows):
    feat = df.copy()
    for lag in lags:
        feat[f"lag_{lag}"] = feat[target_col].shift(lag)
    for w in windows:
        feat[f"roll_mean_{w}"] = feat[target_col].shift(1).rolling(w).mean()
        feat[f"roll_std_{w}"]  = feat[target_col].shift(1).rolling(w).std()
        feat[f"roll_max_{w}"]  = feat[target_col].shift(1).rolling(w).max()
    # Date features
    feat["hour"]    = feat.index.hour
    feat["weekday"] = feat.index.weekday
    feat["month"]   = feat.index.month
    feat["quarter"] = feat.index.quarter
    return feat.dropna()

df_feat = make_lag_features(df, "sales", lags=[1,7,14,28], windows=[7,14,28])

Then fit any regression model (XGBoost, LightGBM are particularly effective).

Time Series Cross-Validation

Standard k-fold is wrong — future data would appear in training. Use walk-forward validation:

from sklearn.model_selection import TimeSeriesSplit

tscv = TimeSeriesSplit(n_splits=5, gap=7)   # gap=7 avoids leakage
for fold, (train_idx, val_idx) in enumerate(tscv.split(X)):
    X_tr, X_val = X[train_idx], X[val_idx]
    y_tr, y_val = y[train_idx], y[val_idx]
    model.fit(X_tr, y_tr)
    preds = model.predict(X_val)
    print(f"Fold {fold}: MAE={mean_absolute_error(y_val, preds):.4f}")

Deep Learning for Time Series

1D Convolutions

Treat time series like a 1D image; learn temporal patterns:

class TCN(nn.Module):
    def __init__(self, n_inputs, n_outputs, kernel_size=3):
        super().__init__()
        self.conv = nn.Sequential(
            nn.Conv1d(n_inputs, 64, kernel_size, padding=kernel_size//2),
            nn.ReLU(),
            nn.Conv1d(64, 64, kernel_size, padding=kernel_size//2),
            nn.ReLU(),
            nn.Conv1d(64, n_outputs, 1)   # output projection
        )
    def forward(self, x):
        # x shape: (batch, channels, seq_len)
        return self.conv(x)[:, :, -1]    # take last time step

Dilated causal convolutions (WaveNet) increase receptive field exponentially without increasing parameters.

LSTM for Forecasting

class LSTMForecaster(nn.Module):
    def __init__(self, n_features, hidden_size, horizon):
        super().__init__()
        self.lstm = nn.LSTM(n_features, hidden_size, batch_first=True, num_layers=2)
        self.head = nn.Linear(hidden_size, horizon)

    def forward(self, x):
        # x: (batch, seq_len, n_features)
        out, _ = self.lstm(x)
        return self.head(out[:, -1, :])  # (batch, horizon)

Temporal Fusion Transformer (TFT)

Combines LSTM with self-attention and explicit variable selection. Handles: - Multiple static covariates - Known future inputs (holidays, promotions) - Unknown future inputs (economic indicators) - Multi-quantile forecasting output

from pytorch_forecasting import TemporalFusionTransformer, TimeSeriesDataSet

# Define dataset
training = TimeSeriesDataSet(
    data, time_idx="time_idx", target="value",
    group_ids=["store_id", "product_id"],
    max_encoder_length=52, max_prediction_length=12,
    static_categoricals=["store_id"],
    time_varying_known_reals=["price", "promotion"],
    time_varying_unknown_reals=["value"]
)

N-BEATS, N-HiTS, PatchTST

Modern pure deep learning forecasters: - N-BEATS: deep stacks of residual blocks with interpretable trend/seasonality components - N-HiTS: hierarchical interpolation — different scales processed separately - PatchTST: Vision-Transformer approach — treat patches of the time series like image patches - iTransformer: invert attention — attend over variates, not time steps

Forecasting Metrics

MASE (Mean Absolute Scaled Error): normalizes by the naive one-step forecast (random walk). MASE < 1 means the model beats the naive baseline.

Handling Special Cases

Intermittent Demand (Many Zeros)

Use Croston's method or ADIDA; or Tweedie/negative binomial distribution with XGBoost.

Multiple Related Series (Hierarchical)

Top-down: forecast aggregate, disaggregate proportionally. Bottom-up: forecast each series, aggregate. Reconciliation (optimal): fit all levels jointly, reconcile for consistency (scikit-hts, statsforecast).

Anomaly Detection in Time Series

• STL decomposition → check residuals against 3σ threshold
• Isolation Forest on lag features
• LSTM reconstruction error (high error = anomaly)
• Streaming methods: CUSUM, ADWIN

Benchmark Comparison