|
| 1 | +"""Forecasting and anomaly detection on country-level emissions. |
| 2 | +
|
| 3 | +Simple, dependency-free approach designed to run inside the Prefect flow |
| 4 | +*after* dbt has built ``mart.mart_country_emissions``: |
| 5 | +
|
| 6 | +* **Forecast:** ordinary least squares on (year, total_emissions) per country, |
| 7 | + projected ``horizon`` years forward, with a +/- 1.96 * residual-stderr band. |
| 8 | +* **Anomalies:** rolling z-score (window=3) on year-over-year change; rows with |
| 9 | + ``|z| > 2.5`` are flagged. |
| 10 | +
|
| 11 | +Both outputs are written to dedicated mart tables so the dashboard can read |
| 12 | +them with the same cached ``query()`` helper as everything else. |
| 13 | +""" |
| 14 | + |
| 15 | +from __future__ import annotations |
| 16 | + |
| 17 | +from dataclasses import dataclass |
| 18 | + |
| 19 | +import numpy as np |
| 20 | +import pandas as pd |
| 21 | +import structlog |
| 22 | + |
| 23 | +log = structlog.get_logger(__name__) |
| 24 | + |
| 25 | + |
| 26 | +@dataclass(frozen=True) |
| 27 | +class ForecastConfig: |
| 28 | + """Configuration for the forecast task.""" |
| 29 | + |
| 30 | + horizon_years: int = 5 |
| 31 | + min_history_years: int = 4 |
| 32 | + anomaly_z_threshold: float = 2.5 |
| 33 | + anomaly_window: int = 3 |
| 34 | + |
| 35 | + |
| 36 | +def _fit_linear(years: np.ndarray, values: np.ndarray) -> tuple[float, float, float]: |
| 37 | + """Return (slope, intercept, residual_stderr) from a 1-D OLS fit.""" |
| 38 | + n = len(years) |
| 39 | + if n < 2: |
| 40 | + return 0.0, float(values.mean() if n else 0.0), 0.0 |
| 41 | + x_mean = years.mean() |
| 42 | + y_mean = values.mean() |
| 43 | + denom = ((years - x_mean) ** 2).sum() |
| 44 | + if denom == 0: |
| 45 | + return 0.0, float(y_mean), 0.0 |
| 46 | + slope = float(((years - x_mean) * (values - y_mean)).sum() / denom) |
| 47 | + intercept = float(y_mean - slope * x_mean) |
| 48 | + fitted = slope * years + intercept |
| 49 | + residuals = values - fitted |
| 50 | + dof = max(n - 2, 1) |
| 51 | + stderr = float(np.sqrt((residuals**2).sum() / dof)) |
| 52 | + return slope, intercept, stderr |
| 53 | + |
| 54 | + |
| 55 | +def forecast_country_emissions( |
| 56 | + history: pd.DataFrame, config: ForecastConfig | None = None |
| 57 | +) -> pd.DataFrame: |
| 58 | + """Project emissions ``horizon_years`` forward for every country. |
| 59 | +
|
| 60 | + Args: |
| 61 | + history: Output of ``mart.mart_country_emissions`` with columns |
| 62 | + ``country_code``, ``year``, ``total_emissions_tonnes``. |
| 63 | + config: Forecast configuration; defaults to ``ForecastConfig()``. |
| 64 | +
|
| 65 | + Returns: |
| 66 | + DataFrame with columns ``country_code``, ``year``, ``forecast_tonnes``, |
| 67 | + ``lower_band``, ``upper_band``, ``model``. |
| 68 | + """ |
| 69 | + cfg = config or ForecastConfig() |
| 70 | + out_rows: list[dict[str, object]] = [] |
| 71 | + |
| 72 | + for country, group in history.groupby("country_code", sort=True): |
| 73 | + g = group.sort_values("year") |
| 74 | + if len(g) < cfg.min_history_years: |
| 75 | + log.warning("forecast_skipped_insufficient_history", country=country, years=len(g)) |
| 76 | + continue |
| 77 | + years = g["year"].to_numpy(dtype=float) |
| 78 | + values = g["total_emissions_tonnes"].to_numpy(dtype=float) |
| 79 | + slope, intercept, stderr = _fit_linear(years, values) |
| 80 | + |
| 81 | + last_year = int(years.max()) |
| 82 | + for h in range(1, cfg.horizon_years + 1): |
| 83 | + yr = last_year + h |
| 84 | + point = slope * yr + intercept |
| 85 | + band = 1.96 * stderr * np.sqrt(h) # widening band |
| 86 | + out_rows.append( |
| 87 | + { |
| 88 | + "country_code": country, |
| 89 | + "year": yr, |
| 90 | + "forecast_tonnes": max(0.0, point), |
| 91 | + "lower_band": max(0.0, point - band), |
| 92 | + "upper_band": max(0.0, point + band), |
| 93 | + "model": "ols_linear", |
| 94 | + } |
| 95 | + ) |
| 96 | + |
| 97 | + return pd.DataFrame(out_rows) |
| 98 | + |
| 99 | + |
| 100 | +def detect_anomalies(history: pd.DataFrame, config: ForecastConfig | None = None) -> pd.DataFrame: |
| 101 | + """Flag country-years whose YoY change is a rolling z-score outlier. |
| 102 | +
|
| 103 | + Args: |
| 104 | + history: ``mart.mart_country_emissions`` rows. |
| 105 | + config: Anomaly window + threshold. |
| 106 | +
|
| 107 | + Returns: |
| 108 | + DataFrame containing only the anomalous rows, with the computed |
| 109 | + ``yoy_pct``, ``z_score``, and a string ``severity`` label. |
| 110 | + """ |
| 111 | + cfg = config or ForecastConfig() |
| 112 | + df = history.sort_values(["country_code", "year"]).copy() |
| 113 | + df["yoy_pct"] = df.groupby("country_code")["total_emissions_tonnes"].pct_change() * 100.0 |
| 114 | + |
| 115 | + # Compute the rolling baseline from prior rows only (shift by 1) so that |
| 116 | + # an anomalous point does not contaminate its own z-score. |
| 117 | + grouped = df.groupby("country_code")["yoy_pct"] |
| 118 | + prior = grouped.shift(1) |
| 119 | + rolling_mean = prior.groupby(df["country_code"]).transform( |
| 120 | + lambda s: s.rolling(cfg.anomaly_window, min_periods=2).mean() |
| 121 | + ) |
| 122 | + rolling_std = prior.groupby(df["country_code"]).transform( |
| 123 | + lambda s: s.rolling(cfg.anomaly_window, min_periods=2).std() |
| 124 | + ) |
| 125 | + df["z_score"] = (df["yoy_pct"] - rolling_mean) / rolling_std.replace(0, np.nan) |
| 126 | + |
| 127 | + anomalies = df[df["z_score"].abs() > cfg.anomaly_z_threshold].copy() |
| 128 | + anomalies["severity"] = np.where(anomalies["z_score"].abs() > 4.0, "critical", "warning") |
| 129 | + return anomalies[ |
| 130 | + ["country_code", "year", "total_emissions_tonnes", "yoy_pct", "z_score", "severity"] |
| 131 | + ].reset_index(drop=True) |
0 commit comments