Time Series Forecasting

Decomposing a Time Series: Trend, Seasonality, and Residuals

A time series $\{X_t\}_{t=1}^T$ is a sequence of observations recorded at equally spaced points in time. The first step in understanding such data is to separate the signal into interpretable components. The classical additive decomposition posits

$$X_t = T_t + S_t + R_t,$$

where $T_t$ denotes the trend (the long-term, slowly varying direction), $S_t$ is the seasonal component (regular, calendar-driven fluctuations), and $R_t$ is the residual (or irregular) component that captures everything else. When the amplitude of seasonal swings grows with the level of the series, a multiplicative decomposition $X_t = T_t \cdot S_t \cdot R_t$ is more appropriate; the multiplicative form can be transformed to an additive one by taking logarithms.

Trend is typically estimated by smoothing the series, for example with a moving average of order $m$:

$$\hat{T}_t = \frac{1}{m} \sum_{j=-(m-1)/2}^{(m-1)/2} X_{t+j} \quad\text{(centered moving average, $m$ odd)}.$$

The choice of $m$ determines how much high-frequency variability is removed. For economic data such as quarterly GDP, a four-quarter moving average eliminates seasonal fluctuations and yields a smooth estimate of the underlying growth path.

Seasonality captures repeating patterns within a fixed period — daily cycles, weekly patterns, or annual seasonality. After the trend is removed, seasonal factors can be estimated by averaging the de-trended values for each season (e.g., all Januaries in monthly data). The seasonal component is then normalized so that it sums to zero over a full cycle (additive) or averages to one (multiplicative).

The residual $R_t$ is defined as what remains after extracting trend and seasonality. Ideally, $R_t$ behaves like a stationary white-noise process — uncorrelated random shocks with constant variance. If the decomposition is adequate, the residuals should show no systematic pattern or temporal dependence. The presence of significant autocorrelation in $R_t$ signals that the model is misspecified.

From an economic standpoint, decomposition mirrors the way analysts think about macroeconomic indicators: policymakers care about the trend (potential output), the seasonal (holiday spending spikes), and the residual (unexpected shocks such as oil price disruptions). Isolating these components allows for more interpretable diagnostics and more defensible forecasts.

Autocorrelation: Measuring Linear Dependence Across Time

The autocorrelation function (ACF) quantifies how a series correlates with its own past values. For a weakly stationary process (constant mean, variance, and autocovariance depending only on lag), the autocovariance at lag $k$ is

$$\gamma_k = \mathrm{Cov}(X_t, X_{t-k}) = \mathbb{E}[(X_t-\mu)(X_{t-k}-\mu)],$$

and the theoretical autocorrelation at lag $k$ is

$$\rho_k = \frac{\gamma_k}{\gamma_0},\qquad \rho_0 = 1.$$

$\rho_k$ lies in $[-1,1]$ and measures the linear dependence between two observations $k$ time steps apart. A value close to $1$ ($-1$) indicates that if $X_{t-k}$ is above the mean, $X_t$ tends to be above (below) the mean by a proportional amount.

In practice we estimate $\rho_k$ from the sample. Let $\bar{x} = \frac{1}{T}\sum_{t=1}^T x_t$. The sample autocorrelation is

$$\hat{\rho}_k = \frac{\sum_{t=k+1}^T (x_t - \bar{x})(x_{t-k} - \bar{x})}{\sum_{t=1}^T (x_t - \bar{x})^2}, \quad k = 1,2,\dots$$

This is a biased but consistent estimator; the denominator uses all $T$ terms to guarantee that the sequence $\{\hat{\rho}_k\}$ forms a positive-semi-definite autocorrelation matrix. Under the null hypothesis that the data are independent white noise, Bartlett's approximation gives $\text{Var}[\hat{\rho}(k)] \approx 1/T$ for $k \ge 1$, yielding the familiar $95\%$ confidence bands at $\pm 1.96/\sqrt{T}$ that appear on every correlogram.

Interpreting the ACF yields insight into the data-generating process:

• A slow, linear decay of the ACF often signals a non-stationary trend (a random walk has $\rho_k\approx 1$ for all $k$).
• A damped sinusoidal decay suggests the presence of a stationary autoregressive component with complex roots, often associated with business-cycle or seasonal fluctuations.
• A sharp cutoff after lag $q$, where $\hat{\rho}_k$ is essentially zero for $k>q$, points to a moving-average model of order $q$, $\mathrm{MA}(q)$.
• A seasonal time series will exhibit large positive spikes at lags that are multiples of the seasonal period (12, 24, …), superimposed on the underlying pattern.

In financial economics, a significantly positive $\rho(1)$ in daily returns would suggest short-term predictability — a finding with profound implications for market efficiency.

Partial Autocorrelation: Isolating Direct Effects

While the ACF at lag $k$ mixes direct and indirect influences (via intermediate lags), the partial autocorrelation function (PACF) isolates the direct linear relationship between $X_t$ and $X_{t-k}$ after controlling for $X_{t-1}, X_{t-2}, \dots, X_{t-k+1}$. Formally, the partial autocorrelation at lag $k$, denoted $\phi_{kk}$, is the last coefficient of a linear regression of $X_t$ on its $k$ most recent values:

$$X_t = \phi_{k1} X_{t-1} + \phi_{k2} X_{t-2} + \dots + \phi_{kk} X_{t-k} + \varepsilon_t.$$

Thus $\phi_{kk}$ measures the additional contribution of $X_{t-k}$ when all lags 1 through $k-1$ are already in the model.

The theoretical PACF can be obtained by solving the Yule–Walker equations for successive autoregressive orders. In practice, the sample PACF is computed recursively via the Durbin–Levinson recursion:

$$\hat{\phi}_{11} = \hat{\rho}_1,$$ $$\hat{\phi}_{kk} = \frac{\hat{\rho}_k - \sum_{j=1}^{k-1} \hat{\phi}_{k-1,j}\,\hat{\rho}_{k-j}}{1 - \sum_{j=1}^{k-1} \hat{\phi}_{k-1,j}\,\hat{\rho}_j}, \quad k\ge 2,$$

with the intermediate coefficients updated as $\hat{\phi}_{kj} = \hat{\phi}_{k-1,j} - \hat{\phi}_{kk}\,\hat{\phi}_{k-1,k-j}$ for $j=1,\dots,k-1$.

Interpretation of the PACF complements that of the ACF:

Model	ACF	PACF
$\mathrm{AR}(p)$	Tails off (exponential/sinusoidal decay)	Cuts off after lag $p$
$\mathrm{MA}(q)$	Cuts off after lag $q$	Tails off
$\mathrm{ARMA}(p,q)$	Tails off	Tails off

Visualizing ACF and PACF with Correlograms in Python

A correlogram is the paired visualization of the sample ACF and PACF as bar plots against lag, overlaid with the confidence envelope. The statsmodels library provides ready-made functions:

import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

df = pd.read_csv("temp.csv", header=0, index_col=0, parse_dates=True)

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
plot_acf(df, lags=40, alpha=0.05, ax=axes[0])
plot_pacf(df, lags=40, alpha=0.05, ax=axes[1], method="ywm")
axes[0].set_title("ACF — Minimum Daily Temperature")
axes[1].set_title("PACF — Minimum Daily Temperature")
plt.tight_layout()
plt.show()

The lags argument controls the maximum horizon inspected, while alpha sets the confidence level (typically $0.05$ for $95\%$ bands). The method="ywm" flag requests Yule–Walker estimates with sample-mean adjustment, the standard choice for PACF computation. The Ljung–Box Q-statistic,

$$Q = T(T+2)\sum_{k=1}^{h}\frac{\hat{\rho}(k)^2}{T-k},$$

provides a joint test of the null that the first $h$ autocorrelations are simultaneously zero; rejection indicates that residual structure remains to be modeled.

Foundations: Stationarity and Identification Tools

Before modelling, we must confront the notion of stationarity. A time series $\{y_t\}$ is strictly stationary if its joint distribution is invariant under time shifts. In practice, we rely on weak stationarity: constant mean $\mathbb{E}[y_t] = \mu$, constant variance $\text{Var}(y_t) = \sigma^2$, and covariance $\text{Cov}(y_t, y_{t-k})$ that depends only on the lag $k$. Most classical linear models assume weak stationarity; non-stationary behaviour must be removed prior to modelling.

To test for a unit root, we employ the Augmented Dickey-Fuller (ADF) test. The null hypothesis $H_0$: the series has a unit root (is non-stationary) versus the alternative $H_1$: stationarity. A $p$-value below $0.05$ rejects $H_0$, indicating that no further differencing is required. If the $p$-value exceeds the threshold, the series is differenced once, $\nabla y_t = y_t - y_{t-1}$, and the test is reapplied.

Autoregressive (AR) Models

An autoregressive model of order $p$, AR($p$), expresses the current observation as a linear combination of its $p$ previous values plus a white-noise shock:

$$X_t = c + \sum_{i=1}^p \phi_i X_{t-i} + \varepsilon_t,$$

where $\varepsilon_t \sim \text{WN}(0,\sigma_\varepsilon^2)$ is a serially uncorrelated error term with zero mean and constant variance. Using the lag operator $L$ defined by $L^k X_t = X_{t-k}$, the AR($p$) model can be compactly written as $\Phi(L)X_t = c + \varepsilon_t$, where $\Phi(L) = 1 - \sum_{i=1}^p \phi_i L^i$ is the autoregressive polynomial.

The constant $c$ is related to the unconditional mean $\mu$ by $c = \mu(1-\phi_1-\cdots-\phi_p)$. Stationarity requires that all roots of the characteristic equation $1 - \phi_1 z - \phi_2 z^2 - \cdots - \phi_p z^p = 0$ lie outside the unit circle. For $p=1$, this reduces to $|\phi_1| < 1$. The interpretation is intuitive: the impact of a shock $\varepsilon_t$ decays geometrically as it propagates through future values.

Economic Intuition: In financial and macroeconomic contexts, AR models capture momentum and persistence. An AR(1) model with $\phi_1 > 0$ implies that a positive shock to GDP growth will positively influence subsequent quarters, reflecting the inertia in economic expansion.

Moving Average (MA) Models

It is crucial to distinguish a Moving Average (MA) model from the simple moving average smoothing technique used in descriptive analytics. An MA model of order $q$ asserts that the current observation is a linear combination of the current and past white noise error terms:

$$X_t = \mu + \varepsilon_t + \sum_{i=1}^q \theta_i \varepsilon_{t-i},$$

where $\mu$ is the mean of the series and $\theta_i$ are the moving average parameters. In lag operator notation, this is $X_t = \mu + \Theta(L)\varepsilon_t$, with $\Theta(L) = 1 + \sum_{i=1}^q \theta_i L^i$.

Economic Intuition: MA models represent the lingering effects of transient, unobserved shocks. An MA(1) process implies that an unexpected policy change or market shock ($\varepsilon_{t-1}$) affects the current period $X_t$ but completely dissipates by period $t+1$. This "finite memory" of shocks is highly applicable in modeling inventory adjustments or short-term market overreactions. An MA($q$) process is always stationary, but it must satisfy the invertibility condition (roots of $\Theta(z)=0$ outside the unit circle) to be expressed as an infinite AR process.

ARMA Models: Unifying AR and MA

When both autoregressive and moving-average dynamics are present, we obtain the ARMA($p,q$) model:

$$X_t = c + \sum_{i=1}^p \phi_i X_{t-i} + \varepsilon_t + \sum_{j=1}^q \theta_j \varepsilon_{t-j}.$$

ARMA models provide a parsimonious description of stationary time series. They are a cornerstone of the Box–Jenkins methodology, which proceeds in three iterative steps:

1. Identification: examine ACF and PACF to propose plausible orders $(p,q)$.
2. Estimation: fit the model via maximum likelihood (or conditional least squares).
3. Diagnostic checking: verify that residuals resemble white noise; if not, refine the model.

ARIMA: Accommodating Non-Stationarity via Integration

When the series itself is not stationary but its $d$-th difference is, we use an Autoregressive Integrated Moving Average model, denoted ARIMA($p,d,q$). Introducing the difference operator $\nabla = 1 - L$, where $L y_t = y_{t-1}$ is the backshift operator, the model can be written as

$$\Phi(L) (1-L)^d X_t = c + \Theta(L) \varepsilon_t.$$

The term $(1-L)^d$ removes $d$ unit roots, rendering the differenced series stationary. Most economic and business time series require $d=1$ or, occasionally, $d=2$.

Economic Intuition: The integration component models random walks, which are foundational in the Efficient Market Hypothesis. A series with $d=1$ implies that shocks have a permanent effect on the level of the series, characteristic of stock prices or aggregate price levels, whereas the differenced series (returns or inflation) remains stationary.

Implementation in Python:

from statsmodels.tsa.arima.model import ARIMA
from sklearn.metrics import mean_squared_error

# Fit MA(1) model: order = (p,d,q) = (0,0,1)
model = ARIMA(traindata, order=(0,0,1))
ma_fit = model.fit()
predictions = ma_fit.predict(start=len(traindata),
                             end=len(traindata)+len(testdata)-1,
                             dynamic=False)
error = mean_squared_error(testdata, predictions)

The dynamic=False argument ensures one-step-ahead forecasts that use actual lagged values where available.

Seasonal ARIMA (SARIMA) and SARIMAX

Many series exhibit repeating seasonal patterns. SARIMA extends ARIMA by incorporating seasonal differencing and seasonal AR and MA components. The general seasonal model is denoted SARIMA$(p,d,q)(P,D,Q)_s$, where $s$ is the seasonal period (e.g., 12 for monthly data). The formulation with backshift operators is:

$$\Phi(B^s) \phi(B) (1-B^s)^D (1-B)^d y_t = c + \Theta(B^s) \theta(B) \varepsilon_t,$$

with seasonal AR polynomial $\Phi(B^s) = 1 - \Phi_1 B^s - \cdots - \Phi_P B^{sP}$ and seasonal MA polynomial $\Theta(B^s) = 1 + \Theta_1 B^s + \cdots + \Theta_Q B^{sQ}$.

The SARIMAX model extends SARIMA by incorporating exogenous regressors:

$$\Phi(B^s) \phi(B) (1-B^s)^D (1-B)^d y_t = c + \beta Z_t + \Theta(B^s) \theta(B) \varepsilon_t.$$

This bridges the gap between pure time series analysis and regression, allowing analysts to control for known external interventions such as holiday effects, marketing campaigns, or macroeconomic indicators.

import statsmodels.api as sm
from statsmodels.tsa.stattools import adfuller

# 1. Stationarity Check
adf_result = adfuller(series.dropna())
print(f'ADF Statistic: {adf_result[0]:.4f}, p-value: {adf_result[1]:.4f}')

# 2. Train-Test Split
train_size = int(len(series) * 0.8)
train, test = series[:train_size], series[train_size:]

# 3. Fit SARIMA Model
model = sm.tsa.SARIMAX(train, order=(1, 1, 1), seasonal_order=(1, 1, 1, 12))
results = model.fit(disp=False)
print(results.summary())

# 4. Forecasting and Visualization
forecast = results.get_forecast(steps=len(test))
pred_mean = forecast.predicted_mean
conf_int = forecast.conf_int(alpha=0.05)

plt.figure(figsize=(12, 6))
plt.plot(train.index, train, label='Training Data', color='blue')
plt.plot(test.index, test, label='Actual Test Data', color='green')
plt.plot(pred_mean.index, pred_mean, label='SARIMA Forecast', color='red', linewidth=2)
plt.fill_between(conf_int.index, conf_int.iloc[:, 0], conf_int.iloc[:, 1],
                 color='pink', alpha=0.3, label='95% Confidence Interval')
plt.title('SARIMA Forecast with Prediction Intervals')
plt.xlabel('Time')
plt.ylabel('Value')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

The Vector Autoregressive (VAR) Foundation

The VAR($p$) model for a $K$-dimensional vector $y_t = (y_{1t}, y_{2t}, \dots, y_{Kt})'$ is written as

$$y_t = c + \Phi_1 y_{t-1} + \Phi_2 y_{t-2} + \dots + \Phi_p y_{t-p} + \varepsilon_t,$$

where $c$ is a $K \times 1$ vector of intercepts, each $\Phi_i$ is a $K \times K$ coefficient matrix, and $\varepsilon_t \sim \text{WN}(0, \Sigma)$ is a $K$-dimensional white noise process. The covariance matrix $\Sigma$ is positive definite, allowing for instantaneous correlations across equations — a crucial feature that distinguishes VARs from separate univariate regressions. Using the lag operator, the compact polynomial form is $\Phi(L) y_t = c + \varepsilon_t$, with $\Phi(L) = I - \Phi_1 L - \dots - \Phi_p L^p$.

The model is stable (stationary) if all roots of the determinantal equation $\det(\Phi(z)) = 0$ lie outside the unit circle. Under stability, the process has a causal infinite-order VMA($\infty$) representation.

Moving Average Components: The VARMA Model

The VARMA($p, q$) model adds moving average terms:

$$y_t = c + \sum_{i=1}^{p} \Phi_i y_{t-i} + \varepsilon_t + \sum_{j=1}^{q} \Theta_j \varepsilon_{t-j},$$

or in lag operator notation, $\Phi(L) y_t = c + \Theta(L) \varepsilon_t$, where $\Theta(L) = I + \Theta_1 L + \dots + \Theta_q L^q$. The process is invertible if all roots of $\det(\Theta(z)) = 0$ lie outside the unit circle. The identification of VARMA models is more delicate than in the univariate case because the matrices $\Phi(\cdot)$ and $\Theta(\cdot)$ can be multiplied by any non-singular rotation without changing the second-order properties. Standard practice imposes a canonical form, such as the echelon form.

Estimation of VARMA Models

Estimation typically proceeds via maximum likelihood under the assumption of Gaussian innovations. The conditional Gaussian log-likelihood is

$$\ell(\theta | y_1,\dots,y_T) = -\frac{TK}{2}\ln(2\pi) - \frac{T}{2}\ln|\Sigma| - \frac{1}{2}\sum_{t=p+1}^{T} \varepsilon_t'\Sigma^{-1}\varepsilon_t,$$

where $\varepsilon_t$ are computed recursively given the parameter matrices. For pure VAR models, equation-by-equation OLS is equivalent to multivariate least squares and provides consistent estimates, but for VARMA the likelihood is inherently non-linear.

Model Selection and Diagnostic Checking

Choosing the orders $(p,q)$ requires a balance between fit and parsimony. Information criteria are computed over a grid of plausible orders:

$$\text{AIC} = -2\ell_{\text{max}} + 2m,\quad \text{BIC} = -2\ell_{\text{max}} + m\ln T,$$

where $m$ is the total number of estimated parameters. After fitting, the residuals $\hat{\varepsilon}_t$ should resemble white noise. Diagnostics include multivariate portmanteau tests such as the Ljung–Box on residual cross-correlations:

$$Q(m) = T(T+2) \sum_{h=1}^{m} \frac{1}{T-h} \text{tr}(\hat{\mathbf{C}}_h' \hat{\mathbf{C}}_0^{-1} \hat{\mathbf{C}}_h \hat{\mathbf{C}}_0^{-1}),$$

where $\hat{\mathbf{C}}_h$ is the sample autocovariance matrix at lag $h$.

Forecasting with VARMA Models

Once an adequate VARMA$(p,q)$ has been estimated, optimal $h$-step-ahead forecasts are obtained by iterating the dynamic equations while setting future innovations to their expected value of zero:

$$y_{T+h|T} = \hat{c} + \sum_{i=1}^{p} \hat{\Phi}_i y_{T+h-i|T} + \sum_{j=h}^{q} \hat{\Theta}_j \hat{\varepsilon}_{T+h-j}.$$

A critical tool for interpretation is the forecast error variance decomposition (FEVD), which quantifies the proportion of the $h$-step forecast variance of each variable attributable to its own innovations versus those of other variables.

Incorporating Exogenous Regressors: VARMAX

The VARMAX$(p,q,b)$ model extends VARMA by including current and lagged values of a $d$-dimensional vector of exogenous variables $x_t$:

$$y_t = c + \sum_{i=1}^{p} \Phi_i y_{t-i} + \sum_{i=1}^{q} \Theta_i \varepsilon_{t-i} + \sum_{i=0}^{b} \Beta_i x_{t-i} + \varepsilon_t.$$

Here each $\Beta_i$ is a $K \times d$ matrix of coefficients. The term $x_t$ is assumed fully predetermined or strictly exogenous. In financial forecasting, a VARMAX model might predict stock returns and trading volume (endogenous variables) using macroeconomic indicators like unemployment rates or monetary policy decisions (exogenous variables).

Visualization Strategies for Multivariate Forecasts

1. Impulse Response Functions (IRFs): Plot the dynamic response of each variable to a one-standard-deviation shock in another variable. IRFs reveal the propagation mechanism and persistence of shocks through the system.
2. Forecast Error Variance Decomposition (FEVD): Visualize the proportion of forecast error variance for each variable attributable to shocks from other variables.
3. Multi-panel forecast plots: Display actual versus predicted values for all variables in a grid layout, with confidence bands indicating forecast uncertainty.
4. Coefficient heatmaps: For large systems, visualize the AR and MA coefficient matrices as heatmaps to identify strong interdependencies.

from statsmodels.tsa.statespace.varmax import VARMAX

# Fit VARMA(1,1) model
model = VARMAX(endog=data, order=(1, 1), trend='c')
results = model.fit(disp=False)

# Forecast
forecast = results.forecast(steps=12)

# Impulse Response Analysis
irf = results.impulse_responses(steps=20, orthogonalized=True)

Exponential Smoothing constitutes a foundational family of time series forecasting methods that generate predictions by applying exponentially decreasing weights to past observations. Unlike simple moving averages, which assign equal weight to a fixed window of historical data, exponential smoothing models operate on the premise that recent observations contain more predictive information than distant ones.

Simple Exponential Smoothing (SES)

Simple Exponential Smoothing is the most basic formulation, designed strictly for time series data that exhibit no discernible trend or seasonality. The model relies on a single parameter, the smoothing factor $\alpha$, which governs the rate at which the influence of past observations decays.

Let $y_t$ denote the observed value at time $t$. The smoothed statistic $s_t$, which represents the level of the series, is updated recursively:

$$s_t = \alpha y_t + (1 - \alpha) s_{t-1},$$

where $0 < \alpha < 1$. To understand the economic and statistical intuition behind this formulation, we can expand the recursive equation through backward substitution:

$$s_t = \alpha y_t + \alpha(1 - \alpha) y_{t-1} + \alpha(1 - \alpha)^2 y_{t-2} + \dots$$

By induction, assuming the process extends infinitely into the past, the smoothed value is an infinite weighted sum of all past observations:

$$s_t = \alpha \sum_{j=0}^{\infty} (1 - \alpha)^j y_{t-j}.$$

Thus, the weight assigned to observation $y_{t-k}$ is $\alpha(1-\alpha)^k$, which decays geometrically with lag $k$. When $\alpha$ is close to 1, the series adapts quickly to new observations; a smaller $\alpha$ produces a smoother level that reacts sluggishly to transient shocks. The forecast for all future periods is constant and equals the most recent level:

$$\hat{y}_{T+h|T} = s_T, \qquad h = 1,2,\dots$$

The parameter $\alpha$ is usually chosen by minimizing the sum of squared one-step-ahead forecast errors over the training set. SES can be cast as a special case of the ARIMA(0,1,1) model, which provides a statistical framework for inference.

Python Implementation:

import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.api import SimpleExpSmoothing

# Sample data
data = [41.7275, 24.0418, 32.3281, 37.3287, 46.2132,
        29.3463, 36.4829, 42.9777, 48.9015, 31.1802,
        37.7179, 40.4202, 51.2069, 31.8872, 40.9783,
        43.7725, 55.5586, 33.8509, 42.0764, 45.6423,
        59.7668, 35.1919, 44.3197, 47.9137]
index = pd.date_range(start="2005", end="2010-Q4", freq="QS-OCT")
series = pd.Series(data, index)

train = series.iloc[:-5]
test  = series.iloc[-5:]

model = SimpleExpSmoothing(train, initialization_method="heuristic")
ses_fit = model.fit(smoothing_level=0.2, optimized=False)
forecast = ses_fit.forecast(steps=5)

plt.figure(figsize=(10,6))
plt.plot(train, label='Training data', color='black')
plt.plot(test, label='Test data', color='blue', marker='o')
plt.plot(forecast, label='SES Forecast (α=0.2)', color='red', linestyle='--', marker='s')
plt.title('Simple Exponential Smoothing')
plt.legend()
plt.xlabel('Date')
plt.ylabel('Value')
plt.grid(alpha=0.3)
plt.show()

Holt's Linear Trend Method

When a time series exhibits a clear trend, SES is inadequate because its forecasts remain constant indefinitely. extended the method by adding a second equation to capture a (local) linear trend:

$$\ell_t = \alpha y_t + (1-\alpha)(\ell_{t-1} + b_{t-1}), \quad\text{(Level)}$$ $$b_t = \beta^*(\ell_t - \ell_{t-1}) + (1-\beta^*)b_{t-1}, \quad\text{(Trend)}$$

where $\alpha \in (0,1)$ controls the level adaptation, and $\beta^* \in (0,1)$ controls the trend adaptation. The $h$-step-ahead forecast is:

$$\hat{y}_{T+h|T} = \ell_T + h \, b_T.$$

A damped trend variant, which prevents the forecast from growing indefinitely, multiplies the trend term by a damping factor $\phi \in (0,1)$, yielding $\hat{y}_{T+h|T} = \ell_T + \sum_{i=1}^{h} \phi^i \, b_T$. This is often more robust in practice, especially for long-range forecasts.

Holt-Winters Seasonal Method

The Holt-Winters method extends Holt's formulation by incorporating a seasonal component. The model exists in two variants depending on the nature of the seasonality:

Additive Holt-Winters (seasonal fluctuations are roughly constant over time):

$$\ell_t = \alpha (y_t - s_{t-m}) + (1-\alpha)(\ell_{t-1} + b_{t-1}),$$ $$b_t = \beta^* (\ell_t - \ell_{t-1}) + (1-\beta^*) b_{t-1},$$ $$s_t = \gamma (y_t - \ell_{t-1} - b_{t-1}) + (1-\gamma) s_{t-m}.$$

The forecast for $h$ periods ahead is $\hat{y}_{T+h|T} = \ell_T + h\, b_T + s_{T+h-m(k+1)}$, where $k = \lfloor (h-1)/m \rfloor$ ensures the most recent relevant seasonal factor is used.

Multiplicative Holt-Winters (seasonal variation scales with level):

$$\ell_t = \alpha \left(\frac{y_t}{s_{t-m}}\right) + (1-\alpha)(\ell_{t-1} + b_{t-1}),$$ $$b_t = \beta^* (\ell_t - \ell_{t-1}) + (1-\beta^*) b_{t-1},$$ $$s_t = \gamma \left(\frac{y_t}{\ell_{t-1} + b_{t-1}}\right) + (1-\gamma) s_{t-m}.$$

The multiplicative forecast is $\hat{y}_{T+h|T} = (\ell_T + h b_T) s_{T+h-m(k+1)}$.

The ETS Framework and Prediction Intervals

Modern implementations of exponential smoothing are grounded in the ETS (Error, Trend, Seasonal) state-space framework. This statistical formulation not only unifies the various exponential smoothing methods but also provides a rigorous basis for calculating prediction intervals. By modeling the error term explicitly (either additive or multiplicative), the ETS framework allows analysts to quantify forecast uncertainty, which widens as the forecast horizon $h$ increases.

from statsmodels.tsa.holtwinters import ExponentialSmoothing

# Fit Multiplicative Holt-Winters model with damped trend
hw_model = ExponentialSmoothing(
    train_hw, 
    trend='add', 
    seasonal='mul', 
    seasonal_periods=12,
    damped_trend=True,
    initialization_method="estimated"
)
hw_fit = hw_model.fit(optimized=True)

# Forecast and extract prediction intervals
hw_forecast = hw_fit.forecast(steps=12)
pred_obj = hw_fit.get_prediction(start=train_hw.index[0], end=test_hw.index[-1])
pred_summary = pred_obj.summary_frame(alpha=0.05)

fig, ax = plt.subplots(figsize=(14, 7))
ax.plot(train_hw.index, train_hw, label='Training Data', color='black', lw=1.5)
ax.plot(test_hw.index, test_hw, label='Actual Test Data', color='blue', marker='o')
ax.plot(test_hw.index, hw_forecast, label='Holt-Winters Forecast', color='red', lw=2)
ax.fill_between(test_hw.index, 
                pred_summary['pi_lower'].iloc[-12:], 
                pred_summary['pi_upper'].iloc[-12:],
                color='pink', alpha=0.3, label='95% Prediction Interval')
ax.set_title("Holt-Winters Exponential Smoothing with Prediction Intervals")
ax.legend()
plt.tight_layout()
plt.show()