The model, formula by formula

def corrected_electricity(
    daily: pd.Series, temp: pd.Series, fit_until: pd.Timestamp, window_years: int | None = None, holidays: bool = False
) -> pd.Series:
    """Monthly mean of log consumption after removing temperature, weekday and calendar-month effects (D4).

    The coefficients are fitted only on days before `fit_until`, so a backtest never sees the future.
    The two options are the variants tested in D24 (not used by the published model):
    window_years: each calendar year's residuals come from a fit on the `window_years` years ending with it.
    holidays: add a French public-holiday dummy.
    """
    df = pd.DataFrame({"y": np.log(daily)}).join(temp.rename("t"), how="inner")
    idx = df.index
    X = pd.concat([
        pd.DataFrame({
            "const": 1.0,
            "hdd": (HDD_BASE - df["t"]).clip(lower=0),
            "cdd": (df["t"] - CDD_BASE).clip(lower=0),
        }, index=idx),
        pd.get_dummies(idx.dayofweek, prefix="dow", drop_first=True, dtype=float).set_index(idx),
        pd.get_dummies(idx.month, prefix="m", drop_first=True, dtype=float).set_index(idx),
        *([pd.DataFrame({"holiday": idx.isin(french_holidays(set(idx.year))).astype(float)}, index=idx)] if holidays else []),
    ], axis=1)
    fit = idx < fit_until
    if window_years is None:
        beta, *_ = np.linalg.lstsq(X[fit].to_numpy(), df["y"][fit].to_numpy(), rcond=None)
        resid = df["y"] - X.to_numpy() @ beta
    else:
        resid = pd.Series(np.nan, index=idx)
        for y in sorted(set(idx.year)):
            window = fit & (idx.year > y - window_years) & (idx.year <= y)
            if window.sum() < 200:  # too few days in the first years: use everything known so far
                window = fit & (idx.year <= y)
            beta, *_ = np.linalg.lstsq(X[window].to_numpy(), df["y"][window].to_numpy(), rcond=None)
            year = idx.year == y
            resid[year] = df["y"][year] - X[year].to_numpy() @ beta
    month = resid.groupby(idx.to_period("M"))
    out = month.mean()[month.count() >= MIN_DAYS]
    out.attrs["weather"] = {"hdd": float(beta[1]), "cdd": float(beta[2]), "days": int(fit.sum())}  # shown on the model page
    return out

RIDGE_ALPHAS = np.logspace(-2, 3, 30)  # ridge penalties tried; RidgeCV picks one by leave-one-out error
MODELS = {
    "ar": (["y_lag2"], LinearRegression),
    "ridge": (FEATURES, lambda: make_pipeline(StandardScaler(), RidgeCV(alphas=RIDGE_ALPHAS))),
    "gbm": (FEATURES, lambda: HistGradientBoostingRegressor(max_depth=3, learning_rate=0.05, max_iter=200, random_state=0)),
}

def training_rows(f: pd.DataFrame, t: pd.Period) -> pd.DataFrame:
    """Months whose IPI was already published when nowcasting t (up to t-2), with complete features.

    COVID months are left out of training (still scored in the backtest): a 25% collapse and
    rebound teaches the models relationships that do not hold in normal months.
    """
    return f.loc[: t - 2].drop(COVID, errors="ignore").dropna(subset=["y", *FEATURES])

def explain(raw) -> dict:
    """The published ridge model taken apart: coefficients, and each pending month's calculation step by step.

    prediction = intercept + sum_j coef_j * (x_j - mean_j) / scale_j. The assert proves the
    decomposition shown on the website is exactly what the model computes.
    """
    f, train, pending = final_sets(raw)
    model = MODELS["ridge"][1]().fit(train[FEATURES], train["y"])
    scaler, ridge = model[0], model[-1]
    z = (pending[FEATURES] - scaler.mean_) / scaler.scale_
    contrib = z * ridge.coef_
    pred = ridge.intercept_ + contrib.sum(axis=1)
    assert np.allclose(pred, model.predict(pending[FEATURES])), "decomposition differs from the model"
    ar = LinearRegression().fit(train[["y_lag2"]], train["y"])
    return {
        "alpha": float(ridge.alpha_),
        "alphas": [float(RIDGE_ALPHAS[0]), float(RIDGE_ALPHAS[-1]), len(RIDGE_ALPHAS)],
        "intercept": float(ridge.intercept_),
        "train": {"n": len(train), "from": str(train.index.min()), "to": str(train.index.max())},
        "features": [
            {"name": c, "mean": float(m), "scale": float(s), "coef": float(w)}
            for c, m, s, w in zip(FEATURES, scaler.mean_, scaler.scale_, ridge.coef_)
        ],
        "months": [
            {"m": str(t), "x": pending.loc[t, FEATURES].astype(float).tolist(), "z": z.loc[t].tolist(),
             "contrib": contrib.loc[t].tolist(), "pred": float(pred[t])}
            for t in pending.index
        ],
        "ar": {"a": float(ar.intercept_), "b": float(ar.coef_[0])},
        "weather": f.attrs["weather"],
    }