Time Series Analysis
Overview
Time series analysis covers the end-to-end workflow for sequential temporal data: data preparation, stationarity testing, decomposition, autocorrelation analysis, model selection, temporally-correct cross-validation, and evaluation. Apply this skill for any dataset where observations are ordered in time and the temporal structure is scientifically meaningful.
When to Use This Skill
Use this skill when:
- The dataset has a datetime index and the ordering of observations matters
- You need to forecast future values of a variable
- You need to detect trends, seasonality, or change points
- You are performing anomaly detection on temporal data
- You need to evaluate a forecasting model with proper temporal splits
- You are fitting ARIMA, Prophet, LSTM, or any other time series model
Data Preparation
Datetime Index Setup
python
import pandas as pd
# Load and parse dates
df = pd.read_csv("data.csv", parse_dates=["timestamp"])
df = df.sort_values("timestamp").reset_index(drop=True)
df = df.set_index("timestamp")
# Verify index is datetime
print(df.index.dtype) # should be datetime64[ns]
print(df.index.is_monotonic_increasing) # must be True
Frequency Detection and Setting
python
# Infer frequency from data
inferred_freq = pd.infer_freq(df.index)
print(f"Inferred frequency: {inferred_freq}")
# Common codes: T (minute), H (hourly), D (daily), W (weekly), M (month-end), Q (quarter), A (annual)
# Explicitly set frequency (required by many statsmodels functions)
df = df.asfreq(inferred_freq) # may introduce NaT rows for missing timestamps
print(df.index.freq)
Handling Missing Timestamps
python
# Check for gaps
full_index = pd.date_range(start=df.index.min(), end=df.index.max(), freq=inferred_freq)
missing_timestamps = full_index.difference(df.index)
print(f"Missing timestamps: {len(missing_timestamps)}")
# Reindex to fill in gaps, then interpolate
df = df.reindex(full_index)
n_missing = df["value"].isna().sum()
print(f"Missing values after reindex: {n_missing}")
# Interpolation options (choose based on expected pattern)
df["value"] = df["value"].interpolate(method="time") # linear in time
# df["value"] = df["value"].interpolate(method="cubic") # cubic spline
# df["value"] = df["value"].ffill() # forward fill (flat)
# df["value"] = df["value"].fillna(df["value"].rolling(7, min_periods=1).mean()) # rolling mean
Resampling
python
# Downsample: daily → weekly mean
df_weekly = df.resample("W").mean()
# Downsample: minute-level → hourly sum
df_hourly = df.resample("H").sum()
# Upsample: daily → hourly (creates NaN, then fill)
df_hourly = df.resample("H").interpolate(method="time")
# Multiple aggregations at once
df_monthly = df.resample("M").agg({
"sales": "sum",
"temperature": "mean",
"visitors": "max"
})
Stationarity Analysis
Why it matters: Most classical models (ARIMA) require stationarity. Non-stationary series have time-varying mean or variance, which violates model assumptions.
ADF Test (Augmented Dickey-Fuller)
Null hypothesis: the series has a unit root (non-stationary). Reject H0 for stationarity.
python
from statsmodels.tsa.stattools import adfuller
def test_stationarity_adf(series, name="Series", alpha=0.05):
"""Run ADF test and print interpretation."""
result = adfuller(series.dropna(), autolag="AIC")
adf_stat, p_value, n_lags, n_obs, critical_values, _ = result
print(f"\nADF Test — {name}")
print(f" ADF Statistic: {adf_stat:.4f}")
print(f" p-value: {p_value:.4f}")
print(f" Lags used: {n_lags}")
print(f" Observations: {n_obs}")
print(f" Critical values:")
for level, crit in critical_values.items():
marker = " <-- ADF statistic is more extreme" if adf_stat < crit else ""
print(f" {level}: {crit:.4f}{marker}")
if p_value < alpha:
print(f" CONCLUSION: Reject H0 (p={p_value:.4f} < {alpha}). Series is STATIONARY.")
else:
print(f" CONCLUSION: Fail to reject H0 (p={p_value:.4f} >= {alpha}). Series is NON-STATIONARY.")
return {"adf_stat": adf_stat, "p_value": p_value, "stationary": p_value < alpha}
result = test_stationarity_adf(df["value"], name="Sales")
KPSS Test
Null hypothesis: the series is stationary (trend-stationary). Reject H0 for non-stationarity. Use alongside ADF — contradictory results suggest difference-stationary behavior.
python
from statsmodels.tsa.stattools import kpss
def test_stationarity_kpss(series, name="Series", regression="c", alpha=0.05):
"""
regression: 'c' (constant) or 'ct' (constant + trend)
"""
kpss_stat, p_value, n_lags, critical_values = kpss(series.dropna(), regression=regression, nlags="auto")
print(f"\nKPSS Test — {name}")
print(f" KPSS Statistic: {kpss_stat:.4f}")
print(f" p-value: {p_value:.4f} (Note: p is bounded [0.01, 0.10])")
print(f" Lags used: {n_lags}")
print(f" Critical values: {critical_values}")
if p_value < alpha:
print(f" CONCLUSION: Reject H0. Series is NON-STATIONARY (unit root present).")
else:
print(f" CONCLUSION: Fail to reject H0. Series is STATIONARY.")
return {"kpss_stat": kpss_stat, "p_value": p_value, "stationary": p_value >= alpha}
# Combining ADF and KPSS
# Both say stationary → Strong evidence of stationarity
# Both say non-stationary → Strong evidence; needs differencing
# ADF: stationary, KPSS: non-stationary → Trend-stationary (remove trend)
# ADF: non-stationary, KPSS: stationary → Difference-stationary (difference it)
Making a Non-Stationary Series Stationary
python
import numpy as np
import matplotlib.pyplot as plt
series = df["value"]
# 1. First-order differencing (removes linear trend)
diff1 = series.diff().dropna()
# 2. Log transform (stabilizes variance for exponential growth)
log_series = np.log(series.replace(0, np.nan))
# 3. Log + differencing (common for financial/economic data)
log_diff = log_series.diff().dropna()
# 4. Seasonal differencing (removes annual seasonality in monthly data)
# For monthly data with annual seasonality:
seasonal_diff = series.diff(12).dropna()
# 5. Double differencing (second-order non-stationarity)
diff2 = series.diff().diff().dropna()
# Test each until stationary — prefer the fewest differences needed
for name, s in [("Original", series), ("Diff(1)", diff1), ("Log", log_series), ("Log+Diff", log_diff)]:
test_stationarity_adf(s, name=name)
# Plot to visually confirm
fig, axes = plt.subplots(4, 1, figsize=(12, 12))
for ax, (name, s) in zip(axes, [("Original", series), ("Diff(1)", diff1),
("Log", log_series), ("Log+Diff", log_diff)]):
ax.plot(s)
ax.set_title(name)
ax.set_xlabel("Time")
plt.tight_layout()
plt.savefig("stationarity_transforms.png", dpi=150)
Decomposition
When to Use Additive vs Multiplicative
- Additive (trend + seasonal + residual): seasonal fluctuations are constant in magnitude. Use when the series is roughly linear.
- Multiplicative (trend × seasonal × residual): seasonal fluctuations grow with the level. Use when the series grows exponentially. Can be applied as
log → additive decomposition → exp.
Classical Decomposition
python
from statsmodels.tsa.seasonal import seasonal_decompose
import matplotlib.pyplot as plt
# period: number of observations per seasonal cycle
# e.g., 12 for monthly data with annual seasonality, 7 for daily with weekly seasonality
result = seasonal_decompose(df["value"], model="additive", period=12)
fig, axes = plt.subplots(4, 1, figsize=(12, 10))
result.observed.plot(ax=axes[0], title="Observed")
result.trend.plot(ax=axes[1], title="Trend")
result.seasonal.plot(ax=axes[2], title="Seasonal")
result.resid.plot(ax=axes[3], title="Residual")
plt.tight_layout()
plt.savefig("decomposition_classical.png", dpi=150)
plt.show()
# Residual diagnostics
resid = result.resid.dropna()
print(f"Residual mean: {resid.mean():.4f} (should be ~0)")
print(f"Residual std: {resid.std():.4f}")
STL Decomposition (Robust, Recommended)
STL (Seasonal-Trend decomposition using LOESS) is more robust to outliers and handles non-constant seasonality.
python
from statsmodels.tsa.seasonal import STL
stl = STL(
df["value"],
period=12, # seasonal period
seasonal=13, # seasonal smoother window (odd, >= 7)
trend=None, # auto-selected based on period
robust=True # downweight outliers in LOESS
)
result = stl.fit()
fig = result.plot()
fig.set_size_inches(12, 10)
plt.tight_layout()
plt.savefig("decomposition_stl.png", dpi=150)
# Extract components
trend = result.trend
seasonal = result.seasonal
residual = result.resid
# Seasonal strength: how much variance is explained by the seasonal component
var_resid = residual.var()
var_seasonal_plus_resid = (seasonal + residual).var()
seasonal_strength = max(0, 1 - var_resid / var_seasonal_plus_resid)
print(f"Seasonal strength: {seasonal_strength:.3f} (>0.6 = strong seasonality)")
Detecting Seasonality: ACF and FFT
python
from statsmodels.graphics.tsaplots import plot_acf
import numpy as np
import matplotlib.pyplot as plt
# ACF approach: peaks at multiples of the seasonal period
fig, ax = plt.subplots(figsize=(12, 4))
plot_acf(df["value"].dropna(), lags=60, ax=ax)
ax.set_title("ACF — look for peaks at lag s, 2s, 3s...")
plt.savefig("acf_seasonality.png", dpi=150)
# FFT periodogram approach
series_values = df["value"].dropna().values
n = len(series_values)
fft_vals = np.abs(np.fft.rfft(series_values - series_values.mean())) ** 2
freqs = np.fft.rfftfreq(n)
# Find dominant frequency
dominant_freq_idx = np.argmax(fft_vals[1:]) + 1 # skip DC component
dominant_period = 1 / freqs[dominant_freq_idx] if freqs[dominant_freq_idx] > 0 else np.inf
print(f"Dominant period from FFT: {dominant_period:.1f} time steps")
plt.figure(figsize=(12, 4))
plt.plot(1/freqs[1:], fft_vals[1:])
plt.xlabel("Period (time steps)")
plt.ylabel("Power")
plt.title("FFT Periodogram")
plt.xlim(0, n // 2)
plt.savefig("fft_periodogram.png", dpi=150)
Autocorrelation Analysis
Autocorrelation analysis reveals the temporal memory structure of the series and guides model order selection.
python
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.stats.diagnostic import acorr_ljungbox
import matplotlib.pyplot as plt
# ACF: lag at which autocorrelation cuts off → MA(q) order
# PACF: lag at which partial autocorrelation cuts off → AR(p) order
fig, axes = plt.subplots(2, 1, figsize=(12, 8))
plot_acf(df["value"].dropna(), lags=40, ax=axes[0], alpha=0.05)
axes[0].set_title("ACF — indicates MA order (q): note where it drops inside confidence band")
plot_pacf(df["value"].dropna(), lags=40, ax=axes[1], alpha=0.05, method="ywmle")
axes[1].set_title("PACF — indicates AR order (p): note where it drops inside confidence band")
plt.tight_layout()
plt.savefig("acf_pacf.png", dpi=150)
plt.show()
# ACF / PACF interpretation guide:
# AR(p): ACF decays exponentially/sinusoidally; PACF cuts off at lag p
# MA(q): PACF decays; ACF cuts off at lag q
# ARMA: Both decay gradually → use AIC/BIC to select (p, q)
# Ljung-Box test: are residuals white noise?
# Apply to residuals AFTER fitting a model
lb_result = acorr_ljungbox(df["value"].dropna(), lags=[10, 20, 30], return_df=True)
print(lb_result)
# If p-value > 0.05 at all lags → residuals are white noise → good model fit
# If p-value < 0.05 → model has not captured all autocorrelation structure
Forecasting Method Selection