🎯 What You'll Learn
- Understand what makes time series data fundamentally different from tabular data — temporal ordering, autocorrelation, and non-i.i.d. structure
- Identify the four components of a time series: trend, seasonality, cyclicality, and noise
- Test for stationarity using the Augmented Dickey-Fuller test and make series stationary through differencing and transformations
- Decompose a time series into its components using
seasonal_decomposefrom statsmodels - Interpret ACF and PACF plots to identify autocorrelation patterns and seasonal periodicity
1 What Is Time Series Data?
A time series is a sequence of observations collected at successive, typically equally-spaced, points in time. The defining characteristic that separates time series from tabular data is that the order of observations matters. The value at time t depends on values at times t−1, t−2, and so on — a property called autocorrelation. This violates the fundamental i.i.d. (independent and identically distributed) assumption that underpins most standard machine learning algorithms.
Examples of time series data that appear throughout ML and data science:
- Finance: stock prices, exchange rates, volatility indices, trading volumes
- Operations: website traffic, API request counts, server CPU/memory utilization
- Retail: daily/weekly/monthly sales, inventory levels, customer counts
- Climate/IoT: temperature, rainfall, energy consumption, sensor readings
- Health: heart rate, EEG signals, blood glucose, hospital admissions
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
# Loading time series data with pandas
# The key: parse_dates and index_col together create a DatetimeIndex
# df = pd.read_csv('sales_data.csv', parse_dates=['date'], index_col='date')
# Generate synthetic monthly sales data to illustrate concepts
np.random.seed(42)
dates = pd.date_range(start='2019-01-01', periods=60, freq='MS') # 5 years, monthly
months = dates.month
# Trend + seasonality + noise
trend = np.linspace(1000, 1800, 60) # growing trend
seasonal = 200 * np.sin(2 * np.pi * (months - 3) / 12) # annual season (peak in summer)
noise = np.random.normal(0, 80, 60)
sales = trend + seasonal + noise
ts = pd.Series(sales, index=dates, name='monthly_sales')
print(f"Series type: {type(ts)}")
print(f"Index type: {type(ts.index)}")
print(f"Frequency: {ts.index.freq}")
print(f"Date range: {ts.index[0].date()} to {ts.index[-1].date()}")
print(f"Observations: {len(ts)}")
print(f"\nFirst 6 observations:")
print(ts.head(6).round(1))
# Key pandas operations for time series
print(f"\nMonthly stats (last 12 months):")
print(ts[-12:].describe().round(1))
# Resample to quarterly (useful for higher-frequency data)
quarterly = ts.resample('QS').sum()
print(f"\nQuarterly total sales (last 4 quarters):")
print(quarterly[-4:].round(0))Standard machine learning algorithms (logistic regression, random forest, SVMs) assume that training samples are drawn independently from the same distribution. Time series observations are neither independent (adjacent values are correlated) nor identically distributed (the distribution may change over time with trend and seasonality). This means: (1) standard cross-validation is wrong — you must never use future data to predict the past; (2) standard feature importance metrics may be misleading; (3) models that ignore temporal structure often underperform specialized time series methods.
2 Components of a Time Series
Every time series can be thought of as a combination of underlying components. Understanding which components are present is the essential first step before choosing a modeling or forecasting approach.
The Four Components
- Trend (T): the long-term direction of the series — upward, downward, or flat. Example: the steady growth in e-commerce sales over 10 years. Trends can be linear, exponential, or non-monotonic.
- Seasonality (S): regular, periodic fluctuations with a fixed and known frequency. Example: retail sales spike every December; website traffic drops every weekend; energy consumption peaks every summer. The period is fixed — 12 months for annual, 7 days for weekly, 24 hours for daily.
- Cyclicality (C): longer-period fluctuations that are irregular and do not have a fixed frequency. Example: economic boom-and-bust cycles, housing market cycles. Cyclicality is distinguished from seasonality by its lack of a fixed period — a recession might last 2 years, 5 years, or 10 years.
- Residuals/Noise (R): the random, unexplained variation that remains after accounting for trend, seasonality, and cycles. Ideally, a well-modeled time series has residuals that look like white noise — no remaining patterns.
Additive vs Multiplicative Models
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
np.random.seed(42)
# Additive model: y = T + S + R (use when seasonal amplitude is constant)
# The seasonal swing stays ~200 units regardless of the trend level
dates = pd.date_range('2018-01-01', periods=72, freq='MS')
trend_additive = np.linspace(500, 900, 72)
seasonal_additive = 150 * np.sin(2 * np.pi * np.arange(72) / 12)
noise = np.random.normal(0, 30, 72)
y_additive = trend_additive + seasonal_additive + noise
# Multiplicative model: y = T × S × R (use when seasonal amplitude grows with trend)
# The seasonal swing is a FRACTION of the trend — grows as trend grows
trend_mult = np.linspace(100, 500, 72)
seasonal_mult = 1 + 0.3 * np.sin(2 * np.pi * np.arange(72) / 12)
noise_mult = np.random.normal(1, 0.05, 72)
y_multiplicative = trend_mult * seasonal_mult * noise_mult
ts_additive = pd.Series(y_additive, index=dates)
ts_mult = pd.Series(y_multiplicative, index=dates)
print("Additive model: seasonal amplitude is CONSTANT regardless of trend")
print(f" Early seasonal range: {ts_additive[:12].max()-ts_additive[:12].min():.0f}")
print(f" Late seasonal range: {ts_additive[-12:].max()-ts_additive[-12:].min():.0f}")
print()
print("Multiplicative model: seasonal amplitude GROWS with trend")
print(f" Early seasonal range: {ts_mult[:12].max()-ts_mult[:12].min():.0f}")
print(f" Late seasonal range: {ts_mult[-12:].max()-ts_mult[-12:].min():.0f}")Plot your time series and look at the seasonal fluctuations over time. If the peaks and troughs stay roughly the same height regardless of the trend level — use an additive model. If the swings get bigger as the series grows (the peaks are proportionally higher, troughs proportionally lower) — use a multiplicative model. A log transform converts a multiplicative series into an additive one: log(T × S × R) = log(T) + log(S) + log(R). This is one of the most useful preprocessing tricks for time series.
3 Visualizing Time Series
Before any modeling, thorough visual exploration reveals structure, outliers, and patterns that guide all subsequent decisions. Key visualization techniques include line plots, rolling statistics, and lag plots.
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
# Reconstruct synthetic monthly sales series
dates = pd.date_range(start='2019-01-01', periods=60, freq='MS')
months = np.arange(1, 61) % 12 + 1
trend = np.linspace(1000, 1800, 60)
seasonal = 200 * np.sin(2 * np.pi * (months - 3) / 12)
noise = np.random.normal(0, 80, 60)
ts = pd.Series(trend + seasonal + noise, index=dates, name='sales')
# 1. Rolling mean and std — reveals trend and changing variance
window = 12 # 12-month rolling window
rolling_mean = ts.rolling(window=window, center=True).mean()
rolling_std = ts.rolling(window=window, center=True).std()
print("Rolling statistics (last 6 months):")
comparison = pd.DataFrame({
'actual': ts[-6:].round(1),
'rolling_mean': rolling_mean[-6:].round(1),
'rolling_std': rolling_std[-6:].round(1),
})
print(comparison)
# 2. Monthly seasonal pattern (seasonal subseries)
ts_df = ts.to_frame()
ts_df['month'] = ts_df.index.month
monthly_avg = ts_df.groupby('month')['sales'].mean()
print("\nAverage sales by month (seasonal pattern):")
month_names = ['Jan','Feb','Mar','Apr','May','Jun',
'Jul','Aug','Sep','Oct','Nov','Dec']
for m, avg in enumerate(monthly_avg, 1):
bar = '█' * int((avg - 900) / 30)
print(f" {month_names[m-1]}: {avg:7.1f} {bar}")
# 3. Year-over-year comparison
ts_df['year'] = ts_df.index.year
print("\nMean annual sales by year (shows trend):")
print(ts_df.groupby('year')['sales'].mean().round(1))
# 4. Lag plot — check for autocorrelation
lag1_corr = ts.corr(ts.shift(1))
lag12_corr = ts.corr(ts.shift(12))
print(f"\nLag-1 autocorrelation (adjacent months): {lag1_corr:.4f}")
print(f"Lag-12 autocorrelation (same month, prev year): {lag12_corr:.4f}")
print("Strong lag-12 correlation confirms annual seasonality")4 Stationarity
A time series is stationary if its statistical properties — mean, variance, and autocorrelation structure — are constant over time. This is not just a theoretical nicety: virtually every classical forecasting model (ARIMA, VAR, exponential smoothing) requires or implicitly assumes stationarity. Applying these models to non-stationary data produces unreliable forecasts and misleading statistical tests.
Stationarity has a formal definition: a series {Yₜ} is weakly stationary if:
- E[Yₜ] = μ (constant mean)
- Var(Yₜ) = σ² (constant variance)
- Cov(Yₜ, Yₜ₋ₖ) = γ(k) — only depends on lag k, not on time t
import pandas as pd
import numpy as np
from statsmodels.tsa.stattools import adfuller
np.random.seed(42)
dates = pd.date_range('2019-01-01', periods=60, freq='MS')
months = np.arange(1, 61) % 12 + 1
trend = np.linspace(1000, 1800, 60)
seasonal = 200 * np.sin(2 * np.pi * (months - 3) / 12)
noise = np.random.normal(0, 80, 60)
ts = pd.Series(trend + seasonal + noise, index=dates, name='sales')
def adf_test(series, series_name="Series"):
"""Run ADF test and print a clean interpretation."""
result = adfuller(series.dropna(), autolag='AIC')
adf_stat, p_value, lags_used, n_obs, crit_values, _ = result
print(f"\n=== ADF Test: {series_name} ===")
print(f"ADF Statistic: {adf_stat:.4f}")
print(f"p-value: {p_value:.6f}")
print(f"Lags used: {lags_used}")
print(f"Num observations:{n_obs}")
print("Critical values:")
for key, val in crit_values.items():
print(f" {key}: {val:.4f}")
if p_value < 0.05:
print("✅ Reject H0: series IS stationary (p < 0.05)")
else:
print("❌ Fail to reject H0: series IS NOT stationary (p >= 0.05)")
print(" → Apply differencing or log transform")
return p_value
# Test 1: Original series (has trend + seasonality → non-stationary)
p1 = adf_test(ts, "Original Sales (with trend + seasonality)")
# Test 2: First difference (removes trend)
ts_diff = ts.diff().dropna()
p2 = adf_test(ts_diff, "First Difference")
# Test 3: Seasonal difference (lag-12, removes annual seasonality)
ts_seasonal_diff = ts.diff(12).dropna()
p3 = adf_test(ts_seasonal_diff, "Seasonal Difference (lag-12)")
print(f"\nSummary:")
print(f" Original: p = {p1:.4f} {'stationary' if p1 < 0.05 else 'NON-stationary'}")
print(f" First diff: p = {p2:.4f} {'stationary' if p2 < 0.05 else 'NON-stationary'}")
print(f" Seasonal diff: p = {p3:.4f} {'stationary' if p3 < 0.05 else 'NON-stationary'}")The Augmented Dickey-Fuller (ADF) test's null hypothesis is that the series has a "unit root" — meaning it is non-stationary. A p-value < 0.05 gives you confidence to reject this null hypothesis, meaning the series is stationary. A p-value ≥ 0.05 means you cannot reject non-stationarity — you need to transform the data. Do not blindly trust the ADF test alone; always combine with visual inspection of rolling mean and rolling std. A visually obvious trend that the ADF test fails to detect (because it needs more data) should still be differenced.
5 Making a Series Stationary
When the ADF test confirms non-stationarity, there are several standard transformations to apply, in order of severity:
import pandas as pd
import numpy as np
from statsmodels.tsa.stattools import adfuller
np.random.seed(42)
dates = pd.date_range('2019-01-01', periods=60, freq='MS')
months = np.arange(1, 61) % 12 + 1
trend = np.linspace(1000, 1800, 60)
seasonal = 200 * np.sin(2 * np.pi * (months - 3) / 12)
noise = np.random.normal(0, 80, 60)
ts = pd.Series(trend + seasonal + noise, index=dates, name='sales')
def quick_adf(series):
p = adfuller(series.dropna(), autolag='AIC')[1]
return f"p={p:.4f} ({'✅ stationary' if p < 0.05 else '❌ not stationary'})"
# Transformation 1: First differencing — removes linear trend
# Interpretation: "month-over-month change"
ts_diff1 = ts.diff(1)
print(f"First difference (Δ1): {quick_adf(ts_diff1)}")
# Transformation 2: Seasonal differencing (lag=12) — removes annual seasonality
# Interpretation: "same month vs same month last year"
ts_sdiff = ts.diff(12)
print(f"Seasonal diff (Δ12): {quick_adf(ts_sdiff)}")
# Transformation 3: Both first + seasonal differencing (for trend + seasonality)
ts_combined = ts.diff(12).diff(1)
print(f"Seasonal + First diff: {quick_adf(ts_combined)}")
# Transformation 4: Log transform — stabilises growing variance (multiplicative to additive)
ts_log = np.log(ts)
ts_log_diff = ts_log.diff(1)
print(f"Log + First diff: {quick_adf(ts_log_diff)}")
# Transformation 5: Box-Cox transform — generalises log (lambda=0), sqrt (lambda=0.5)
from scipy.stats import boxcox
ts_boxcox, lambda_bc = boxcox(ts.values)
ts_bc = pd.Series(ts_boxcox, index=ts.index)
print(f"Box-Cox (λ={lambda_bc:.3f}) + diff: {quick_adf(ts_bc.diff(1))}")
print("\n=== Practical guidance ===")
print("1. Try first diff → check ADF → if still non-stationary:")
print("2. Try seasonal diff (if seasonal) → check ADF → if still non-stationary:")
print("3. Try log transform first, then diff")
print("4. Rarely need to diff more than twice (d=2); if so, suspect the data")
print("5. Always verify the differenced series looks like stationary white noise")6 Time Series Decomposition
Decomposition formally separates a time series into its trend, seasonal, and residual components. It serves two purposes: (1) understanding the data — how much of the variation is due to trend vs seasonality vs noise? and (2) improving forecasting — some methods forecast each component separately and then recombine.
import pandas as pd
import numpy as np
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.stattools import adfuller
np.random.seed(42)
dates = pd.date_range('2019-01-01', periods=60, freq='MS')
months = np.arange(1, 61) % 12 + 1
trend_component = np.linspace(1000, 1800, 60)
seasonal_component = 200 * np.sin(2 * np.pi * (months - 3) / 12)
noise = np.random.normal(0, 80, 60)
ts = pd.Series(trend_component + seasonal_component + noise, index=dates)
# Additive decomposition (seasonal amplitude is roughly constant)
result_add = seasonal_decompose(
ts,
model='additive', # or 'multiplicative'
period=12, # seasonality period in observations
extrapolate_trend='freq' # fill NaN at edges using trend frequency
)
# Access each component
trend_est = result_add.trend
seasonal_est = result_add.seasonal
residual_est = result_add.resid
print("=== Additive Decomposition Results ===")
print(f"\nTrend component (first 6 months):")
print(trend_est[:6].round(1))
print(f"\nSeasonal component (first 12 months = one full cycle):")
print(seasonal_est[:12].round(1))
print(f"\nResidual component (last 6 months):")
print(residual_est[-6:].round(1))
# Verify: trend + seasonal + residual ≈ original
reconstruction = trend_est + seasonal_est + residual_est
reconstruction_error = (ts - reconstruction).abs().mean()
print(f"\nReconstruction MAE: {reconstruction_error:.6f} (should be ~0)")
# Analyze each component
print(f"\n=== Component Analysis ===")
# Trend: what is the monthly growth rate?
trend_growth = (trend_est.iloc[-1] / trend_est.iloc[0] - 1) * 100
print(f"Total trend growth over period: {trend_growth:.1f}%")
monthly_growth = trend_est.pct_change().mean() * 100
print(f"Average monthly trend growth: {monthly_growth:.2f}%")
annual_growth = (1 + monthly_growth/100)**12 - 1
print(f"Implied annual growth rate: {annual_growth*100:.1f}%")
# Seasonality: which months peak/trough?
seasonal_by_month = seasonal_est[:12]
peak_month = seasonal_by_month.idxmax().strftime('%B')
trough_month = seasonal_by_month.idxmin().strftime('%B')
seasonal_range = seasonal_by_month.max() - seasonal_by_month.min()
print(f"\nSeasonal peak: {peak_month}")
print(f"Seasonal trough: {trough_month}")
print(f"Seasonal range: {seasonal_range:.1f} units")
# Residuals: should look like white noise
resid_adf_p = adfuller(residual_est.dropna(), autolag='AIC')[1]
print(f"\nResidual ADF p-value: {resid_adf_p:.4f}")
print(f"Residuals {'✅ stationary (good — no remaining pattern)' if resid_adf_p < 0.05 else '❌ NOT stationary (model missing a pattern)'}")
print(f"Residual mean: {residual_est.mean():.2f} (should be ~0)")
print(f"Residual std: {residual_est.std():.2f}")The print statements above tell the story in numbers — but decomposition is fundamentally a visual technique. Stacking the four series on a shared time axis makes it immediately obvious how much of the original signal is explained by trend, how much by the repeating annual cycle, and how much is left as unstructured noise:
Additive decomposition of the lesson's synthetic monthly sales series (5 years, upward trend + annual seasonality + noise). Original = Trend + Seasonal + Residual at every point. Hover any panel to read exact values; drag to zoom (double-click to reset).
After decomposition, inspect the residuals carefully. If they look like white noise (random scatter around zero, stationary, no visible patterns) — the decomposition has successfully captured the structure. If the residuals show clear patterns — spikes, changing variance, or autocorrelation — something systematic was missed. Possible causes: the trend is non-linear, there are multiple periodicities (daily + weekly + annual), or there are structural breaks (COVID, a product launch, a competitor exit). Good residuals mean a good model; patterned residuals mean opportunity for improvement.
7 Autocorrelation (ACF) and Partial Autocorrelation (PACF)
The ACF and PACF plots are the fundamental diagnostic tools for time series analysis. They tell you how a series relates to its own past values — information that guides the selection of AR and MA orders for ARIMA models (covered in Lesson 36).
- ACF (Autocorrelation Function): the correlation between the series and its lagged copy at each lag. A spike at lag k means the series at time t is correlated with the series at time t−k. Includes both direct and indirect effects of all intermediate lags.
- PACF (Partial Autocorrelation Function): the correlation at lag k after removing the effect of all shorter lags. Isolates the direct relationship between Yₜ and Yₜ₋ₖ, controlling for lags 1 through k−1.
import pandas as pd
import numpy as np
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.stattools import acf, pacf
np.random.seed(42)
dates = pd.date_range('2019-01-01', periods=60, freq='MS')
months = np.arange(1, 61) % 12 + 1
trend = np.linspace(1000, 1800, 60)
seasonal = 200 * np.sin(2 * np.pi * (months - 3) / 12)
noise = np.random.normal(0, 80, 60)
ts = pd.Series(trend + seasonal + noise, index=dates)
# Compute ACF values manually
nlags = 24 # show up to 24 lags (2 years for monthly data)
acf_values = acf(ts, nlags=nlags, fft=True)
pacf_values = pacf(ts, nlags=nlags)
print("ACF values (with confidence interval ~ ±2/√n):")
n = len(ts)
conf_band = 2 / np.sqrt(n)
print(f"Confidence band: ±{conf_band:.3f}")
print(f"\nLag ACF Significant?")
print("-" * 35)
for lag in [0, 1, 2, 3, 6, 11, 12, 13, 23, 24]:
acf_val = acf_values[lag]
sig = "✅ significant" if abs(acf_val) > conf_band else " not significant"
print(f" {lag:2d} {acf_val:+.4f} {sig}")
# ACF on first-differenced series (stationarized)
ts_diff = ts.diff().dropna()
acf_diff = acf(ts_diff, nlags=24, fft=True)
print(f"\nACF of first-differenced series:")
print(f" Lag 1: {acf_diff[1]:+.4f}")
print(f" Lag 12: {acf_diff[12]:+.4f}")
print(f" Lag 24: {acf_diff[24]:+.4f}")
print("\nReading the patterns:")
print(" ACF cuts off after lag q → MA(q) process")
print(" ACF decays gradually → AR process")
print(" ACF has spikes at 12, 24 → annual seasonality present")The pattern described in the printout — strong, slowly-decaying correlation with spikes that re-emerge every 12 lags — is much easier to spot visually than in a column of numbers. The chart below plots both functions for the same synthetic sales series, with the ±2/√n significance band shaded:
ACF of the raw sales series — slow decay plus spikes at lags 12 and 24 signal trend and annual seasonality. Bars outside the shaded band are statistically significant at the ±2/√n threshold.
# Practical guide to reading ACF / PACF for ARIMA order selection
print("=== ACF/PACF Pattern Recognition Guide ===\n")
patterns = {
"Pure AR(p)": {
"ACF": "Decays gradually (exponential or sinusoidal damping)",
"PACF": f"Cuts off sharply after lag p (spikes 1..p, then ~0)",
"Example": "Stock price returns often show AR(1) or AR(2) patterns"
},
"Pure MA(q)": {
"ACF": "Cuts off sharply after lag q (spikes 1..q, then ~0)",
"PACF": "Decays gradually",
"Example": "Moving average smoothed series; error autocorrelation"
},
"ARIMA(p,d,q)": {
"ACF": "Slow decay — apply differencing first, then look at ACF/PACF",
"PACF": "Slow decay — same as above",
"Example": "Most real-world economic series with trend"
},
"Seasonal AR(1)": {
"ACF": "Significant spikes at lags s, 2s, 3s (e.g. 12, 24, 36 for monthly)",
"PACF": "Single significant spike at lag s",
"Example": "Monthly retail sales with annual seasonality"
},
}
for model, info in patterns.items():
print(f"Model: {model}")
for k, v in info.items():
print(f" {k:8s}: {v}")
print()ACF and PACF are only interpretable on a stationary series. On a raw series with trend or seasonality, both ACF and PACF will show very high correlation at all lags — this isn't informative, it just reflects the non-stationarity. Always difference (and log-transform if needed) until the ADF test confirms stationarity, then examine the ACF and PACF of the transformed series to guide ARIMA order selection.
8 Rolling Statistics & Moving Averages
Rolling statistics are among the most practically useful tools in a time series analyst's toolkit. They smooth noise, reveal underlying trends, and detect changing volatility — all without any formal modeling assumptions.
import pandas as pd
import numpy as np
np.random.seed(42)
dates = pd.date_range('2019-01-01', periods=60, freq='MS')
months = np.arange(1, 61) % 12 + 1
trend = np.linspace(1000, 1800, 60)
seasonal = 200 * np.sin(2 * np.pi * (months - 3) / 12)
noise = np.random.normal(0, 80, 60)
ts = pd.Series(trend + seasonal + noise, index=dates, name='sales')
# ─── Simple Moving Average (SMA) ───
# Equal weight to all observations in the window
sma_3 = ts.rolling(window=3, center=False).mean() # 3-month SMA
sma_12 = ts.rolling(window=12, center=False).mean() # 12-month SMA (annual)
print("Simple Moving Averages — last 6 periods:")
print(pd.DataFrame({
'actual': ts[-6:].round(1),
'SMA-3': sma_3[-6:].round(1),
'SMA-12': sma_12[-6:].round(1),
}))
# ─── Exponentially Weighted Moving Average (EWMA) ───
# More weight to recent observations — adapts faster to changes
# span parameter: effectively the number of periods to give most weight to
ewma_3 = ts.ewm(span=3, adjust=False).mean()
ewma_12 = ts.ewm(span=12, adjust=False).mean()
print("\nEWMA vs SMA comparison (last 3 months):")
comparison = pd.DataFrame({
'actual': ts[-3:].round(1),
'SMA-3': sma_3[-3:].round(1),
'EWMA-3': ewma_3[-3:].round(1),
})
print(comparison)
print("EWMA-3 reacts faster to recent changes than SMA-3")
# ─── Rolling Standard Deviation — detecting changing volatility ───
rolling_std_12 = ts.rolling(window=12).std()
print(f"\nRolling 12-month standard deviation (shows heteroscedasticity):")
print(rolling_std_12[11::12].round(1)) # once per year
# ─── Detecting Heteroscedasticity ───
# If rolling_std grows over time → variance is NOT constant → series is heteroscedastic
# This suggests using a log transform before modeling
early_std = rolling_std_12[:30].mean()
late_std = rolling_std_12[30:].mean()
ratio = late_std / early_std
print(f"\nEarly period rolling std: {early_std:.1f}")
print(f"Late period rolling std: {late_std:.1f}")
print(f"Ratio (late/early): {ratio:.2f}")
print(f"{'⚠️ Heteroscedastic — consider log transform' if ratio > 1.2 else '✅ Roughly homoscedastic'}")Real-World Spotlight: Retail Sales Decomposition
A supermarket chain wants to understand its monthly sales pattern for staffing, inventory planning, and financial forecasting. We'll apply the full time series analysis pipeline: load data, visualize, test stationarity, decompose, and extract actionable business insights from each component.
import pandas as pd
import numpy as np
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.stattools import adfuller, acf
from statsmodels.graphics.tsaplots import plot_acf
np.random.seed(42)
# Simulate 4 years of monthly supermarket sales with realistic structure
n_months = 48
dates = pd.date_range('2020-01-01', periods=n_months, freq='MS')
# Underlying structure
annual_trend = 1.032 # 3.2% annual growth
base_sales = 2_000_000 # £2M monthly baseline
trend_component = np.array([
base_sales * (annual_trend ** (m/12)) for m in range(n_months)
])
# Multiplicative seasonal pattern: Christmas (Dec) is 2.4x average
month_factors = {
1: 0.82, 2: 0.78, 3: 0.88, 4: 0.92, 5: 0.96, 6: 0.98,
7: 0.99, 8: 1.00, 9: 0.97, 10: 1.02, 11: 1.18, 12: 2.40
}
months = dates.month
seasonal_component = np.array([month_factors[m] for m in months])
# Small noise
noise_component = np.random.normal(1.0, 0.03, n_months)
sales = trend_component * seasonal_component * noise_component
ts = pd.Series(sales / 1e6, index=dates, name='sales_millions') # in £M
print("=== Supermarket Monthly Sales Analysis ===")
print(f"Period: {dates[0].strftime('%b %Y')} to {dates[-1].strftime('%b %Y')}")
print(f"Total sales (4 years): £{ts.sum():.1f}M")
print(f"Average monthly sales: £{ts.mean():.2f}M")
print(f"Min: £{ts.min():.2f}M ({ts.idxmin().strftime('%b %Y')})")
print(f"Max: £{ts.max():.2f}M ({ts.idxmax().strftime('%b %Y')})")
# ADF test on raw series
_, p_raw, _, _, _, _ = adfuller(ts, autolag='AIC')
print(f"\nADF p-value (raw): {p_raw:.4f} → {'non-stationary' if p_raw >= 0.05 else 'stationary'}")
print("Log-transforming (multiplicative → additive)...")
ts_log = np.log(ts)
_, p_log, _, _, _, _ = adfuller(ts_log.diff().dropna(), autolag='AIC')
print(f"ADF p-value (log + diff): {p_log:.4f} → {'✅ stationary' if p_log < 0.05 else 'non-stationary'}")
# Multiplicative decomposition (correct for this data)
result = seasonal_decompose(ts, model='multiplicative', period=12,
extrapolate_trend='freq')
# Extract components
trend_est = result.trend
seasonal_est = result.seasonal
resid_est = result.resid
# Annual growth rate from trend
start_trend = trend_est.iloc[0]
end_trend = trend_est.iloc[-1]
total_growth = (end_trend / start_trend - 1) * 100
years_elapsed = n_months / 12
annual_growth_est = ((end_trend / start_trend) ** (1/years_elapsed) - 1) * 100
print(f"\n=== Trend Analysis ===")
print(f"Total growth (4 years): {total_growth:.1f}%")
print(f"Implied annual growth: {annual_growth_est:.1f}% (true: 3.2%)")
# Seasonal pattern
print(f"\n=== Seasonal Analysis ===")
seasonal_df = pd.DataFrame({
'month': dates.month,
'seasonal_factor': seasonal_est
}).groupby('month').mean().round(3)
seasonal_df.index = ['Jan','Feb','Mar','Apr','May','Jun',
'Jul','Aug','Sep','Oct','Nov','Dec']
print("Monthly seasonal factors (multiplicative, mean=1.0):")
for m, row in seasonal_df.iterrows():
factor = row['seasonal_factor']
bar = '█' * int(factor * 10)
change_pct = (factor - 1) * 100
print(f" {m:3s}: {factor:.3f} ({change_pct:+.0f}%) {bar}")
# Residual analysis
_, p_resid, _, _, _, _ = adfuller(resid_est.dropna(), autolag='AIC')
resid_mean = resid_est.mean()
print(f"\n=== Residual Diagnostics ===")
print(f"Residual mean: {resid_mean:.4f} (should be ~1.0 for multiplicative)")
print(f"ADF p-value: {p_resid:.4f} → {'✅ stationary residuals (good)' if p_resid < 0.05 else '❌ residuals still show pattern'}")
# ACF of residuals
acf_vals = acf(resid_est.dropna(), nlags=13, fft=True)
conf = 2 / np.sqrt(len(resid_est.dropna()))
print(f"\nACF of residuals (significant if |value| > {conf:.3f}):")
for lag in [1, 2, 3, 6, 12]:
flag = "⚠️ significant" if abs(acf_vals[lag]) > conf else "✅ not significant"
print(f" Lag {lag:2d}: {acf_vals[lag]:+.4f} {flag}")
print("\n=== Business Insights ===")
print(f"1. Year-over-year sales growth: ~{annual_growth_est:.1f}% annually")
dec_factor = seasonal_df.loc['Dec', 'seasonal_factor']
feb_factor = seasonal_df.loc['Feb', 'seasonal_factor']
print(f"2. December is {dec_factor:.1f}x average monthly sales")
print(f"3. February is the weakest month ({feb_factor:.2f}x average)")
print(f"4. Holiday season (Nov-Dec) = ~{(seasonal_df.loc['Nov','seasonal_factor']+dec_factor)*100/12*100:.0f}% of annual sales")
print(f"5. Residuals are stationary → decomposition captured all systematic patterns")These decomposition insights directly inform business decisions: staffing up 2.4× in December, negotiating supplier contracts based on the trend growth rate, scheduling store renovations in February (trough month), and using the trend + seasonal components as a baseline forecast. The stationary residuals confirm that the decomposition successfully extracted the systematic patterns — only genuine random variation remains. ACF of residuals showing no significant autocorrelation means no remaining predictable signal has been left on the table.
✍️ Practice Exercises
- Download the Air Passengers dataset (classic monthly airline passenger data, available as
seaborn.load_dataset('flights')). Apply both additive and multiplicative decomposition. Which model fits better — and how can you tell by looking at the residuals? - The ADF test gives different results on the original series vs the differenced series. Run ADF on: (a) the raw Air Passengers series, (b) after log transform, (c) after log + first diff, (d) after log + seasonal diff (lag=12). At which step does it first become stationary?
- Using the Air Passengers dataset, compute the ACF and PACF of the log-differenced series. Identify: (a) which lags are significant, (b) whether the pattern suggests AR, MA, or ARMA structure, (c) whether seasonal AR or MA terms are needed.
- Generate a time series with a structural break: 3 years of normal growth, then an abrupt step change (e.g., a new competitor enters the market). Run decomposition and ADF. Does the decomposition correctly capture the break? What happens to the residuals?
📚 Primary Source for This Lesson
Forecasting: Principles and Practice (3rd ed.) — Hyndman & Athanasopoulos
This free online textbook is the gold standard for learning time series analysis and forecasting. Chapters 1–3 cover everything in this lesson in greater depth: time series patterns, decomposition (STL vs classical), and transformations. The examples use R, but the concepts translate directly to Python. Highly recommended for Chapter 3 (time series decomposition) and Chapter 9 (ARIMA) as preparation for Lesson 36.