Will the US fall into a recession in the next 12 months?

Walk-forward: each month is forecast using only recessions the NBER had announced by then. Before about 2010 most inputs come from later-revised data, so the early record is an optimistic ceiling.

1. Collect the data as it was known. Every series is stored with each version ever published (FRED/ALFRED), so the model only sees figures released by the decision date. Before a series' first archived vintage, its earliest vintage is used with a conservative release delay.
2. Transform each series into a change or spread (table below). Monthly changes are averaged over 3 months, because a single month is mostly noise. Daily and weekly series are averaged to months first.
   View the Python · recession.features.apply_transform 25 lines

   def apply_transform(x: pd.Series, tx: R.Tx, freq: str) -> pd.Series:
       ppy = PERIODS_PER_YEAR[freq]
       if tx is R.Tx.NONE:
           y = x
       elif tx is R.Tx.DIFF:
           y = x.diff()
       elif tx is R.Tx.DIFF2:
           y = x.diff().diff()
       elif tx is R.Tx.LOG_DIFF:
           y = 100 * _log(x).diff()
       elif tx is R.Tx.LOG_DIFF2:
           y = 100 * _log(x).diff().diff()
       elif tx is R.Tx.YOY:
           y = 100 * (x / x.shift(ppy) - 1)
       elif tx is R.Tx.PCT_6M_ANN:
           y = 100 * ((x / x.shift(max(ppy // 2, 1))) ** 2 - 1)
       elif tx is R.Tx.MAX_3Y_INCREASE:
           lp = 100 * _log(x)
           y = (lp - lp.shift(1).rolling(3 * ppy, min_periods=3 * ppy).max()).clip(lower=0)
       else:
           raise NotImplementedError(f"transform {tx.value} is not used by any Stage 3 feature")
       if tx in DIFF_LIKE and freq == "M":
           y = y.rolling(SMOOTH_MONTHS, min_periods=SMOOTH_MONTHS).mean()
       return y
   

3. Standardize each one against its own history up to that month only: z = (value − median) ÷ scale, where scale = interquartile range ÷ 1.349 (the standard deviation's equivalent for normal data), capped at ±4. Both are shown for every factor below.
   View the Python · recession.features.expanding_robust_z 16 lines

   def expanding_robust_z(paths: pd.DataFrame) -> pd.DataFrame:
       """
       z(T) from the median and IQR of the feature's own real-time path up to T.
   
       Event-type features (the Hamilton oil shock is zero most months; oil was
       price-controlled before 1973) have an IQR of zero. There the scale falls
       back to the expanding standard deviation, so a rare non-zero value still
       registers instead of dividing by zero.
       """
       ex = paths.expanding(MIN_HISTORY_MONTHS)
       med = ex.median()
       scale = (ex.quantile(0.75) - ex.quantile(0.25)) / 1.349
       scale = scale.where(scale > 0, ex.std())
       z = (paths - med) / scale.where(scale > 0)
       return z.replace([np.inf, -np.inf], np.nan).clip(-Z_CLIP, Z_CLIP)
   

4. Orient and average: multiply by −sign so negative always means recession-like, then average within each block. Indexes are left out before the date they were first published.
   View the Python · recession.features.block_scores 15 lines

   def block_scores(z: pd.DataFrame, specs: list[Spec], included: list[str]) -> tuple[pd.DataFrame, pd.DataFrame]:
       """Equal-weight mean of signed z within each block. Negative = recessionary."""
       sign = {s.key: s.sign for s in specs}
       z = mask_unpublished(z)
       signed = pd.DataFrame({k: -sign[k] * z[k] for k in included})
       block_of = {s.key: s.block for s in specs}
       names = {b.number: b.name for b in R.BLOCKS}
       scores, counts = {}, {}
       for b in sorted({block_of[k] for k in included}):
           cols = [k for k in included if block_of[k] == b]
           label = f"B{b:02d} {names[b]}"
           scores[label] = signed[cols].mean(axis=1)
           counts[label] = signed[cols].notna().sum(axis=1)
       return pd.DataFrame(scores), pd.DataFrame(counts)
   

5. Combine with a probit: probability = Φ(intercept + Σ weight × block score), where Φ is the standard normal distribution.
   View the Python · recession.model.Probit 75 lines

   class Probit:
       """
       Probit with an L2 penalty on slopes. Fisher-information covariance for intervals.
   
       sign_constrained=True forces every slope <= 0: inputs are block scores
       oriented so that NEGATIVE = recessionary (registry sign convention), so a
       weaker block may only RAISE the probability. The constraint encodes the
       registry's declared signs; it is not estimated from the data.
       """
   
       def __init__(self, lam: float, sign_constrained: bool = False):
           self.lam = lam
           self.sign_constrained = sign_constrained
           self.beta: np.ndarray | None = None
           self.cov: np.ndarray | None = None
   
       @staticmethod
       def _X(x: np.ndarray) -> np.ndarray:
           return np.column_stack([np.ones(len(x)), x])
   
       def fit(self, x: np.ndarray, y: np.ndarray) -> "Probit":
           X = self._X(x)
           pen = np.full(X.shape[1], self.lam)
           pen[0] = 0.0
   
           def nll(b):
               eta = X @ b
               ll = y * norm.logcdf(eta) + (1 - y) * norm.logcdf(-eta)
               # d/deta of the log-likelihood, computed stably
               g_pos = np.exp(norm.logpdf(eta) - norm.logcdf(eta))
               g_neg = np.exp(norm.logpdf(eta) - norm.logcdf(-eta))
               grad = -(X.T @ (y * g_pos - (1 - y) * g_neg)) + pen * b
               return -ll.sum() + 0.5 * (pen * b * b).sum(), grad
   
           b0 = np.zeros(X.shape[1])
           b0[0] = norm.ppf(np.clip(y.mean(), 0.01, 0.99))
           bounds = ([(None, None)] + [(None, 0.0)] * (X.shape[1] - 1)) if self.sign_constrained else None
           res = minimize(nll, b0, jac=True, method="L-BFGS-B", bounds=bounds)
           self.converged = bool(res.success)
           if not res.success:
               import warnings
               warnings.warn(f"Probit did not converge: {res.message}")
           self.beta = res.x
           eta = X @ self.beta
           w = norm.pdf(eta) ** 2 / np.clip(norm.cdf(eta) * norm.cdf(-eta), 1e-12, None)
           H = X.T @ (X * w[:, None]) + np.diag(pen)
           self.H = H
           self.cov = np.linalg.pinv(H)
           return self
   
       def predict(self, x: np.ndarray) -> np.ndarray:
           return norm.cdf(self._X(x) @ self.beta)
   
       def interval(self, x: np.ndarray, q=(0.1, 0.9), draws: int = 4000, seed: int = 0) -> np.ndarray:
           """
           Probability quantiles from parameter uncertainty (Laplace approximation).
   
           With the sign constraint, slopes the fit holds at exactly zero (binding) are
           kept at zero, and the uncertainty of the others is the curvature of the fit in
           the free directions only. Drawing every slope from the full covariance (as
           before 5 Oct 2026) gave half the zero-weight blocks positive weights, which the
           model forbids, and let blocks with no weight move the range (external review).
           """
           rng = np.random.default_rng(seed)
           if not self.sign_constrained or getattr(self, "H", None) is None:
               B = rng.multivariate_normal(self.beta, self.cov, size=draws)
           else:
               free = np.r_[True, self.beta[1:] < -1e-8]          # intercept and the non-zero slopes
               cov_free = np.linalg.pinv(self.H[np.ix_(free, free)])
               B = np.zeros((draws, len(self.beta)))
               B[:, free] = rng.multivariate_normal(self.beta[free], cov_free, size=draws)
               B[:, 1:] = np.minimum(B[:, 1:], 0.0)               # a weaker block may only raise the probability
           p = norm.cdf(self._X(x) @ B.T)
           return np.quantile(p, q, axis=1).T
   

   View the Python · recession.live.forecast 91 lines

   def forecast(scores: Optional[pd.DataFrame] = None, paths: Optional[pd.DataFrame] = None,
                counts: Optional[pd.DataFrame] = None, decision: Optional[str] = None,
                as_of: Optional[str] = None) -> dict:
       """
       decision="YYYY-MM": the official month-end forecast for that month.
       as_of="YYYY-MM-DD": an in-month update using everything published by that day
       (requires the Stage 3 paths to include that date: run_stage3(extra_date=...)).
       """
       fd = features.FEATURE_DIR
       scores = pd.read_parquet(fd / "block_scores.parquet") if scores is None else scores
       paths = pd.read_parquet(fd / "feature_paths.parquet") if paths is None else paths
       counts = pd.read_parquet(fd / "block_counts.parquet") if counts is None else counts
       panel = _panel(scores, paths)
       blocks = model.block_columns(panel)
       present = panel[blocks].notna().sum(axis=1) >= model.MIN_BLOCKS_PRESENT
       if as_of is not None:
           T = pd.Timestamp(as_of).normalize()                     # in-month (event) update
       elif decision:
           T = pd.Timestamp(decision) + pd.offsets.MonthEnd(0)     # official month-end
       else:
           T = panel.index[present & panel.index.is_month_end][-1]
       earlier = panel.index[present & (panel.index < T)]
       prev = earlier[-1]
       # An update for a given day is in-month even on the last day of a month; only decision= (or the default)
       # makes the official month-end forecast. (A 30 Sep update used to become September's official entry.)
       kind = "in-month" if as_of is not None else "month-end"
       Xb = panel[blocks].fillna(0.0)
   
       tr = model._training_rows(panel, T) & present
       y = panel.loc[tr, "y"].values
       v2 = Probit(model.RIDGE_LAMBDA, sign_constrained=True).fit(Xb.loc[tr].values, y)
       ts_rows = tr & panel["term_spread"].notna()
       ts = Probit(1e-6).fit(panel.loc[ts_rows, ["term_spread"]].values, panel.loc[ts_rows, "y"].values)
       te_rows = ts_rows & panel["ebp"].notna()
       tsebp = Probit(1e-6).fit(panel.loc[te_rows, ["term_spread", "ebp"]].values, panel.loc[te_rows, "y"].values)
   
       x, xp = Xb.loc[[T]].values, Xb.loc[[prev]].values
       p, p_prev = float(v2.predict(x)[0]), float(v2.predict(xp)[0])
       lo, hi = (float(v) for v in v2.interval(x)[0])
       beta = v2.beta[1:]
   
       def bench(fit, cols, when):
           row = panel.loc[[when], cols]
           return float(fit.predict(row.values)[0]) if row.notna().all(axis=None) else None
   
       ind = model.indicator_benchmarks(pd.DatetimeIndex([T])).iloc[0]
       block_rows = {}
       for j, b in enumerate(blocks):
           now, before = float(Xb.loc[T, b]), float(Xb.loc[prev, b])
           recent = counts[b].loc[:T].iloc[-25:-1]
           block_rows[b] = {
               "score": round(now, 3), "score_prev": round(before, 3),
               "weight": round(float(beta[j]), 6),
               "contribution": round(float(beta[j] * now), 4),          # probit index units; + = more risk
               "change_contribution": round(float(beta[j] * (now - before)), 4),
               "factors": int(counts.loc[T, b]) if T in counts.index and pd.notna(counts.loc[T, b]) else 0,
               "usual_factors": int(recent.max()) if len(recent) else 0,
           }
   
       out = {
           "decision_month": f"{T:%Y-%m}",
           "kind": kind,
           "as_of": f"{T:%Y-%m-%d}",
           "compared_with": f"{prev:%Y-%m-%d}",
           "probability": round(p, 4), "probability_exact": p,
           "intercept": round(float(v2.beta[0]), 6),
           "ridge_lambda": model.RIDGE_LAMBDA,
           "interval_80": [round(lo, 4), round(hi, 4)],
           "probability_prev_month_same_model": round(p_prev, 4),
           "alarm_threshold": DECISIONS.alarm_threshold,
           "alarm_model": p >= DECISIONS.alarm_threshold,   # the 12-block model's own reading
           "benchmarks": {
               "Term spread + EBP probit": bench(tsebp, ["term_spread", "ebp"], T),
               "Term spread probit": bench(ts, ["term_spread"], T),
               "Fed Board TS+EBP probability (published)": _num(ind["Fed Board TS+EBP probability (published)"]),
               "Chauvet-Piger probability": _num(ind["Chauvet-Piger probability"]),
               "Sahm rule (real-time)": _num(ind["Sahm rule (real-time)"]),
           },
           "blocks": block_rows,
           "features_included": None,   # filled by the caller from the Stage 3 report
           "train_months": int(tr.sum()),
           "base_rate": round(float(y.mean()), 4),   # what "no skill" would have forecast (K3)
           "recessions_in_training": int(panel.loc[tr & (panel["y"] == 1), "peak"].nunique()),
       }
       out["flags"] = flags(out)
       out["kill_criteria"] = kill_status(out)
       # The alarm follows the HEADLINE (the benchmark while K3/K4/K5 hold), not the demoted model.
       out["alarm"] = (headline(out)["probability"] or 0) >= DECISIONS.alarm_threshold
       out["explanation"] = explain(out)
       return out