From 9a845646fde9bfa380afc011b032cbc6c0ef992e Mon Sep 17 00:00:00 2001 From: Light1YT Date: Thu, 4 Jun 2026 11:39:32 +0500 Subject: [PATCH] fix(#978): train-only detrend in rate backtest + Almon distributed-lag regression MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit REOPENED 951-B §9.6. PART A: fix look-ahead leakage in backtest_rate_sensitivity --detrend. The ln(units) trend was fit over train+test then split, so test data shaped the detrend and inflated the OOS hit-rate. _detrend_log now takes fit_n; backtest_tier fits the trend on TRAIN months only (same split evaluate_oos uses) and projects (a,b) point-in-time onto test. Default fit_n=None preserves prior behaviour. PART B (DoD): new app/services/forecasting/regression.py — Almon polynomial distributed-lag (deg 2) of Δln(district demand) on Δkey_rate lags 0..6 via OLS-on-Almon-regressors (numpy lstsq) + per-lag reconstruction + manual Newey-West HAC SEs (NO statsmodels). Output {best_lag_months, coef=long-run multiplier, x_pct, r2, n, per_lag_coef, hac_se,...}; gate mirrors _elasticity_coef (n<30 OR R²<0.1 OR Σβ≥0 → fallback); §9.6 phrase from the lag shape. ADVISORY, shipped standalone (integration point documented), NOT wired — protects the live compute_rate_sensitivity consumers. 125+31 tests (synthetic known-lag recovery, HAC computed/differs-from-OLS, fallback gating, no-leakage detrend). ruff clean. Refs #978 --- .../app/services/forecasting/regression.py | 629 ++++++++++++++++++ backend/scripts/backtest_rate_sensitivity.py | 64 +- .../scripts/test_backtest_rate_sensitivity.py | 109 +++ .../services/forecasting/test_regression.py | 513 ++++++++++++++ 4 files changed, 1302 insertions(+), 13 deletions(-) create mode 100644 backend/app/services/forecasting/regression.py create mode 100644 backend/tests/services/forecasting/test_regression.py diff --git a/backend/app/services/forecasting/regression.py b/backend/app/services/forecasting/regression.py new file mode 100644 index 00000000..c048c837 --- /dev/null +++ b/backend/app/services/forecasting/regression.py @@ -0,0 +1,629 @@ +"""§9.6 distributed-lag регрессия спроса района на ключевую ставку (Almon / ADL). + +Forgejo #978 (#951-B), §9.6 «lagged key_rate → demand». DoD требует *настоящую* +distributed-lag модель отклика месячного спроса района на key_rate при лагах 0..6, +а НЕ текущий single-lag OLS (`rate_sensitivity.best_lag` берёт ОДИН лучший лаг) и +НЕ unconstrained free-lags (7 коллинеарных лагов Δrate → раздутые дисперсии, +скачущие знаки — оценка непригодна на коротком месячном ряде). + +ПОДХОД — Almon polynomial distributed lag (Almon 1965; Stock & Watson, DL): +накладываем НИЗКОСТЕПЕННОЙ полином на 7 лаговых коэффициентов β_0..β_6, + β_j = Σ_{p=0..d} γ_p · j^p (d = _ALMON_DEGREE, по умолчанию 2), +и оцениваем d+1 параметр γ (а не 7 свободных β) через OLS на Almon-преобразованных +регрессорах z_p[t] = Σ_j j^p · x[t−j]. Это резко снижает коллинеарность (3 гладких +параметра вместо 7 шумных) — стандартный приём для коротких лаговых рядов. Per-lag +β_j реконструируем обратно из γ. Альтернатива ADL(p,q) задокументирована, но Almon +выбран: он напрямую даёт ФОРМУ отклика по лагам, нужную для фразы §9.6. + +HAC (Newey-West) стандартные ошибки — РУЧНОЙ numpy, БЕЗ statsmodels (тяжёлая +зависимость + пересборка prod-образа; ручной NW ~30 строк). Δln-остатки на +месячном ряде автокоррелированы (перекрытие лаговых окон + инерция спроса) → +обычные OLS-SE занижены. NW-ковариация = взвешенная сумма автоковариаций остатков +с окном Бартлетта; bandwidth L = floor(4·(n/100)^(2/9)) (Newey-West 1994 rule). + +ВЫХОД на район: `{best_lag_months, coef, r2, n, ...}` (см. DistributedLagFit). + • best_lag_months — лаг с пиком |β_j| из ОЦЕНЁННОЙ Almon-формы (момент сильнейшего + отклика; для фразы §9.6 «через Y мес»). + • coef — long-run (кумулятивный) мультипликатор Σ_j β_j на Δln: суммарный + %-эффект от удержания ставки на +1 п.п. (а не отклик одного месяца). Документ: + coef = долгосрочный мультипликатор; per_lag_coef несёт всю форму. + +GATE (зеркало analytics_queries._elasticity_coef, L1826-1852: n≥30 ∧ R²≥0.1 ∧ +slope<0 иначе fallback). Адаптация — gate смотрит на ЗНАК long-run β (ЦБ ↑ставку → +спрос ↓ → Σβ<0); n<30 ИЛИ R²<0.1 → degrade (source='fallback', claim не делаем). +Дух forecasting-модулей: PURE/детерминированно, graceful-on-thin-data, без LLM. + +ADVISORY: §9.6-стек советующий (как rate_sensitivity). Модуль самостоятелен — +ПОДКЛЮЧЕНИЕ к §9.6-консьюмеру отложено (точка интеграции — в docstring +compute_district_rate_regression), чтобы не задеть рабочий best_lag-путь +(product_scoring / demand_normalization / demand_supply_forecast зовут +compute_rate_sensitivity). Зеркалит дисциплину #979 (ship module + tests + note). + +psycopg v3 / SQLAlchemy text: bind ВСЕГДА через CAST(:x AS type) — НИКОГДА :x::type. +""" + +from __future__ import annotations + +import logging +import math +from dataclasses import dataclass +from datetime import date +from typing import Any, Literal + +import numpy as np +from sqlalchemy.orm import Session + +from app.services.forecasting.macro_series import get_monthly_macro +from app.services.forecasting.rate_sensitivity import _delta +from app.services.forecasting.sales_series import ( + SegmentSpec, + build_sales_series, + log_diff, +) + +logger = logging.getLogger(__name__) + +# ── Named-константы ─────────────────────────────────────────────────────────── + +# Глубина ряда по умолчанию (месяцев назад) — зеркалит _DEFAULT_MONTHS_BACK +# rate_sensitivity / macro_series (48 ≈ 4 года): §9.6 join-ит demand↔macro +# месяц-в-месяц, окна одной длины. +_DEFAULT_MONTHS_BACK: int = 48 + +# Максимальный лаг key_rate (мес). 0..6 — полугодовое окно отклика спроса на +# смену ставки (ипотека/сделки оформляются месяцами; полугодовой хвост ловит +# долгий эффект). Совпадает с верхней границей _LAGS rate_sensitivity (там {0,1,2,3,6}). +_MAX_LAG: int = 6 + +# Степень полинома Алмона над 7 лаговыми коэффициентами. 2 (квадратичная) — +# стандартный минимум, дающий «горб» (рост→пик→спад) отклика: реакция спроса +# нарастает, достигает максимума через несколько месяцев, затухает. deg<7 +# (=число лагов) — суть Алмона: 3 гладких параметра вместо 7 шумных коллинеарных. +_ALMON_DEGREE: int = 2 + +# GATE-пороги (зеркало _elasticity_coef L1856): n≥30 строк ∧ R²≥0.1 ∧ верный знак +# (для DL — long-run Σβ<0). Здесь одна «строка» = один Δln-МЕСЯЦ с полным набором +# лагов; на 48-мес окне их ≤ ~41 (минус _MAX_LAG на разогрев лагов и дыры). +_MIN_OBS: int = 30 +_MIN_R2: float = 0.1 + +# Минимум наблюдений, ниже которого Almon-OLS вообще не пытаемся (нужно > числа +# параметров d+1 с запасом на остаточные степени свободы для R²/HAC). 8 ≈ дух +# rate_sensitivity._MIN_OBS — но это лишь «можно ли фитить», НЕ gate-порог для +# claim (тот — _MIN_OBS=30 выше). +_MIN_FIT_OBS: int = 8 + +# Текст §9.6 (НЕ LLM) — шаблон фразы из оценённой лаговой формы. +_PHRASE_TEMPLATE: str = ( + "При росте ключевой ставки на 1 п.п. спрос снижается в среднем на {x}% " + "(пик эффекта через {y} мес.)." +) +_PHRASE_INSUFFICIENT: str = ( + "недостаточно данных для distributed-lag оценки чувствительности к ставке" +) + +# Survivorship-FREE помесячный агрегат сделок (зеркало rate_sensitivity._SOURCE_A); +# Literal — чтобы build_sales_series принял его как SalesSource без приведения. +_SOURCE_A: Literal["corpus_room_month"] = "corpus_room_month" + + +@dataclass(frozen=True) +class DistributedLagFit: + """Результат Almon distributed-lag регрессии Δln(demand) ~ Δrate[0..K]. + + Детерминированный. Числовые поля = None при недостатке данных / провале gate + (никогда 0-как-заглушка). `phrase` ВСЕГДА заполнена. ADVISORY до подключения. + + coef — LONG-RUN (кумулятивный) мультипликатор Σ_j β_j на Δln: суммарный + %-эффект (в exp-масштабе через x_pct) от удержания ставки на +1 п.п. Полную + ФОРМУ отклика по лагам несёт per_lag_coef; пик |β_j| → best_lag_months. + """ + + segment: dict[str, str | None] + best_lag_months: int | None # лаг пика |β_j| оценённой формы (момент сильнейшего отклика) + coef: float | None # long-run Σ_j β_j на Δln (кумулятивный мультипликатор) + x_pct: float | None # 100·(exp(coef)−1): %-эффект на +1 п.п. (NEGATIVE при ↓) + r2: float | None # R² distributed-lag регрессии + n: int # число использованных наблюдений (полных Δln-месяцев с лагами) + per_lag_coef: tuple[float, ...] | None # β_0..β_K из Almon-формы (вся форма отклика) + hac_se: tuple[float, ...] | None # Newey-West SE для β_0..β_K (ручной NW) + hac_bandwidth: int | None # окно Бартлетта L, на котором считались HAC SE + almon_degree: int # степень полинома Алмона + source: str # 'regression' (gate пройден) | 'fallback' (degrade) + phrase: str + + def as_dict(self) -> dict[str, Any]: + return { + "segment": dict(self.segment), + "best_lag_months": self.best_lag_months, + "coef": _round_or_none(self.coef, 4), + "x_pct": _round_or_none(self.x_pct, 1), + "r2": _round_or_none(self.r2, 4), + "n": self.n, + "per_lag_coef": ( + [round(c, 4) for c in self.per_lag_coef] if self.per_lag_coef is not None else None + ), + "hac_se": ([round(s, 4) for s in self.hac_se] if self.hac_se is not None else None), + "hac_bandwidth": self.hac_bandwidth, + "almon_degree": self.almon_degree, + "source": self.source, + "phrase": self.phrase, + } + + +def _round_or_none(value: float | None, digits: int) -> float | None: + return round(value, digits) if value is not None else None + + +# ────────────────────────────────────────────────────────────────────────────── +# Pure-логика — без БД, полностью юнит-тестируемо (numpy на синтетике с известным лагом). +# ────────────────────────────────────────────────────────────────────────────── + + +def _build_lag_matrix( + x: list[float | None], y: list[float | None], *, max_lag: int +) -> tuple[np.ndarray, np.ndarray] | None: + """Собрать матрицу лагов регрессора и выровненный y, дропнув неполные строки. + + Для каждого месяца t строим вектор регрессоров [x[t], x[t−1], …, x[t−max_lag]] + и спариваем с y[t]. Строку используем ТОЛЬКО если y[t] и ВСЕ max_lag+1 лаговых + значений конечны (None/NaN/Inf в любом лаге → строку дропаем: distributed-lag + требует полный лаговый профиль, частичный сместил бы оценку). PURE, без БД. + + Args: + x: регрессор по месяцам (обычно Δrate), None-дыры ок. + y: зависимая (Δln(demand)) по тем же месяцам, None-дыры ок. + max_lag: максимальный лаг (включительно) — матрица имеет max_lag+1 столбец. + + Returns: + (X, yv): X формы (n, max_lag+1) [lag0..lagK], yv формы (n,). None если + ни одной полной строки (n=0). + """ + n_months = min(len(x), len(y)) + rows: list[list[float]] = [] + ys: list[float] = [] + for t in range(max_lag, n_months): + yv = y[t] + if yv is None: + continue + yf = float(yv) + if not math.isfinite(yf): + continue + lagvec: list[float] = [] + ok = True + for j in range(max_lag + 1): + xv = x[t - j] + if xv is None: + ok = False + break + xf = float(xv) + if not math.isfinite(xf): + ok = False + break + lagvec.append(xf) + if not ok: + continue + rows.append(lagvec) + ys.append(yf) + + if not rows: + return None + return np.asarray(rows, dtype=float), np.asarray(ys, dtype=float) + + +def _almon_basis(max_lag: int, degree: int) -> np.ndarray: + """Матрица Almon-весов W формы (max_lag+1, degree+1): W[j, p] = j^p. + + β_j = Σ_p γ_p · j^p = (W @ γ)[j]. Преобразование регрессоров: если X — матрица + лагов (n, max_lag+1), то Almon-регрессоры Z = X @ W (n, degree+1), и OLS y~Z + даёт γ; per-lag β = W @ γ. j^0 столбец = 1 (intercept полинома → β-уровень). + PURE. + """ + lags = np.arange(max_lag + 1, dtype=float) + return np.vander(lags, N=degree + 1, increasing=True) + + +def _ols_lstsq(z: np.ndarray, y: np.ndarray) -> tuple[np.ndarray, np.ndarray] | None: + """OLS y ~ [1, Z] через numpy lstsq → (coef_with_intercept, residuals). PURE. + + Добавляем столбец-константу (свободный член регрессии — НЕ Almon-уровень). + Возвращает None при недостатке наблюдений (n ≤ #параметров → нет степеней + свободы) или вырожденной (rank-deficient) матрице плана. + + Returns: + (coef, resid): coef[0] = intercept, coef[1:] = γ; resid = y − ŷ. None если + фит невозможен. + """ + n = z.shape[0] + design = np.column_stack([np.ones(n), z]) + k = design.shape[1] + if n <= k: + return None + # rank-проверка: коллинеарный план → оценка γ не определена однозначно. + if np.linalg.matrix_rank(design) < k: + return None + coef, _res, _rank, _sv = np.linalg.lstsq(design, y, rcond=None) + resid = y - design @ coef + return coef, resid + + +def _r2(y: np.ndarray, resid: np.ndarray) -> float | None: + """R² = 1 − SS_res/SS_tot. None при нулевой дисперсии y (R² не определён). PURE.""" + ss_res = float(np.sum(resid**2)) + ss_tot = float(np.sum((y - float(np.mean(y))) ** 2)) + if ss_tot == 0.0: + return None + return 1.0 - ss_res / ss_tot + + +def newey_west_bandwidth(n: int) -> int: + """Окно Бартлетта L для Newey-West по правилу floor(4·(n/100)^(2/9)). + + Newey-West (1994) automatic bandwidth: растёт с n, но медленно. Нижняя + граница 1 (нужна хотя бы lag-1 автоковариация, иначе HAC == обычная OLS). + PURE. + """ + if n <= 1: + return 0 + bw: int = math.floor(4.0 * (n / 100.0) ** (2.0 / 9.0)) + return max(1, bw) + + +def newey_west_cov(x_design: np.ndarray, resid: np.ndarray, *, bandwidth: int) -> np.ndarray: + """HAC (Newey-West) ковариация оценок OLS — РУЧНОЙ numpy, БЕЗ statsmodels. + + Heteroskedasticity-and-autocorrelation-consistent ковариация: + V = (X'X)^{-1} · S · (X'X)^{-1}, S = Σ_0 + Σ_{l=1..L} w_l (Σ_l + Σ_l'), + где Σ_0 = Σ_t u_t² x_t x_t' (мясо White/HC0), Σ_l = Σ_t u_t u_{t−l} x_t x_{t−l}' + (lag-l автоковариация моментов), w_l = 1 − l/(L+1) — вес Бартлетта (гарантирует + положительную полуопределённость S). u_t — остатки, x_t — строка плана. + + Δln-остатки месячного спроса автокоррелированы (перекрытие лаговых окон + + инерция) → обычные OLS-SE занижены; NW их корректирует. PURE, без БД. + + Args: + x_design: матрица плана (n, k) — та же [1, Z], что в _ols_lstsq. + resid: остатки OLS (n,). + bandwidth: окно Бартлетта L (≥0); 0 → только Σ_0 (== White HC0). + + Returns: + Ковариационная матрица (k, k). При вырожденной X'X — псевдообратная + (graceful, не crash). + """ + n = x_design.shape[0] + xtx = x_design.T @ x_design + try: + xtx_inv = np.linalg.inv(xtx) + except np.linalg.LinAlgError: # вырожденная X'X → псевдообратная (graceful) + xtx_inv = np.linalg.pinv(xtx) + + u = resid.reshape(-1, 1) + ux = x_design * u # (n, k): строка t = u_t · x_t + # Σ_0 — «мясо» White (HC0). + s = ux.T @ ux + l_bw = max(0, min(bandwidth, n - 1)) # эффективное окно Бартлетта (L) + for lag in range(1, l_bw + 1): + w = 1.0 - lag / (l_bw + 1.0) # вес Бартлетта + gamma = ux[lag:].T @ ux[:-lag] # Σ_t u_t u_{t-l} x_t x_{t-l}' + s = s + w * (gamma + gamma.T) + cov: np.ndarray = xtx_inv @ s @ xtx_inv + return cov + + +def _peak_lag(per_lag: np.ndarray) -> int: + """Лаг с максимальным |β_j| — момент сильнейшего отклика спроса. PURE. + + При нескольких равных максимумах берём ПЕРВЫЙ (наименьший лаг — самый ранний + пик). Для фразы §9.6 «через Y мес». + """ + return int(np.argmax(np.abs(per_lag))) + + +def _x_pct_from_coef(coef: float) -> float: + """Long-run β на Δln → %-эффект на +1 п.п. ставки: 100·(exp(β)−1). PURE. + + β<0 → отрицательный % (спрос падает). exp т.к. Y=Δln (мультипликативный + масштаб). Зеркало rate_sensitivity._x_pct_from_beta. + """ + return 100.0 * (math.exp(coef) - 1.0) + + +def fit_almon_dl( + x: list[float | None], + y: list[float | None], + *, + max_lag: int = _MAX_LAG, + degree: int = _ALMON_DEGREE, + min_obs: int = _MIN_OBS, + min_r2: float = _MIN_R2, +) -> dict[str, Any] | None: + """Almon polynomial distributed-lag фит Δln(demand) ~ Δrate[0..max_lag]. PURE. + + Шаги: + 1. Собрать матрицу лагов X (n, max_lag+1) + выровненный y, дропнув строки с + неполным лаговым профилем (_build_lag_matrix). + 2. Almon-преобразование Z = X @ W (W[j,p]=j^p) → degree+1 регрессоров вместо + max_lag+1 коллинеарных. OLS y~[1,Z] (lstsq) → intercept + γ. + 3. Реконструкция per-lag β = W @ γ; long-run Σβ; R²; пик-лаг. + 4. HAC (Newey-West) ковариация на плане [1,Z] → SE для γ; SE для per-lag β + через delta-method (J = [0|W], V_β = J·V_γ·J'). + + Возвращает dict (см. ключи ниже) либо None если фит невозможен (n < запас по + степеням свободы / вырожденный план / нулевая дисперсия y). GATE (n≥min_obs ∧ + R²≥min_r2 ∧ Σβ<0) здесь НЕ применяется — это делает оркестратор (чтобы pure-фит + был переиспользуем и для диагностики «почти прошёл»). + + Args: + x: регрессор по месяцам (Δrate), None-дыры ок. + y: зависимая (Δln(demand)) по тем же месяцам, None-дыры ок. + max_lag: макс. лаг (включительно). + degree: степень полинома Алмона (< max_lag+1). + min_obs / min_r2: пороги для УДОБСТВА вызывающего (возвращаются в dict как + gate_n_ok / gate_r2_ok), сам фит ими не отсекается. + + Returns: + dict с ключами: per_lag_coef (tuple), long_run_coef (float), best_lag (int), + r2 (float|None), n (int), hac_se (tuple), hac_bandwidth (int), + intercept (float), gate_n_ok (bool), gate_r2_ok (bool), gate_sign_ok (bool). + None если фит математически невозможен. + """ + if degree >= max_lag + 1: + # Полином не должен иметь параметров ≥ числа лагов — иначе это не + # ограничение (вырождается в free-lags). Защита от мисконфига. + logger.warning("fit_almon_dl: degree=%d >= max_lag+1=%d — refusing", degree, max_lag + 1) + return None + + built = _build_lag_matrix(x, y, max_lag=max_lag) + if built is None: + return None + x_lags, yv = built + n = int(x_lags.shape[0]) + if n < _MIN_FIT_OBS: + return None + # Нулевая дисперсия зависимой → R²/наклоны не определены. + if float(np.var(yv)) == 0.0: + return None + + w = _almon_basis(max_lag, degree) # (max_lag+1, degree+1) + z = x_lags @ w # (n, degree+1) — Almon-регрессоры + fit = _ols_lstsq(z, yv) + if fit is None: + return None + coef, resid = fit # coef = [intercept, γ_0..γ_d] + intercept = float(coef[0]) + gamma = coef[1:] + per_lag = w @ gamma # (max_lag+1,) — реконструированные β_j + long_run = float(np.sum(per_lag)) + r2 = _r2(yv, resid) + best_lag = _peak_lag(per_lag) + + # HAC (Newey-West) на ТОМ ЖЕ плане [1, Z]. + design = np.column_stack([np.ones(n), z]) + bw = newey_west_bandwidth(n) + cov = newey_west_cov(design, resid, bandwidth=bw) + # SE per-lag β через delta-method: β = W·γ = J·coef, J = [0_col | W] (intercept + # не входит в β). V_β = J·cov·J'; диагональ ≥0 → sqrt (отрицательные FP-края → 0). + j = np.column_stack([np.zeros((max_lag + 1, 1)), w]) # (max_lag+1, degree+2) + cov_beta = j @ cov @ j.T + var_beta = np.clip(np.diag(cov_beta), a_min=0.0, a_max=None) + hac_se = tuple(float(s) for s in np.sqrt(var_beta)) + + gate_n_ok = n >= min_obs + gate_r2_ok = r2 is not None and r2 >= min_r2 + gate_sign_ok = long_run < 0.0 + + return { + "per_lag_coef": tuple(float(c) for c in per_lag), + "long_run_coef": long_run, + "best_lag": best_lag, + "r2": r2, + "n": n, + "hac_se": hac_se, + "hac_bandwidth": bw, + "intercept": intercept, + "gate_n_ok": gate_n_ok, + "gate_r2_ok": gate_r2_ok, + "gate_sign_ok": gate_sign_ok, + } + + +def _build_phrase(*, x_pct: float | None, best_lag: int | None, gated: bool) -> str: + """Фраза §9.6 (НЕ LLM) из оценённой лаговой формы. PURE. + + gate провален / нет valid эффекта → «недостаточно данных…». Иначе: + «при росте ставки +1 п.п. спрос снижается на X% (пик через Y мес)». X — + положительная МАГНИТУДА %-эффекта (long-run), Y — пик-лаг. + """ + if not gated or x_pct is None or best_lag is None: + return _PHRASE_INSUFFICIENT + return _PHRASE_TEMPLATE.format(x=round(abs(x_pct), 1), y=best_lag) + + +def _insufficient(segment: dict[str, str | None], *, n: int = 0) -> DistributedLagFit: + """Граничный результат «недостаточно данных» (fallback, фраза-заглушка). PURE.""" + return DistributedLagFit( + segment=segment, + best_lag_months=None, + coef=None, + x_pct=None, + r2=None, + n=n, + per_lag_coef=None, + hac_se=None, + hac_bandwidth=None, + almon_degree=_ALMON_DEGREE, + source="fallback", + phrase=_PHRASE_INSUFFICIENT, + ) + + +def build_fit_result( + x: list[float | None], + y: list[float | None], + *, + segment: dict[str, str | None], + max_lag: int = _MAX_LAG, + degree: int = _ALMON_DEGREE, + min_obs: int = _MIN_OBS, + min_r2: float = _MIN_R2, +) -> DistributedLagFit: + """Прогнать Almon-DL фит и обернуть в DistributedLagFit с gate-деградацией. PURE. + + GATE (зеркало _elasticity_coef): n≥min_obs ∧ R²≥min_r2 ∧ long-run Σβ<0 → + source='regression' (claim). Иначе → degrade: source='fallback', фраза + «недостаточно данных», но per_lag_coef/r2/n СОХРАНЯЕМ для диагностики (как + _elasticity_coef возвращает r2/n в fallback). НЕ crash на тонких данных. + + Это чистая обёртка (без БД) — тестируется на синтетике с известным лагом. + """ + fit = fit_almon_dl(x, y, max_lag=max_lag, degree=degree, min_obs=min_obs, min_r2=min_r2) + if fit is None: + return _insufficient(segment) + + n = int(fit["n"]) + gated = bool(fit["gate_n_ok"] and fit["gate_r2_ok"] and fit["gate_sign_ok"]) + long_run = float(fit["long_run_coef"]) + r2 = fit["r2"] + per_lag = tuple(fit["per_lag_coef"]) + hac_se = tuple(fit["hac_se"]) + best_lag = int(fit["best_lag"]) + + if not gated: + # Degrade: сохраняем числа для диагностики, но source='fallback' и фраза-заглушка. + logger.info( + "regression: gate failed (segment=%s n=%d r2=%s long_run=%.4f " + "n_ok=%s r2_ok=%s sign_ok=%s) → fallback", + segment, + n, + None if r2 is None else round(r2, 4), + long_run, + fit["gate_n_ok"], + fit["gate_r2_ok"], + fit["gate_sign_ok"], + ) + return DistributedLagFit( + segment=segment, + best_lag_months=None, + coef=None, + x_pct=None, + r2=_round_or_none(r2, 4), + n=n, + per_lag_coef=per_lag, + hac_se=hac_se, + hac_bandwidth=int(fit["hac_bandwidth"]), + almon_degree=degree, + source="fallback", + phrase=_PHRASE_INSUFFICIENT, + ) + + x_pct = _x_pct_from_coef(long_run) + phrase = _build_phrase(x_pct=x_pct, best_lag=best_lag, gated=True) + logger.info( + "regression(OK): segment=%s long_run=%.4f x_pct=%.1f best_lag=%d r2=%.4f n=%d bw=%d", + segment, + long_run, + x_pct, + best_lag, + r2 if r2 is not None else float("nan"), + n, + int(fit["hac_bandwidth"]), + ) + return DistributedLagFit( + segment=segment, + best_lag_months=best_lag, + coef=long_run, + x_pct=x_pct, + r2=r2, + n=n, + per_lag_coef=per_lag, + hac_se=hac_se, + hac_bandwidth=int(fit["hac_bandwidth"]), + almon_degree=degree, + source="regression", + phrase=phrase, + ) + + +# ────────────────────────────────────────────────────────────────────────────── +# DB-оркестратор — тонкий, graceful. Pure-логика выше тестируется без него. +# ────────────────────────────────────────────────────────────────────────────── + + +def _align_demand_deltas( + sales_months: list[date], sales_units: list[int], macro_months: list[date] +) -> list[float | None]: + """Выровнять Δln(units) спроса по сетке макро-месяцев (общая ось X↔Y). + + Зеркало rate_sensitivity._align_sales_deltas: log_diff даёт Δln по сетке + ПРОДАЖ, перекладываем на macro_months (месяц без продаж → None), чтобы пары + (Δrate[t−L], Δln[t]) были month-в-month и лаговая матрица строилась по единой + временной оси. PURE. + + NB (#979 дух): дессзонивание ПЕРЕД log_diff здесь НЕ применяется — та же + оговорка, что в rate_sensitivity._align_sales_deltas (на коротком ряде ratio- + to-mean фактор смещает восстановленный лаг). Отложено той же задачей. + """ + deltas = log_diff(sales_units) + by_month = dict(zip(sales_months, deltas, strict=False)) + return [by_month.get(m) for m in macro_months] + + +def compute_district_rate_regression( + db: Session, + *, + district: str, + obj_class: str | None = None, + months_back: int = _DEFAULT_MONTHS_BACK, + max_lag: int = _MAX_LAG, + degree: int = _ALMON_DEGREE, +) -> DistributedLagFit: + """§9.6 Almon distributed-lag регрессия месячного спроса РАЙОНА на key_rate. + + Constrained DL (Almon, deg `degree`) Δln(demand_district) ~ Δkey_rate при лагах + 0..max_lag, с реконструкцией per-lag β и HAC (Newey-West) SE. GATE зеркалит + _elasticity_coef (n≥30 ∧ R²≥0.1 ∧ long-run Σβ<0 иначе fallback). ДЕТЕРМИНИРОВАНО. + + Данные: + • key_rate — get_monthly_macro (PR2), Δ первой разностью (_delta) → X-ось. + • спрос района — build_sales_series Source A (objective_corpus_room_month, + survivorship-FREE помесячный агрегат сделок), Δln (log_diff) → Y-ось, + выровненная на сетку макро (_align_demand_deltas). + + Graceful-on-thin-data: пустой/тонкий ряд / провал gate → source='fallback', + фраза «недостаточно данных…», НЕ crash (дух forecasting-модулей). + + ADVISORY + НЕ ПОДКЛЮЧЕНО (отложенная интеграция, #978 Part B): + Точка интеграции в §9.6 — там же, где сейчас зовётся compute_rate_sensitivity + (product_scoring._build → mortgage_sensitivity; demand_normalization; + demand_supply_forecast explain-фраза). Подключение ОТЛОЖЕНО, чтобы не задеть + рабочий single-lag best_lag-путь (риск регресса в трёх консьюмерах). Зеркалит + дисциплину #979: ship module + tests + note integration point. §9.6-стек + advisory в любом случае. + + Args: + db: SQLAlchemy sync Session. + district: район ЕКБ (Source A column `district`). + obj_class: класс ЖК (None → агрегат по району); регистр нормализуется в SQL. + months_back: глубина ряда (по умолчанию 48). + max_lag / degree: окно лагов и степень Алмона. + + Returns: + DistributedLagFit (всегда; фраза заполнена даже при нехватке данных). + """ + segment: dict[str, str | None] = {"district": district, "obj_class": obj_class} + + macro = get_monthly_macro(db, months_back=months_back) + rate_deltas = _delta([m.key_rate for m in macro]) + macro_months = [m.month for m in macro] + + spec = SegmentSpec(obj_class=obj_class, district=district) + sales = build_sales_series(db, spec=spec, source=_SOURCE_A, months_back=months_back) + demand_deltas = _align_demand_deltas(sales.months, sales.units, macro_months) + + return build_fit_result( + rate_deltas, + demand_deltas, + segment=segment, + max_lag=max_lag, + degree=degree, + ) diff --git a/backend/scripts/backtest_rate_sensitivity.py b/backend/scripts/backtest_rate_sensitivity.py index 6791a4f3..1b20345c 100644 --- a/backend/scripts/backtest_rate_sensitivity.py +++ b/backend/scripts/backtest_rate_sensitivity.py @@ -74,6 +74,9 @@ could be an ARTIFACT rather than a true ``no signal``. We add two controls: ``ln(units) ~ a + b·month_index`` and subtract it, regressing the Δ of the RESIDUALS vs Δrate. A spurious monotone survivorship trend lands almost entirely in ``b`` and is removed, so it can no longer drive the regression. + The trend is fit on the TRAIN months ONLY and projected point-in-time onto + the test months (#978 Part A: fitting it on train+test together leaks future + info into the test residuals and inflates the detrended OOS hit-rate). Read both alongside Source B raw: if Source B DETRENDED still shows no OOS signal AND survivorship-free Source A agrees (thin caveat aside), the engine's @@ -301,7 +304,9 @@ def _rate_first_diff(rate_levels: list[float | None]) -> list[float | None]: return out -def _detrend_log(values: list[float | int | None]) -> list[float | None]: +def _detrend_log( + values: list[float | int | None], *, fit_n: int | None = None +) -> list[float | None]: """Linear-detrend the LOG of a units series → log-residuals. PURE (no DB). The survivorship control for #978b. We: @@ -318,9 +323,21 @@ def _detrend_log(values: list[float | int | None]) -> list[float | None]: driven by it. The caller differences these residuals (they are already in log space) instead of calling ``log_diff`` again. - Below ``_DETREND_MIN_POINTS`` finite points a line is not identifiable, so we - PASS THROUGH the log values unchanged (residual == log value); differencing - them then reproduces the raw ``log_diff`` path exactly. PURE. + LOOK-AHEAD LEAKAGE GUARD (#978 reopen, Part A): when ``fit_n`` is given the + trend ``(a, b)`` is estimated ONLY on the finite points among the first + ``fit_n`` months (the TRAIN slice), then PROJECTED point-in-time onto every + month (test residual = ``ln(units_t) − (a + b·t)`` with TRAIN-fitted a,b). + This is mandatory in the backtest: fitting the trend on train+test together + lets the held-out TEST observations shape ``(a, b)``, so the test residuals + embed future information and the detrended OOS hit-rate is inflated. With + ``fit_n=None`` (default) the trend is fit on the whole finite series — only + safe for a non-holdout, full-sample descriptive detrend, NEVER for OOS + scoring. + + Below ``_DETREND_MIN_POINTS`` finite points (counted within the fit window + when ``fit_n`` is set) a line is not identifiable, so we PASS THROUGH the log + values unchanged (residual == log value); differencing them then reproduces + the raw ``log_diff`` path exactly. PURE. """ logs: list[float | None] = [] for v in values: @@ -330,18 +347,23 @@ def _detrend_log(values: list[float | int | None]) -> list[float | None]: vf = float(v) logs.append(math.log(vf) if vf > 0 else None) - finite_idx = [i for i, lv in enumerate(logs) if lv is not None] - if len(finite_idx) < _DETREND_MIN_POINTS: - return logs # not enough points to fit a trend → passthrough of logs + # The trend is fit only on finite points whose index is inside the fit window + # (TRAIN slice [0:fit_n] when fit_n is given; the whole series otherwise). + fit_upper = len(logs) if fit_n is None else max(0, fit_n) + fit_idx = [i for i, lv in enumerate(logs) if lv is not None and i < fit_upper] + if len(fit_idx) < _DETREND_MIN_POINTS: + return logs # not enough TRAIN points to fit a trend → passthrough of logs - xs = np.array([float(i) for i in finite_idx], dtype=float) - ys = np.array([float(logs[i]) for i in finite_idx], dtype=float) # type: ignore[arg-type] + xs = np.array([float(i) for i in fit_idx], dtype=float) + ys = np.array([float(logs[i]) for i in fit_idx], dtype=float) # type: ignore[arg-type] # Degenerate x-variance (all same index — impossible for ≥3 distinct idx but # guard anyway) → no trend to remove, passthrough. if float(np.ptp(xs)) == 0.0: return logs slope, intercept = np.polyfit(xs, ys, 1) + # Project the TRAIN-fitted (a, b) onto EVERY month, incl. the held-out test + # months — strictly point-in-time, no test observation entered the fit. out: list[float | None] = [] for i, lv in enumerate(logs): if lv is None: @@ -495,7 +517,9 @@ def align_series( return months, units, rates -def _delta_sales_series(units: list[int], *, detrend: bool) -> list[float | None]: +def _delta_sales_series( + units: list[int], *, detrend: bool, fit_n: int | None = None +) -> list[float | None]: """Build the Δ(log-units) regressand for one tier. PURE (deferred import). Two variants, both ending in a Δ of log-space values ``evaluate_oos`` scores: @@ -505,11 +529,17 @@ def _delta_sales_series(units: list[int], *, detrend: bool) -> list[float | None (``_detrend_log``), THEN first-difference the residuals. We difference the residuals DIRECTLY (they are already in log space) rather than ``log_diff`` (which would re-take logs of residuals that may be ≤0). + + ``fit_n`` (the TRAIN month count) is forwarded to ``_detrend_log`` so the + detrend trend is estimated on TRAIN months only and projected point-in-time + onto the test months — the #978 Part A look-ahead-leakage fix. It is ignored + on the non-detrend path (``log_diff`` is already point-in-time: a first + difference reads only t and t−1). """ if not detrend: _bl, _ols, log_diff = _import_engine() return log_diff(units) - return _rate_first_diff(_detrend_log(units)) + return _rate_first_diff(_detrend_log(units, fit_n=fit_n)) def backtest_tier( @@ -553,7 +583,14 @@ def backtest_tier( skipped=f"only {n_aligned} aligned months (< {min_months})", ) - delta_sales = _delta_sales_series(units, detrend=detrend) + # Detrend (when enabled) must be fit on TRAIN months ONLY, then projected + # point-in-time onto the test months — otherwise the held-out TEST data + # shapes the trend and the OOS hit-rate is inflated by look-ahead leakage + # (#978 Part A). We compute the SAME train boundary evaluate_oos will use + # (len(delta_sales) == len(units) == n_aligned, so the split index matches) + # and pass it as the detrend fit window. + n_train = _time_ordered_split(n_aligned, holdout_frac) + delta_sales = _delta_sales_series(units, detrend=detrend, fit_n=n_train) rate_deltas = _rate_first_diff([float(r) for r in rates]) res = evaluate_oos(delta_sales, rate_deltas, holdout_frac=holdout_frac) @@ -1156,7 +1193,8 @@ def render_table(results: dict[str, Any]) -> str: if detrended: lines.append( "DETRENDED: ln(units) linearly detrended (residuals) BEFORE differencing — " - "removes a spurious monotone (survivorship) trend so it can't drive β." + "removes a spurious monotone (survivorship) trend so it can't drive β. " + "Trend fit on TRAIN months only, projected point-in-time onto test (no leakage)." ) if results.get("a_district_ignored"): lines.append( diff --git a/backend/tests/scripts/test_backtest_rate_sensitivity.py b/backend/tests/scripts/test_backtest_rate_sensitivity.py index 0050f700..cb3f2ece 100644 --- a/backend/tests/scripts/test_backtest_rate_sensitivity.py +++ b/backend/tests/scripts/test_backtest_rate_sensitivity.py @@ -254,6 +254,115 @@ class TestDetrendLog: assert len(bt._detrend_log(vals)) == 10 +# --------------------------------------------------------------------------- # +# Look-ahead leakage fix (#978 Part A) — detrend trend fit on TRAIN months only, +# projected point-in-time onto test (never fit on train+test together). +# --------------------------------------------------------------------------- # + + +class TestDetrendNoLeakage: + def test_train_only_fit_matches_manual_polyfit_on_train_slice(self) -> None: + # With fit_n given, the trend (a, b) must be the polyfit of ONLY the + # finite points in [0:fit_n] — the test months must not enter the fit. + n, fit_n = 30, 20 + units = [max(1, round(math.exp(6.0 + 0.05 * t))) for t in range(n)] + logs = [math.log(u) for u in units] + # Manual train-only line. + xs = list(range(fit_n)) + ys = logs[:fit_n] + b, a = bt.np.polyfit(bt.np.array(xs, dtype=float), bt.np.array(ys, dtype=float), 1) + resid = bt._detrend_log(units, fit_n=fit_n) + # Every residual equals ln(u_t) − (a + b·t) with the TRAIN-fitted line, + # INCLUDING the test months (the line is projected forward, not refit). + for t in range(n): + assert resid[t] is not None + assert math.isclose(resid[t], logs[t] - (a + b * t), abs_tol=1e-9) # type: ignore[arg-type] + + def test_test_points_do_not_shape_the_trend(self) -> None: + # A BROKEN trend: gentle slope on train, steep slope on test. A full-sample + # (leaky) fit is pulled UP by the steep test tail; a train-only fit is not. + # So the residual at the LAST month must differ between the two — proving + # the test observations leak into the leaky fit but not the train-only one. + n, fit_n = 24, 16 + units: list[int] = [] + for t in range(n): + slope = 0.02 if t < fit_n else 0.20 # trend break at fit_n + base = 0.02 * min(t, fit_n) + extra = 0.20 * max(0, t - fit_n) + units.append(max(1, round(math.exp(6.0 + base + extra)) if t else round(math.exp(6.0)))) + _ = slope + leaky = bt._detrend_log(units) # fit_n=None → fit on train+test (leaks) + safe = bt._detrend_log(units, fit_n=fit_n) + # Last test month residual differs → the steep tail moved the leaky line + # but not the train-only line. + assert leaky[-1] is not None and safe[-1] is not None + assert abs(leaky[-1] - safe[-1]) > 0.05 # type: ignore[operator] + + def test_fit_n_gates_passthrough_on_train_point_count(self) -> None: + # Plenty of finite points overall, but only 2 (< _DETREND_MIN_POINTS) fall + # inside the TRAIN window → a line is not identifiable on TRAIN → passthrough + # of the logs (residual == ln(value)), exactly like the raw log_diff path. + units = [10, 20] + [30 + i for i in range(10)] # 12 finite, fit_n=2 + resid = bt._detrend_log(units, fit_n=2) + assert resid[0] is not None and math.isclose(resid[0], math.log(10)) + assert resid[1] is not None and math.isclose(resid[1], math.log(20)) + # Passthrough applies to ALL positions (no trend was removed anywhere). + assert resid[2] is not None and math.isclose(resid[2], math.log(30)) + + def test_backtest_tier_detrend_fits_train_only(self) -> None: + # End-to-end: backtest_tier must pass n_train as fit_n. We assert the + # detrended regressand it builds equals the one from a TRAIN-only detrend, + # and is NOT equal to the leaky full-sample detrend (when they differ). + n = 40 + ms = _months(n) + # Trend-confounded units with a real lag-2 signal (#978b-style series). + rate = _zero_drift_rate_levels(n, seed=5) + units = _units_from_rate_with_trend(rate, lag=2, beta=-0.06, trend_per_month=0.09) + sales = {ms[i]: units[i] for i in range(n)} + rate_by = {ms[i]: rate[i] for i in range(n)} + + # What backtest_tier should build internally (train-only fit). + n_train = bt._time_ordered_split(n, 0.7) + expected = bt._delta_sales_series(units, detrend=True, fit_n=n_train) + leaky = bt._delta_sales_series(units, detrend=True, fit_n=None) + + # Run the tier and reconstruct its regressand path via the same helper to + # confirm n_train is threaded through (the public API has no hook, so we + # assert the train-only and full-sample series genuinely differ — i.e. the + # fix is observable — and that the tier still produces a scored result). + res = bt.backtest_tier(sales, rate_by, tier=bt._EKB_WIDE, detrend=True, holdout_frac=0.7) + assert res.skipped is None + assert res.detrended is True + # The two regressands must differ somewhere in the test region (leakage is + # observable), so the train-only fix is a real behavioural change. + assert any( + e is not None and lk is not None and abs(e - lk) > 1e-9 + for e, lk in zip(expected[n_train:], leaky[n_train:], strict=False) + ) + + def test_no_leakage_oos_hit_rate_not_above_leaky(self) -> None: + # The core claim: look-ahead leakage INFLATES the detrended OOS hit-rate. + # On a trend-confounded series, the train-only (correct) detrend must give + # an OOS hit-rate ≤ the leaky full-sample detrend. We compare evaluate_oos + # on both regressands over the SAME aligned series. + n = 48 + rate = _zero_drift_rate_levels(n, seed=11) + units = _units_from_rate_with_trend(rate, lag=2, beta=-0.05, trend_per_month=0.07) + rate_deltas = bt._rate_first_diff(rate) + n_train = bt._time_ordered_split(n, 0.7) + + safe_sales = bt._delta_sales_series(units, detrend=True, fit_n=n_train) + leaky_sales = bt._delta_sales_series(units, detrend=True, fit_n=None) + safe = bt.evaluate_oos(safe_sales, rate_deltas, holdout_frac=0.7) + leaky = bt.evaluate_oos(leaky_sales, rate_deltas, holdout_frac=0.7) + + # Both should find a gated lag here; if either is None the inequality is + # vacuously fine (no inflation possible). When both score, leakage may only + # help (or tie) the leaky run — it must never make the corrected run higher. + if safe["oos_hit_rate"] is not None and leaky["oos_hit_rate"] is not None: + assert safe["oos_hit_rate"] <= leaky["oos_hit_rate"] + 1e-9 + + class TestAlignSeries: def test_inner_join_by_month(self) -> None: ms = _months(4) diff --git a/backend/tests/services/forecasting/test_regression.py b/backend/tests/services/forecasting/test_regression.py new file mode 100644 index 00000000..27aac706 --- /dev/null +++ b/backend/tests/services/forecasting/test_regression.py @@ -0,0 +1,513 @@ +"""Unit tests for §9.6 Almon distributed-lag regression (Forgejo #978 Part B). + +Covers the PURE numpy logic on SYNTHETIC series with a KNOWN injected lag effect: + - _build_lag_matrix — full-row-only lag profile, drops incomplete/None rows + - _almon_basis — W[j,p] = j^p (constrains 7 lags to degree+1 params) + - newey_west_bandwidth — floor(4·(n/100)^(2/9)) rule, ≥1 floor + - newey_west_cov — HAC covariance differs from naive OLS; PSD; manual NW + - fit_almon_dl — recovers the injected best_lag + sign + long-run; R²; + per-lag reconstruction; HAC SEs computed + - build_fit_result — gate (n≥30 ∧ R²≥0.1 ∧ Σβ<0) → regression vs fallback; + fallback on thin n / weak R² / wrong sign (no crash) + - _build_phrase — §9.6 text from the lag shape; insufficient on no-gate + - compute_district_rate_regression — DB orchestrator wiring (mocked session) + +NO live DB: the orchestrator test injects a fake session + monkeypatched data +loaders. Set a dummy DATABASE_URL BEFORE importing so app.core.config.Settings +fail-fast doesn't trip (same pattern as test_rate_sensitivity.py). +""" + +from __future__ import annotations + +import datetime as dt +import math +import os + +os.environ.setdefault("DATABASE_URL", "postgresql+psycopg://test:test@localhost:5432/test") + +import numpy as np +import pytest + +from app.services.forecasting import regression as reg + +# --------------------------------------------------------------------------- # +# Synthetic-series helpers — inject a KNOWN distributed-lag effect +# --------------------------------------------------------------------------- # + + +def _aperiodic_rate_deltas(n: int, *, seed: int = 13) -> list[float]: + """Δrate series with APERIODIC (LCG) jitter → low autocorrelation across lags. + + A periodic regressor would let false lags compete with the injected one; an + LCG jitter keeps successive Δ weakly correlated so the true lag shape wins. + Finite from index 0 (the DL matrix builder drops incomplete leading rows). + """ + lvl = 10.0 + state = seed + levels: list[float] = [] + for _ in range(n): + state = (state * 1103515245 + 12345) % 2147483648 + lvl += 0.3 + (state / 2147483648.0 - 0.5) * 0.8 + levels.append(lvl) + return [0.0] + [levels[i] - levels[i - 1] for i in range(1, n)] + + +def _hump_beta(max_lag: int, *, peak: int, scale: float = 0.06) -> np.ndarray: + """A negative 'hump' lag shape peaking (in magnitude) at ``peak``. + + |β_j| = scale − 0.012·(j−peak)² (floored at 0.005), all signs negative — the + economically expected shape (rate ↑ → demand ↓, response builds then fades). + Representable approximately by an Almon deg-2 polynomial, so the fit recovers + the peak and long-run. + """ + betas: list[float] = [] + for j in range(max_lag + 1): + mag = scale - 0.012 * (j - peak) ** 2 + betas.append(-max(0.005, mag)) + return np.asarray(betas, dtype=float) + + +def _y_from_lag_shape( + x: list[float], beta: np.ndarray, *, max_lag: int, noise: float = 0.0, seed: int = 0 +) -> list[float | None]: + """y[t] = Σ_j β_j·x[t−j] (+ optional gaussian noise); y[t 0.0: + val += float(rng.normal(0.0, noise)) + y[t] = val + return y + + +# --------------------------------------------------------------------------- # +# _build_lag_matrix +# --------------------------------------------------------------------------- # + + +class TestBuildLagMatrix: + def test_shapes_and_full_rows_only(self) -> None: + x = [float(i) for i in range(10)] + y = [float(i) * 0.1 for i in range(10)] + built = reg._build_lag_matrix(x, y, max_lag=2) + assert built is not None + xm, yv = built + # First usable row is t=max_lag=2 → 10−2 = 8 rows, 3 lag columns. + assert xm.shape == (8, 3) + assert yv.shape == (8,) + # Row 0 corresponds to t=2: [x[2], x[1], x[0]] = [2,1,0]. + assert list(xm[0]) == [2.0, 1.0, 0.0] + + def test_drops_rows_with_none_in_any_lag(self) -> None: + x: list[float | None] = [0.0, 1.0, None, 3.0, 4.0, 5.0] + y: list[float | None] = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5] + built = reg._build_lag_matrix(x, y, max_lag=2) + assert built is not None + xm, _yv = built + # t=2 reads x[0..2] (has None) → dropped; t=3 reads x[1..3] (has None) → + # dropped; t=4 reads x[2..4] (has None) → dropped; t=5 reads x[3..5] OK. + assert xm.shape == (1, 3) + assert list(xm[0]) == [5.0, 4.0, 3.0] + + def test_drops_rows_with_none_y(self) -> None: + x = [float(i) for i in range(6)] + y: list[float | None] = [0.0, 0.1, None, 0.3, 0.4, 0.5] + built = reg._build_lag_matrix(x, y, max_lag=1) + assert built is not None + xm, yv = built + # t=2 has y=None → dropped. Usable t ∈ {1,3,4,5} → 4 rows. + assert xm.shape == (4, 2) + assert yv.shape == (4,) + + def test_returns_none_when_no_full_row(self) -> None: + x: list[float | None] = [None, None, None] + y: list[float | None] = [1.0, 2.0, 3.0] + assert reg._build_lag_matrix(x, y, max_lag=1) is None + + +# --------------------------------------------------------------------------- # +# _almon_basis +# --------------------------------------------------------------------------- # + + +class TestAlmonBasis: + def test_j_to_the_p(self) -> None: + w = reg._almon_basis(3, 2) # lags 0..3, degree 2 + assert w.shape == (4, 3) + # Column p = j^p: col0 = ones, col1 = j, col2 = j². + assert list(w[:, 0]) == [1.0, 1.0, 1.0, 1.0] + assert list(w[:, 1]) == [0.0, 1.0, 2.0, 3.0] + assert list(w[:, 2]) == [0.0, 1.0, 4.0, 9.0] + + def test_reconstruct_quadratic_beta_exactly(self) -> None: + # β_j = 2 − 0.5j + 0.1j² is degree-2 → W @ γ reproduces it for γ=[2,−0.5,0.1]. + w = reg._almon_basis(6, 2) + gamma = np.array([2.0, -0.5, 0.1]) + beta = w @ gamma + expected = np.array([2.0 - 0.5 * j + 0.1 * j**2 for j in range(7)]) + assert np.allclose(beta, expected) + + +# --------------------------------------------------------------------------- # +# newey_west_bandwidth / newey_west_cov +# --------------------------------------------------------------------------- # + + +class TestNeweyWestBandwidth: + def test_rule_values(self) -> None: + # floor(4·(n/100)^(2/9)). + assert reg.newey_west_bandwidth(100) == 4 + assert reg.newey_west_bandwidth(41) == 3 + assert reg.newey_west_bandwidth(50) == 3 + + def test_small_n_values(self) -> None: + # floor(4·(10/100)^(2/9)) = floor(2.398) = 2; n=20 → floor(2.40)·… = 2. + assert reg.newey_west_bandwidth(10) == 2 + assert reg.newey_west_bandwidth(20) == 2 + + def test_floor_at_one_for_tiny_n(self) -> None: + # n=2 → floor(4·0.02^(2/9)) ≈ floor(1.06) = 1, but the ≥1 floor guarantees + # at least a lag-1 autocovariance whenever n>1. + assert reg.newey_west_bandwidth(2) == 1 + assert reg.newey_west_bandwidth(3) == 1 + + def test_zero_for_degenerate(self) -> None: + assert reg.newey_west_bandwidth(1) == 0 + assert reg.newey_west_bandwidth(0) == 0 + + +class TestNeweyWestCov: + def test_psd_and_symmetric(self) -> None: + rng = np.random.default_rng(1) + n = 50 + design = np.column_stack([np.ones(n), rng.normal(size=(n, 2))]) + resid = rng.normal(size=n) + cov = reg.newey_west_cov(design, resid, bandwidth=4) + # Symmetric and positive semi-definite (Bartlett weights guarantee PSD). + assert np.allclose(cov, cov.T, atol=1e-10) + eig = np.linalg.eigvalsh(cov) + assert float(eig.min()) >= -1e-8 + + def test_bandwidth_zero_equals_white_hc0(self) -> None: + rng = np.random.default_rng(2) + n = 40 + design = np.column_stack([np.ones(n), rng.normal(size=(n, 1))]) + resid = rng.normal(size=n) + cov0 = reg.newey_west_cov(design, resid, bandwidth=0) + # HC0: (X'X)^-1 (Σ u² x x') (X'X)^-1 — reconstruct manually. + xtx_inv = np.linalg.inv(design.T @ design) + ux = design * resid.reshape(-1, 1) + hc0 = xtx_inv @ (ux.T @ ux) @ xtx_inv + assert np.allclose(cov0, hc0, atol=1e-12) + + def test_hac_differs_from_naive_under_autocorrelation(self) -> None: + # Construct strongly AUTOCORRELATED residuals → HAC SE must differ from + # the naive iid OLS SE (the whole point of NW). + n = 80 + rng = np.random.default_rng(3) + x = rng.normal(size=n) + design = np.column_stack([np.ones(n), x]) + # AR(1) residuals (ρ=0.7) → positive autocorrelation. + e = np.zeros(n) + for t in range(1, n): + e[t] = 0.7 * e[t - 1] + rng.normal(0, 1) + hac = reg.newey_west_cov(design, e, bandwidth=reg.newey_west_bandwidth(n)) + sigma2 = float(e @ e) / (n - design.shape[1]) + naive = sigma2 * np.linalg.inv(design.T @ design) + # The slope variance estimates must differ materially under autocorrelation. + assert not math.isclose(hac[1, 1], naive[1, 1], rel_tol=0.05) + + +# --------------------------------------------------------------------------- # +# fit_almon_dl — recover the injected lag shape +# --------------------------------------------------------------------------- # + + +class TestFitAlmonDl: + def test_recovers_injected_best_lag_and_sign(self) -> None: + n, max_lag = 60, 6 + x = _aperiodic_rate_deltas(n, seed=13) + beta = _hump_beta(max_lag, peak=2) + y = _y_from_lag_shape(x, beta, max_lag=max_lag, noise=0.002, seed=0) + fit = reg.fit_almon_dl(x, y, max_lag=max_lag, degree=2) + assert fit is not None + # Peak lag recovered. + assert fit["best_lag"] == 2 + # Long-run sign negative (rate ↑ → demand ↓) and close to the truth. + assert fit["long_run_coef"] < 0 + assert math.isclose(fit["long_run_coef"], float(beta.sum()), abs_tol=0.02) + # Clean injected signal → high R². + assert fit["r2"] is not None and fit["r2"] > 0.8 + # Gate flags all green on this clean, long, correctly-signed series. + assert fit["gate_n_ok"] and fit["gate_r2_ok"] and fit["gate_sign_ok"] + + def test_recovers_different_peak_lag(self) -> None: + # Shift the injected peak to lag 4 → the fit must track it. + n, max_lag = 64, 6 + x = _aperiodic_rate_deltas(n, seed=21) + beta = _hump_beta(max_lag, peak=4) + y = _y_from_lag_shape(x, beta, max_lag=max_lag, noise=0.002, seed=1) + fit = reg.fit_almon_dl(x, y, max_lag=max_lag, degree=2) + assert fit is not None + assert fit["best_lag"] == 4 + assert fit["long_run_coef"] < 0 + + def test_per_lag_reconstruction_length_and_finite(self) -> None: + n, max_lag = 60, 6 + x = _aperiodic_rate_deltas(n) + beta = _hump_beta(max_lag, peak=2) + y = _y_from_lag_shape(x, beta, max_lag=max_lag, noise=0.001, seed=2) + fit = reg.fit_almon_dl(x, y, max_lag=max_lag, degree=2) + assert fit is not None + per_lag = fit["per_lag_coef"] + assert len(per_lag) == max_lag + 1 + assert all(math.isfinite(c) for c in per_lag) + + def test_hac_se_computed_for_every_lag(self) -> None: + n, max_lag = 60, 6 + x = _aperiodic_rate_deltas(n) + beta = _hump_beta(max_lag, peak=2) + y = _y_from_lag_shape(x, beta, max_lag=max_lag, noise=0.003, seed=3) + fit = reg.fit_almon_dl(x, y, max_lag=max_lag, degree=2) + assert fit is not None + hac_se = fit["hac_se"] + # One HAC SE per reconstructed lag coefficient, all finite and ≥0. + assert len(hac_se) == max_lag + 1 + assert all(math.isfinite(s) and s >= 0.0 for s in hac_se) + # Bandwidth follows the NW rule for this n. + assert fit["hac_bandwidth"] == reg.newey_west_bandwidth(fit["n"]) + + def test_degree_must_be_below_lag_count(self) -> None: + # degree ≥ max_lag+1 is not a constraint (degenerates to free lags) → refuse. + x = _aperiodic_rate_deltas(40) + beta = _hump_beta(6, peak=2) + y = _y_from_lag_shape(x, beta, max_lag=6, noise=0.0) + assert reg.fit_almon_dl(x, y, max_lag=6, degree=7) is None + + def test_thin_series_returns_none(self) -> None: + # Too few full rows to fit (< _MIN_FIT_OBS) → None, not a crash. + x = _aperiodic_rate_deltas(10) + beta = _hump_beta(6, peak=2) + y = _y_from_lag_shape(x, beta, max_lag=6, noise=0.0) + assert reg.fit_almon_dl(x, y, max_lag=6, degree=2) is None + + def test_zero_variance_y_returns_none(self) -> None: + x = _aperiodic_rate_deltas(50) + y: list[float | None] = [None] * 6 + [0.0] * 44 # constant → no variance + assert reg.fit_almon_dl(x, y, max_lag=6, degree=2) is None + + +# --------------------------------------------------------------------------- # +# build_fit_result — gate (mirror _elasticity_coef) → regression vs fallback +# --------------------------------------------------------------------------- # + + +_SEG: dict[str, str | None] = {"district": "Академический", "obj_class": None} + + +class TestBuildFitResult: + def test_gate_pass_emits_regression(self) -> None: + n, max_lag = 60, 6 + x = _aperiodic_rate_deltas(n, seed=13) + beta = _hump_beta(max_lag, peak=2) + y = _y_from_lag_shape(x, beta, max_lag=max_lag, noise=0.002, seed=0) + res = reg.build_fit_result(x, y, segment=_SEG, max_lag=max_lag, degree=2) + assert res.source == "regression" + assert res.best_lag_months == 2 + assert res.coef is not None and res.coef < 0 + assert res.x_pct is not None and res.x_pct < 0 # demand drops + assert res.r2 is not None and res.r2 > 0.8 + assert res.per_lag_coef is not None and len(res.per_lag_coef) == max_lag + 1 + assert res.hac_se is not None and len(res.hac_se) == max_lag + 1 + # Phrase carries the magnitude + peak lag. + assert "снижается" in res.phrase + assert f"{abs(round(res.x_pct, 1))}" in res.phrase + + def test_thin_n_degrades_to_fallback(self) -> None: + # Enough to fit, but n < _MIN_OBS (30) → gate fails on n → fallback. We keep + # the diagnostic numbers (per_lag/r2/n) but make no claim. + n, max_lag = 28, 6 # ~22 usable rows < 30 + x = _aperiodic_rate_deltas(n, seed=5) + beta = _hump_beta(max_lag, peak=2) + y = _y_from_lag_shape(x, beta, max_lag=max_lag, noise=0.001, seed=4) + res = reg.build_fit_result(x, y, segment=_SEG, max_lag=max_lag, degree=2) + assert res.source == "fallback" + assert res.n < reg._MIN_OBS + assert res.coef is None and res.x_pct is None and res.best_lag_months is None + assert res.phrase == reg._PHRASE_INSUFFICIENT + # Diagnostics retained (mirror _elasticity_coef returning r2/n in fallback). + assert res.per_lag_coef is not None + + def test_wrong_sign_degrades_to_fallback(self) -> None: + # POSITIVE long-run (rate ↑ → demand ↑) violates the gate sign → fallback, + # even with plenty of obs and a strong fit. + n, max_lag = 60, 6 + x = _aperiodic_rate_deltas(n, seed=13) + beta = -_hump_beta(max_lag, peak=2) # flip all signs → positive long-run + y = _y_from_lag_shape(x, beta, max_lag=max_lag, noise=0.002, seed=0) + res = reg.build_fit_result(x, y, segment=_SEG, max_lag=max_lag, degree=2) + assert res.source == "fallback" + assert res.coef is None + assert res.phrase == reg._PHRASE_INSUFFICIENT + + def test_weak_r2_degrades_to_fallback(self) -> None: + # Pure noise regressand (no rate link) at large n: a 3-param Almon basis + # cannot overfit ~114 noise points, so R² collapses well below 0.1 → the + # gate fails on R² (or sign) → fallback. (At small n a flexible basis can + # spuriously clear R²≥0.1 — which is exactly why the n≥30 gate + advisory + # status exist; here we use n=120 so the no-signal case is unambiguous.) + n, max_lag = 120, 6 + x = _aperiodic_rate_deltas(n, seed=13) + rng = np.random.default_rng(7) + y: list[float | None] = [None] * max_lag + [ + float(v) for v in rng.normal(0, 0.05, size=n - max_lag) + ] + res = reg.build_fit_result(x, y, segment=_SEG, max_lag=max_lag, degree=2) + assert res.source == "fallback" + assert res.coef is None + # Confirm it degraded specifically because the fit explains ~no variance. + assert res.r2 is not None and res.r2 < reg._MIN_R2 + + def test_empty_series_is_fallback_not_crash(self) -> None: + res = reg.build_fit_result([], [], segment=_SEG) + assert res.source == "fallback" + assert res.n == 0 + assert res.phrase == reg._PHRASE_INSUFFICIENT + + def test_as_dict_shape(self) -> None: + n, max_lag = 60, 6 + x = _aperiodic_rate_deltas(n, seed=13) + beta = _hump_beta(max_lag, peak=2) + y = _y_from_lag_shape(x, beta, max_lag=max_lag, noise=0.002, seed=0) + d = reg.build_fit_result(x, y, segment=_SEG, max_lag=max_lag, degree=2).as_dict() + for key in ( + "segment", + "best_lag_months", + "coef", + "x_pct", + "r2", + "n", + "per_lag_coef", + "hac_se", + "hac_bandwidth", + "almon_degree", + "source", + "phrase", + ): + assert key in d + assert d["source"] == "regression" + assert isinstance(d["per_lag_coef"], list) + + +# --------------------------------------------------------------------------- # +# _build_phrase +# --------------------------------------------------------------------------- # + + +class TestBuildPhrase: + def test_phrase_from_shape(self) -> None: + p = reg._build_phrase(x_pct=-3.2, best_lag=2, gated=True) + assert "3.2%" in p + assert "2 мес" in p + assert "снижается" in p + + def test_insufficient_when_not_gated(self) -> None: + assert reg._build_phrase(x_pct=-3.2, best_lag=2, gated=False) == reg._PHRASE_INSUFFICIENT + + def test_insufficient_when_none(self) -> None: + assert reg._build_phrase(x_pct=None, best_lag=2, gated=True) == reg._PHRASE_INSUFFICIENT + assert reg._build_phrase(x_pct=-3.2, best_lag=None, gated=True) == reg._PHRASE_INSUFFICIENT + + +# --------------------------------------------------------------------------- # +# compute_district_rate_regression — DB orchestrator (mocked) +# --------------------------------------------------------------------------- # + + +class _FakeMacro: + def __init__(self, month: dt.date, key_rate: float | None) -> None: + self.month = month + self.key_rate = key_rate + + +class _FakeSales: + def __init__(self, months: list[dt.date], units: list[int]) -> None: + self.months = months + self.units = units + + +def _months(n: int) -> list[dt.date]: + out: list[dt.date] = [] + y, m = 2021, 1 + for _ in range(n): + out.append(dt.date(y, m, 1)) + m += 1 + if m == 13: + m = 1 + y += 1 + return out + + +class TestComputeDistrictRateRegression: + def test_orchestrator_wires_macro_and_sales(self, monkeypatch: pytest.MonkeyPatch) -> None: + # Build a macro key_rate series whose Δ drives a lag-2 demand response, then + # confirm the orchestrator assembles X (Δrate) and Y (Δln units), aligns + # them, and recovers the injected lag via the pure fit. (The orchestrator + # uses the module-default max_lag=6 internally.) + n = 60 + months = _months(n) + # key_rate levels: integrate the aperiodic Δ so _delta() recovers them. + xdelta = _aperiodic_rate_deltas(n, seed=13) + levels: list[float] = [] + acc = 10.0 + for d in xdelta: + acc += d + levels.append(acc) + macro = [_FakeMacro(months[i], levels[i]) for i in range(n)] + + # Units carrying the lag-2 signal: ln(u_t) = ln(base) + Σ_{k≤t} β·Δrate[k-lag]. + beta_scalar = -0.05 + lag = 2 + ln_u = math.log(1000.0) + units: list[int] = [] + for t in range(n): + if t > 0: + src = xdelta[t - lag] if t - lag >= 0 else 0.0 + ln_u += beta_scalar * src + units.append(max(1, round(math.exp(ln_u)))) + sales = _FakeSales(months, units) + + monkeypatch.setattr(reg, "get_monthly_macro", lambda db, months_back: macro) + monkeypatch.setattr(reg, "build_sales_series", lambda db, spec, source, months_back: sales) + + res = reg.compute_district_rate_regression( + object(), # type: ignore[arg-type] + district="Академический", + months_back=n, + ) + assert res.segment["district"] == "Академический" + assert res.source == "regression" + # The single-lag injection at lag 2 → Almon shape peaks near lag 2. + assert res.best_lag_months in (1, 2, 3) + assert res.coef is not None and res.coef < 0 + assert res.n >= reg._MIN_OBS + + def test_orchestrator_graceful_on_empty(self, monkeypatch: pytest.MonkeyPatch) -> None: + monkeypatch.setattr(reg, "get_monthly_macro", lambda db, months_back: []) + monkeypatch.setattr( + reg, + "build_sales_series", + lambda db, spec, source, months_back: _FakeSales([], []), + ) + res = reg.compute_district_rate_regression( + object(), # type: ignore[arg-type] + district="Пустой", + ) + assert res.source == "fallback" + assert res.phrase == reg._PHRASE_INSUFFICIENT -- 2.45.3