Time Series Analysis

1. Time Series Components

from statsmodels.tsa.seasonal import seasonal_decompose
import matplotlib.pyplot as plt

# model='additive' when variation is constant;
# model='multiplicative' when variance grows with the level
result = seasonal_decompose(df['spend'], model='additive', period=12)
result.plot()
plt.tight_layout()
plt.show()

2. Stationarity

A stationary series has constant mean, variance, and autocovariance over time. ARIMA requires stationarity because its AR and MA components model fixed relationships between past and present values — relationships that cannot be fixed if the series is trending.

Augmented Dickey-Fuller (ADF) Test

• H₀: Series has a unit root (non-stationary)
• Decision: Reject H₀ if p-value < 0.05 → series is stationary

from statsmodels.tsa.stattools import adfuller

result = adfuller(df['spend'])
print(f'ADF Statistic: {result[0]:.4f}')
print(f'p-value:       {result[1]:.4f}')
# p < 0.05 → reject unit root → stationary
# p ≥ 0.05 → fail to reject → non-stationary, needs differencing

Making a Series Stationary

import numpy as np

# First-order differencing — removes linear trend (d=1)
df['diff1'] = df['spend'].diff()

# Second-order differencing — removes quadratic trend (d=2)
df['diff2'] = df['spend'].diff().diff()

# Seasonal differencing — removes seasonal pattern (period s=12)
df['seasonal_diff'] = df['spend'].diff(12)

# Log transform first to stabilize variance, then difference
df['log_diff'] = np.log(df['spend']).diff()

# Always re-test with ADF after transforming
result = adfuller(df['diff1'].dropna())
print(f'After differencing — p-value: {result[1]:.4f}')

3. ACF and PACF

The ACF and PACF plots are used to identify the correct ARIMA order after making the series stationary.

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

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
plot_acf( df['diff1'].dropna(), ax=axes[0], lags=30, title='ACF')
plot_pacf(df['diff1'].dropna(), ax=axes[1], lags=30, title='PACF')
plt.tight_layout()
plt.show()
# Shaded region = 95% confidence interval
# Bars outside shading are statistically significant autocorrelations

4. ARIMA and SARIMA Models

ARIMA(p, d, q): $p$ = number of AR lags, $d$ = differencing order, $q$ = number of MA lags.

from statsmodels.tsa.arima.model import ARIMA

model  = ARIMA(df['spend'], order=(2, 1, 2))   # AR=2, d=1, MA=2
result = model.fit()
print(result.summary())

# Forecast 12 steps ahead
forecast = result.forecast(steps=12)
conf_int = result.get_forecast(steps=12).conf_int()

# Check residuals — should look like white noise
result.plot_diagnostics(figsize=(12, 8))
plt.show()

SARIMA adds seasonal terms: SARIMA(p,d,q)(P,D,Q)[s] where P, D, Q are the seasonal AR, differencing, and MA orders and s is the period (12 for monthly data).

from statsmodels.tsa.statespace.sarimax import SARIMAX

model = SARIMAX(df['spend'],
                order=(1, 1, 1),
                seasonal_order=(1, 1, 1, 12))
result = model.fit(disp=False)
print(result.summary())

Model Selection with AIC/BIC

best_aic = np.inf
best_order = None

for p in range(3):
    for q in range(3):
        try:
            m = ARIMA(df['spend'], order=(p, 1, q)).fit()
            if m.aic < best_aic:
                best_aic   = m.aic
                best_order = (p, 1, q)
        except Exception:
            continue

print(f"Best order: {best_order}  AIC: {best_aic:.2f}")

5. Train/Test Split & Cross-Validation

Never shuffle time series data. K-Fold cross-validation randomly reorders observations, which means a model can be trained on future data and tested on past data — leaking future information and producing overly optimistic estimates.

from sklearn.model_selection import TimeSeriesSplit
import numpy as np

# Simple hold-out: last N periods as test
n_test = 12
train = df.iloc[:-n_test]
test  = df.iloc[-n_test:]

# Walk-forward cross-validation (forward-chaining)
tscv = TimeSeriesSplit(n_splits=5)
for fold, (train_idx, test_idx) in enumerate(tscv.split(df)):
    train_fold = df.iloc[train_idx]
    test_fold  = df.iloc[test_idx]
    # train_fold always ends BEFORE test_fold begins — no leakage
    print(f"Fold {fold+1}: train={len(train_fold)}, test={len(test_fold)}")

6. Rolling Features for ML-Based Forecasting

For gradient boosting or neural network forecasters, transform the time series into a supervised learning problem using lag and rolling features. Always sort by date first.

import pandas as pd

df = df.sort_values('date').reset_index(drop=True)

# Lag features — capture recent history
df['lag_1']  = df['spend'].shift(1)    # yesterday
df['lag_7']  = df['spend'].shift(7)    # same day last week
df['lag_30'] = df['spend'].shift(30)   # same day last month
df['lag_365'] = df['spend'].shift(365) # year-over-year (if data available)

# Rolling aggregations — capture trend and variance
df['roll_mean_7']  = df['spend'].rolling(7).mean()
df['roll_std_7']   = df['spend'].rolling(7).std()
df['roll_mean_30'] = df['spend'].rolling(30).mean()
df['roll_max_30']  = df['spend'].rolling(30).max()

# Exponentially weighted mean — more weight to recent observations
df['ewm_7'] = df['spend'].ewm(span=7).mean()

# Date features
df['day_of_week']  = df['date'].dt.dayofweek   # 0=Mon, 6=Sun
df['is_weekend']   = (df['day_of_week'] >= 5).astype(int)
df['is_month_end'] = df['date'].dt.is_month_end.astype(int)
df['month']        = df['date'].dt.month

# Drop NaN rows created by shift/rolling
df = df.dropna().reset_index(drop=True)

7. Forecasting Evaluation Metrics

import numpy as np

def mae(y_true, y_pred):
    return np.mean(np.abs(y_true - y_pred))

def rmse(y_true, y_pred):
    return np.sqrt(np.mean((y_true - y_pred) ** 2))

def mape(y_true, y_pred):
    # UNDEFINED when y_true contains 0 — add guard
    mask = y_true != 0
    return np.mean(np.abs((y_true[mask] - y_pred[mask]) / y_true[mask])) * 100

def wape(y_true, y_pred):
    """Weighted Absolute Percentage Error — safe for zeros."""
    return np.sum(np.abs(y_true - y_pred)) / np.sum(np.abs(y_true)) * 100

y_true = np.array([100, 150, 200, 175, 0, 125])
y_pred = np.array([ 95, 160, 195, 180, 5, 120])

print(f"MAE:  {mae(y_true, y_pred):.2f}")
print(f"RMSE: {rmse(y_true, y_pred):.2f}")
print(f"WAPE: {wape(y_true, y_pred):.2f}%")
# Note: MAPE would be undefined here because y_true contains 0

7. Exponential Smoothing and Prophet

Exponential Smoothing Family

Exponential smoothing assigns exponentially decreasing weights to past observations. Unlike ARIMA, it doesn't require stationarity and handles trends and seasonality through additive components.

Holt-Winters update equations (additive seasonality):

$$L_t = \alpha(y_t - S_{t-m}) + (1-\alpha)(L_{t-1} + T_{t-1})$$

$$T_t = \beta(L_t - L_{t-1}) + (1-\beta)T_{t-1}$$

$$S_t = \gamma(y_t - L_t) + (1-\gamma)S_{t-m}$$

$$\hat{y}_{t+h} = L_t + h \cdot T_t + S_{t-m+h}$$

where $m$ is the seasonal period, $L_t$ is level, $T_t$ is trend, $S_t$ is seasonal component. All smoothing parameters $\alpha, \beta, \gamma \in (0,1)$ are fit by minimizing SSE.

from statsmodels.tsa.holtwinters import ExponentialSmoothing
import pandas as pd

# Monthly data with yearly seasonality
model = ExponentialSmoothing(
    train,
    trend='add',         # 'add' or 'mul'
    seasonal='add',      # 'add' or 'mul' or None
    seasonal_periods=12  # m = 12 for monthly data
)
fit = model.fit(optimized=True)  # auto-optimize alpha, beta, gamma
forecast = fit.forecast(steps=12)

print(f"Alpha (level): {fit.params['smoothing_level']:.3f}")
print(f"Beta (trend):  {fit.params['smoothing_trend']:.3f}")
print(f"Gamma (season):{fit.params['smoothing_seasonal']:.3f}")

Additive vs. multiplicative seasonality: use additive when seasonal swings are roughly constant in magnitude regardless of level (e.g., fixed dollar seasonality). Use multiplicative when seasonal swings scale with the level (e.g., percentage-based seasonality like retail — the holiday spike is larger when baseline sales are higher).

Prophet

Prophet (Meta, 2017) decomposes a time series as $y(t) = g(t) + s(t) + h(t) + \epsilon_t$: trend + seasonality + holidays + noise. Designed for business time series with strong seasonal patterns, irregular holidays, and missing data. Robust to outliers; handles changepoints automatically.

from prophet import Prophet
import pandas as pd

# Prophet requires columns: 'ds' (datetime) and 'y' (value)
df = pd.DataFrame({'ds': date_index, 'y': values})
df_train = df[df['ds'] < '2024-01-01']

model = Prophet(
    seasonality_mode='multiplicative',  # or 'additive'
    yearly_seasonality=True,
    weekly_seasonality=True,
    daily_seasonality=False,
    changepoint_prior_scale=0.05   # regularization on trend changepoints
)

# Add custom seasonality (e.g., monthly)
model.add_seasonality(name='monthly', period=30.5, fourier_order=5)

# Add holiday effects
holidays = pd.DataFrame({'ds': ['2023-12-25', '2024-12-25'], 'holiday': 'christmas'})
model = Prophet(holidays=holidays)

model.fit(df_train)

future = model.make_future_dataframe(periods=90)  # 90-day forecast
forecast = model.predict(future)

# Key columns in forecast:
# yhat: point forecast, yhat_lower/yhat_upper: 80% credible interval
model.plot(forecast)
model.plot_components(forecast)  # shows trend, yearly, weekly components

When to use Prophet: multiple seasonality periods, irregular holidays/events, business stakeholders who need interpretable components, datasets with missing values or outliers. Limitations: weaker than ARIMA/LightGBM at capturing autocorrelation; not ideal for high-frequency (sub-daily) data; probabilistic intervals are credible intervals, not prediction intervals.

8. Time Series Anomaly Detection

Anomaly detection in time series is harder than in i.i.d. data because normal values are autocorrelated — what's unusual is context-dependent. An approach that works well on cross-sectional data may perform poorly on sequential data.

Statistical Methods

Rolling Z-Score

Flag points more than $k$ standard deviations from a rolling mean. Simple and interpretable.

import pandas as pd
import numpy as np

def rolling_zscore_anomalies(series, window=30, threshold=3.0):
    rolling_mean = series.rolling(window).mean()
    rolling_std  = series.rolling(window).std()
    z_scores = (series - rolling_mean) / rolling_std
    return z_scores.abs() > threshold

anomalies = rolling_zscore_anomalies(df['value'], window=30, threshold=3.0)

STL Decomposition + Residual Thresholding

Decompose the series into trend + seasonality + residual. Flag anomalies in the residual component (which should be white noise if the model is correct).

from statsmodels.tsa.seasonal import STL

stl = STL(df['value'], seasonal=13)  # seasonal=13 for monthly data
result = stl.fit()

residuals = result.resid
threshold = 3 * residuals.std()
anomalies = residuals.abs() > threshold

# Plot
result.plot()  # shows trend, seasonal, residual components

Isolation Forest on Time Series Features

Extract temporal features (lags, rolling statistics) and apply Isolation Forest. This converts a time series problem into a standard anomaly detection problem on a feature matrix.

from sklearn.ensemble import IsolationForest
import pandas as pd
import numpy as np

def make_time_features(series, lags=[1,2,3,7], windows=[7,30]):
    df = series.to_frame(name='value')
    for lag in lags:
        df[f'lag_{lag}'] = series.shift(lag)
    for w in windows:
        df[f'roll_mean_{w}'] = series.rolling(w).mean()
        df[f'roll_std_{w}']  = series.rolling(w).std()
    return df.dropna()

features = make_time_features(df['value'])

iso = IsolationForest(contamination=0.05, random_state=42)
labels = iso.fit_predict(features)  # -1 = anomaly, 1 = normal

anomaly_idx = features.index[labels == -1]

Choosing an Approach

Evaluation challenge: anomaly labels are rare and expensive. Use precision at top-k (rank anomalies by score, check top k) and time to detection (how many timesteps before the anomaly is flagged). If labels exist, use F1-score but account for temporal proximity — a flag 1 timestep off from the true anomaly may still be acceptable.