Time Series Analysis skill for Codex

Sound time series work starts from a regular, time-zone-correct index with gaps made explicit, separates trend, seasonality and noise before interpreting anything, compares like periods with like, and judges any forecast against a seasonal-naive baseline on data it did not see.

Inspect first

1. Read AGENTS.md for the stack (pandas 2.x, Polars, SQL warehouse) and where the data lives. Load a sample with the timestamp column parsed.
2. Establish the basics and write them down:

import pandas as pd

s = df.set_index("ts")["value"].sort_index()
print(s.index.tz, s.index.min(), s.index.max(), len(s))
print(s.index.to_series().diff().value_counts().head())   # actual spacing
print(s.index.duplicated().sum(), "duplicate timestamps")

3. Decide: is the value a level (temperature, balance), a flow over an interval (orders per hour), or an event stream (individual transactions)? It decides how to aggregate.
4. Ask which time zone the business thinks in; daily and weekly cuts must use it.

Make the index regular

Value type	Resample with	Missing intervals mean
Event stream (each row an event)	.resample("h").size() or .sum() of amounts	Zero, filled by resample automatically
Flow already aggregated	.resample("D").sum(min_count=1)	Unknown: keep NaN unless the source guarantees zero
Level / gauge	.resample("h").mean() or .last()	Unknown: interpolate only short gaps

s = s.tz_convert("Europe/Berlin")                 # if stored as UTC-aware
daily = s.resample("D").sum(min_count=1)          # NaN, not 0, for days with no data
gaps = daily[daily.isna()]
level = gauge.resample("15min").mean().interpolate(limit=4, limit_area="inside")

• pandas 2.2+ aliases: "h", "min", "s", "ME" (month end), "MS" (month start), "QE", "YE". The old "H", "T", "M" raise deprecation warnings.
• Naive timestamps: tz_localize("Europe/Berlin", ambiguous="infer", nonexistent="shift_forward") handles DST transitions; check the 23- and 25-hour days afterwards.
• asfreq selects existing points at a frequency; resample aggregates. Do not use asfreq on event data.

Look before modelling

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

fig, axes = plt.subplots(3, 1, figsize=(10, 8), layout="constrained")
daily.plot(ax=axes[0], title="Daily")
daily.rolling(28, min_periods=21).mean().plot(ax=axes[0])
daily.groupby(daily.index.dayofweek).mean().plot.bar(ax=axes[1], title="Mean by weekday")
daily.groupby(daily.index.month).mean().plot.bar(ax=axes[2], title="Mean by month")
fig.savefig("figures/ts_overview.png")

Check for level shifts (tracking changes, product launches), outliers, changing variance (log transform if variance grows with level), and the seasonal periods present (weekly for daily data; daily and weekly for hourly data).

Decomposition

import numpy as np
from statsmodels.tsa.seasonal import STL

y = np.log1p(daily.ffill(limit=2).dropna())
res = STL(y, period=7).fit()   # with heavy outliers, enable its outlier-resistant fitting option
trend, seasonal, resid = res.trend, res.seasonal, res.resid
strength_seasonal = max(0, 1 - resid.var() / (seasonal + resid).var())

• STL requires no NaN; fill short gaps or drop the leading section.
• Multiple seasonalities (hourly data with daily and weekly cycles): statsmodels.tsa.seasonal.MSTL with periods=(24, 168).
• Seasonal strength above about 0.6 means comparisons must be seasonally aligned.

Comparing periods

• Compare the same weekday mix: week over week, or same period last year aligned by weekday (364 days back), not by date.
• Use trailing 7- or 28-day windows to remove weekday effects in trends.
• Calendar effects: holidays, month length, Easter moving between March and April, paydays. Flag them before attributing changes to anything else.

Anomaly detection

Start simple and explainable:

def anomalies(y: pd.Series, window: int = 28, k: float = 4.0) -> pd.DataFrame:
    med = y.rolling(window, min_periods=window // 2, center=False).median().shift(1)
    mad = (y - med).abs().rolling(window, min_periods=window // 2).median().shift(1)
    score = 0.6745 * (y - med) / mad.replace(0, np.nan)
    return pd.DataFrame({"y": y, "expected": med, "score": score, "flag": score.abs() > k})

• Run it on STL residuals when the series is seasonal, otherwise every Monday is an anomaly.
• shift(1) keeps the window strictly in the past, so the method works in production where the future is unknown.
• Tune k against labelled past incidents if any exist; report precision on them.
• Change points (persistent level shifts) are different from spikes: compare rolling medians before and after, or use a change-point library (ruptures) if installed.

Baseline forecasts and backtesting

Always compute baselines before anything sophisticated:

from statsmodels.tsa.holtwinters import ExponentialSmoothing

def backtest(y: pd.Series, horizon: int = 28, folds: int = 4, season: int = 7) -> pd.DataFrame:
    rows = []
    for k in range(folds, 0, -1):
        cut = len(y) - k * horizon
        train, test = y.iloc[:cut], y.iloc[cut:cut + horizon]
        snaive = np.resize(train.iloc[-season:].to_numpy(), len(test))
        ets = ExponentialSmoothing(train, trend="add", seasonal="add",
                                   seasonal_periods=season).fit().forecast(len(test))
        scale = np.mean(np.abs(np.diff(train.to_numpy(), n=season)))   # MASE denominator
        for name, pred in {"seasonal_naive": snaive, "ets": ets.to_numpy()}.items():
            rows.append({"fold": cut, "model": name,
                         "mae": np.mean(np.abs(test.to_numpy() - pred)),
                         "mase": np.mean(np.abs(test.to_numpy() - pred)) / scale})
    return pd.DataFrame(rows).groupby("model")[["mae", "mase"]].mean()

• Never shuffle or use random K-fold on time series; train only on the past.
• MASE below 1 means better than the in-sample seasonal naive; a model that does not beat seasonal naive on the backtest is not worth shipping.
• Report forecasts with intervals, not single lines.
• For production forecasting, feature-rich models or hierarchical forecasts, hand off to $time-series-forecasting.

Gotchas

• Summing hourly data in UTC to “days” for a business in UTC-8 shifts a third of each day; convert time zones before resampling.
• The last, partial period (today, this month) looks like a drop. Exclude or label it.
• Interpolating long gaps invents data; limit interpolation and report gap counts.
• Rolling means on data with NaN: set min_periods, or the window silently shrinks.
• Autocorrelation makes naive t-tests on before/after means far too confident; compare against the seasonal pattern or use the backtest error as the noise level.
• Ratios (conversion) must be resampled as sum of numerator / sum of denominator.

Verification

• Index is monotonic, unique and at the declared frequency (pd.infer_freq returns it).
• Resampled totals equal raw totals for the same window (for flows and events).
• Decomposition components add back to the series (trend + seasonal + resid).
• Backtest table shows each model’s error versus seasonal naive over at least 3 folds.

Report template

**Series:** daily orders, 2023-01-01 to 2025-09-30, Europe/Berlin days; 6 missing days (pipeline outages) left NaN
**Seasonality:** strong weekly (strength 0.81), mild yearly peak in Nov-Dec
**Trend:** +14% year over year on 28-day trailing average; level shift +9% on 2025-03-04 (pricing change)
**Anomalies:** 5 flagged (|score| > 4 on STL residuals); 3 match known incidents, 2 unexplained (listed)
**Baseline forecast:** next 28 days, ETS MASE 0.78 vs seasonal naive 1.00 over 4 folds; 80% intervals in figures/forecast.png
**Caveats:** Black Friday not modelled explicitly; forecast assumes no further price changes