3.2 Fama-French Factors and Valuation Data

# Core libraries
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from datetime import datetime
import warnings
warnings.filterwarnings('ignore')

# Data fetching
import pandas_datareader as pdr
from fredapi import Fred

# Statistical modeling
from sklearn.decomposition import PCA, SparsePCA
from sklearn.preprocessing import StandardScaler
from sklearn.covariance import LedoitWolf
from hmmlearn import hmm
from statsmodels.tsa.api import VAR
from statsmodels.tsa.stattools import grangercausalitytests, adfuller, kpss
from statsmodels.stats.diagnostic import acorr_ljungbox
from scipy import stats
from scipy.stats import jarque_bera, shapiro

# Plotting
plt.style.use('seaborn-v0_8-darkgrid')
sns.set_palette("husl")
plt.rcParams['figure.figsize'] = (14, 8)
plt.rcParams['font.size'] = 11

print("✅ All libraries imported successfully")
print(f"Analysis date: {datetime.now().strftime('%Y-%m-%d')}")
✅ All libraries imported successfully
Analysis date: 2026-01-26
start_date = "1990-01-01"
end_date   = "2025-12-31"

# --- FRED series codes ---
FRED_SERIES = {
    "M2":   "M2SL",      # M2 money stock (monthly, SA)
    "FED_BS":"WALCL",     # Fed balance sheet total assets (weekly)
    "FFR": "DFF",  
    "Real_Rate": "REAINTRATREARAT1YE",
    "T3M":   "DTB3",     # 3-Month T-Bill rate (monthly)
    "T10Y":  "DGS10",     # 10-Year Treasury yield (daily -> resample monthly)
    "BAA":    "DBAA",       # Moody's Baa corporate yield (monthly)
    "AAA":    "DAAA",       # Moody's Aaa corporate yield (monthly)
    "CPI":    "CPIAUCSL",  # CPI index (monthly, SA)
    "CORE_CPI": "CPILFESL",    # Core CPI
    "GDP":    "GDP",       # Nominal GDP (quarterly, SAAR)
    
    'IORB': 'IORB',                    # Interest on Reserve Balances
    'RRP': 'RRPONTSYD',                # Overnight Reverse Repo
    'TGA': 'WTREGEN',                  # Treasury General Account
    
    'Bank_Reserves': 'TOTRESNS',       # Total Reserves
    'CI_Loans': 'TOTBKCR',             # Bank Credit, All Commercial Banks
    'Consumer_Credit': 'TOTALSL',      # Consumer Credit Outstanding
    'Mortgage_Rate': 'MORTGAGE30US',   # 30Y Mortgage Rate
    'Prime_Rate': 'DPRIME',            # Bank Prime Loan Rate
    'STLFSI': 'STLFSI',                # St. Louis Fed Financial Stress
    'Dollar_Index': 'DTWEXBGS',        # # Global Liquidity: Trade-Weighted Dollar
    'Household_Debt_Service': 'TDSP',  # Household Liquidity: Household Debt Service Ratio
}

# For equity factors: Fama-French via pandas_datareader (famafrench)
FF_FACTORS_DATASET = "F-F_Research_Data_5_Factors_2x3"

def download_fred_series(series_dict, start, end):
    """
    Download FRED series into a single monthly DataFrame.

    Parameters
    ----------
    series_dict : dict
        Mapping logical_name -> FRED code.
    """
    dfs = []
    for name, code in series_dict.items():
        print(f"Downloading {name} ({code}) from FRED...")
        s = pdr.DataReader(code, "fred", start, end)
        s = s.rename(columns={code: name})
        dfs.append(s)

    df = pd.concat(dfs, axis=1)
    # Ensure monthly freq by end-of-month sampling
    df = df.resample("M").last()
    return df
macro_data = download_fred_series(FRED_SERIES, start_date, end_date)

# 1. Start from the original GDP values only (drop NaNs)
gdp_series = macro_data["GDP"].dropna()
# 2. Collapse to unique quarter-end values
#    - If GDP is timestamped at quarter *start* (e.g. 2025-04-01),
#      this will move it to the quarter end (2025-06-30) and give unique labels.
gdp_q = gdp_series.resample("Q").last()   # quarterly, at quarter-end (e.g. 2025-03-31, 2025-06-30, ...)

# 3. Downsample to monthly and forward-fill within the quarter
gdp_m = gdp_q.resample("M").ffill()

# 4. Assign back into macro_raw, aligning on the monthly index
macro_data["GDP"] = gdp_m

macro_data.tail()
Downloading M2 (M2SL) from FRED...
Downloading FED_BS (WALCL) from FRED...
Downloading FFR (DFF) from FRED...
Downloading Real_Rate (REAINTRATREARAT1YE) from FRED...
Downloading T3M (DTB3) from FRED...
Downloading T10Y (DGS10) from FRED...
Downloading BAA (DBAA) from FRED...
Downloading AAA (DAAA) from FRED...
Downloading CPI (CPIAUCSL) from FRED...
Downloading CORE_CPI (CPILFESL) from FRED...
Downloading GDP (GDP) from FRED...
Downloading IORB (IORB) from FRED...
Downloading RRP (RRPONTSYD) from FRED...
Downloading TGA (WTREGEN) from FRED...
Downloading Bank_Reserves (TOTRESNS) from FRED...
Downloading CI_Loans (TOTBKCR) from FRED...
Downloading Consumer_Credit (TOTALSL) from FRED...
Downloading Mortgage_Rate (MORTGAGE30US) from FRED...
Downloading Prime_Rate (DPRIME) from FRED...
Downloading STLFSI (STLFSI) from FRED...
Downloading Dollar_Index (DTWEXBGS) from FRED...
Downloading Household_Debt_Service (TDSP) from FRED...
M2 FED_BS FFR Real_Rate T3M T10Y BAA AAA CPI CORE_CPI ... RRP TGA Bank_Reserves CI_Loans Consumer_Credit Mortgage_Rate Prime_Rate STLFSI Dollar_Index Household_Debt_Service
DATE
2025-08-31 22108.3 6603384.0 4.33 1.523005 4.05 4.23 6.03 5.42 323.364 329.793 ... 77.898 589998.0 3281.9 18720.1524 5059896.38 6.56 7.50 NaN 120.6028 NaN
2025-09-30 22212.4 6608395.0 4.09 1.149868 3.86 4.16 5.83 5.22 324.368 330.542 ... 49.071 804856.0 3068.1 18759.8099 5071365.99 6.30 7.25 NaN 120.5624 NaN
2025-10-31 22298.0 6587034.0 3.86 0.973381 3.73 4.11 5.80 5.22 NaN NaN ... 51.802 957990.0 2944.9 18820.7752 5080601.87 6.17 7.00 NaN 121.7715 NaN
2025-11-30 22322.4 6552419.0 3.89 1.001027 3.73 4.02 5.80 5.18 325.031 331.068 ... 7.561 903394.0 2879.3 18925.1502 5084831.24 6.23 7.00 NaN 121.4288 NaN
2025-12-31 NaN 6640618.0 3.64 -0.135174 3.57 4.18 5.90 5.35 326.030 331.860 ... 105.993 837306.0 NaN 18934.1126 NaN 6.15 6.75 NaN 120.1166 NaN

5 rows × 22 columns

def download_ff_factors(start, end):
    """
    Download monthly Fama-French 5 factors (2x3) using pandas_datareader.
    """
    print("Downloading Fama-French factors (famafrench)...")
    ff_raw = pdr.DataReader(
        name=FF_FACTORS_DATASET,
        data_source="famafrench",
        start=start,
        end=end
    )[0]  # table [0] contains the data

    ff = (ff_raw
          .divide(100)  # convert from % to decimal
          .reset_index(names="date")
          .assign(date=lambda x: pd.to_datetime(x["date"].astype(str)))
          .rename(str.lower, axis="columns")
          .rename(columns={"mkt-rf": "mkt_excess"}))

    ff = ff.set_index("date")
    return ff

ff_factors = download_ff_factors(start_date, end_date)
ff_adj = ff_factors.copy()
ff_adj.index = ff_adj.index.to_period("M").to_timestamp("M")
ff_adj.tail()
Downloading Fama-French factors (famafrench)...
mkt_excess smb hml rmw cma rf
date
2025-07-31 0.0198 -0.0015 -0.0127 -0.0029 -0.0207 0.0034
2025-08-31 0.0185 0.0488 0.0442 -0.0067 0.0208 0.0038
2025-09-30 0.0339 -0.0218 -0.0105 -0.0203 -0.0222 0.0033
2025-10-31 0.0196 -0.0131 -0.0309 -0.0522 -0.0403 0.0037
2025-11-30 -0.0013 0.0147 0.0376 0.0143 0.0068 0.0030

3.3 Shiller CAPE Data

https://shillerdata.com/

import pandas as pd
import numpy as np

def parse_yyyy_mm_to_date(s: str):
    """
    Convert 'YYYY.MM' to a proper datetime.
    Handles the special case where October appears as YYYY.1 instead of YYYY.10.
    """
    s = str(s).strip()
    parts = s.split(".")
    if len(parts) != 2:
        raise ValueError(f"Unexpected date format: {s}")
    
    year = int(parts[0])
    mm = parts[1]

    # Fix: "1" should be "10" (October)
    if mm == "1":
        month = 10
    else:
        month = int(mm)

    return pd.Timestamp(year=year, month=month, day=1)

# --- Load CAPE CSV ---
sp500_cape_raw = pd.read_csv("data/sp500-cape-1990-2025.csv")

# Clean + parse
sp500_cape_raw["date"] = sp500_cape_raw["Date"].apply(parse_yyyy_mm_to_date)

sp500_cape_raw = (
    sp500_cape_raw[["date", "SP500", "CAPE"]]
    .set_index("date")
    .sort_index()
)

# Convert to month-end aligned CAPE series
sp500_cape_m = sp500_cape_raw.resample("M").last()
sp500_cape_m.rename(columns={"CAPE": "CAPE", "SP500": "P"}, inplace=True)

sp500_cape_m.tail()
P CAPE
date
2025-08-31 6408.95 37.85
2025-09-30 6584.02 38.58
2025-10-31 6735.69 39.21
2025-11-30 6740.89 39.12
2025-12-31 6812.63 39.42
import pandas_datareader as pdr

# Fetch the VIXCLS series from FRED
vix_data = pdr.get_data_fred('VIXCLS', start=start_date, end=end_date)

# Display the results
print(vix.tail())
            VIXCLS
DATE              
2025-12-25     NaN
2025-12-26   13.60
2025-12-29   14.20
2025-12-30   14.33
2025-12-31   14.95
# Combine all data
combined_df = pd.concat([
    macro_data,
    ff_adj,
    sp500_cape_m,
    vix_data.rename(columns={"VIXCLS": "VIX"})
], axis=1)

# Drop rows with too many missing values
combined_df = combined_df.dropna(thresh=len(combined_df.columns) * 0.7)
combined_df.tail()
M2 FED_BS FFR Real_Rate T3M T10Y BAA AAA CPI CORE_CPI ... Household_Debt_Service mkt_excess smb hml rmw cma rf P CAPE VIX
2025-07-31 22028.7 6642578.0 4.33 0.838904 4.24 4.37 6.04 5.41 322.132 328.656 ... 11.256338 0.0198 -0.0015 -0.0127 -0.0029 -0.0207 0.0034 6296.50 37.47 16.72
2025-08-31 22108.3 6603384.0 4.33 1.523005 4.05 4.23 6.03 5.42 323.364 329.793 ... NaN 0.0185 0.0488 0.0442 -0.0067 0.0208 0.0038 6408.95 37.85 NaN
2025-09-30 22212.4 6608395.0 4.09 1.149868 3.86 4.16 5.83 5.22 324.368 330.542 ... NaN 0.0339 -0.0218 -0.0105 -0.0203 -0.0222 0.0033 6584.02 38.58 16.28
2025-10-31 22298.0 6587034.0 3.86 0.973381 3.73 4.11 5.80 5.22 NaN NaN ... NaN 0.0196 -0.0131 -0.0309 -0.0522 -0.0403 0.0037 6735.69 39.21 17.44
2025-11-30 22322.4 6552419.0 3.89 1.001027 3.73 4.02 5.80 5.18 325.031 331.068 ... NaN -0.0013 0.0147 0.0376 0.0143 0.0068 0.0030 6740.89 39.12 NaN

5 rows × 31 columns

4. Feature Engineering

4.1 Construct Liquidity Variables

def engineer_liquidity_features(df):
    """
    Create liquidity-related features from raw macro data
    """
    df = df.copy()
    
    # 1. Growth rates
    df['M2_growth'] = df['M2'].pct_change(12) * 100  # YoY % change
    df['FED_BS_growth'] = df['FED_BS'].pct_change(12) * 100
    
    # 2. Term spread
    df['Term_Spread'] = df['T10Y'] - df['FFR']
    
    # 3. Credit spread
    df['Credit_Spread'] = df['BAA'] - df['AAA']
    
    # 4. Real interest rate (ex-post using CPI)
    df['CPI_inflation'] = df['CPI'].pct_change(12) * 100  # YoY inflation
    #df['Real_Rate'] = df['FFR'] - df['CPI_inflation']
    
    # 5. Valuation metrics
    df['log_CAPE'] = np.log(df['CAPE'])
    df['CAPE_zscore'] = (df['CAPE'] - df['CAPE'].rolling(60).mean()) / df['CAPE'].rolling(60).std()
    
    # 6. Risk-adjusted term spread
    df['Risk_Adj_Spread'] = df['Term_Spread'] / (df['VIX'] / 100)

    #7 Net liquidity
    df['Net_Liq'] = df['FED_BS'] - (df['TGA'] + df['RRP'])
    
    return df

# Apply feature engineering
df_features = engineer_liquidity_features(combined_df)

print("✅ Feature engineering complete")
print(f"\nNew features created:")
new_cols = [c for c in df_features.columns if c not in combined_df.columns]
for col in new_cols:
    print(f"  - {col}")

df_features.tail(5)
✅ Feature engineering complete

New features created:
  - M2_growth
  - FED_BS_growth
  - Term_Spread
  - Credit_Spread
  - CPI_inflation
  - log_CAPE
  - CAPE_zscore
  - Risk_Adj_Spread
  - Net_Liq
M2 FED_BS FFR Real_Rate T3M T10Y BAA AAA CPI CORE_CPI ... VIX M2_growth FED_BS_growth Term_Spread Credit_Spread CPI_inflation log_CAPE CAPE_zscore Risk_Adj_Spread Net_Liq
2025-07-31 22028.7 6642578.0 4.33 0.838904 4.24 4.37 6.04 5.41 322.132 328.656 ... 16.72 4.480153 -7.464249 0.04 0.63 2.731801 3.623541 1.242869 0.239234 6271856.555
2025-08-31 22108.3 6603384.0 4.33 1.523005 4.05 4.23 6.03 5.42 323.364 329.793 ... NaN 4.427283 -7.298001 -0.10 0.61 2.939220 3.633631 1.306865 NaN 6013308.102
2025-09-30 22212.4 6608395.0 4.09 1.149868 3.86 4.16 5.83 5.22 324.368 330.542 ... 16.28 4.492553 -6.661865 0.07 0.61 3.022700 3.652734 1.461587 0.429975 5803489.929
2025-10-31 22298.0 6587034.0 3.86 0.973381 3.73 4.11 5.80 5.22 NaN NaN ... 17.44 4.645651 -6.080511 0.25 0.58 2.789925 3.668932 1.575675 1.433486 5628992.198
2025-11-30 22322.4 6552419.0 3.89 1.001027 3.73 4.02 5.80 5.18 325.031 331.068 ... NaN 4.271300 -5.108093 0.13 0.62 2.711969 3.666634 1.490940 NaN 5649017.439

5 rows × 40 columns

4.2 Construct Valuation Spread Measure

We define valuation spread as the return differential between cheap (high B/M) and expensive (low B/M) stocks. This is proxied by the HML (High Minus Low) factor from Fama-French.

def construct_valuation_spread(df):
    """
    Construct valuation spread measure
    
    Primary measure: HML (Fama-French) represents value premium
    Secondary: CAPE z-score as market-level valuation
    """
    df = df.copy()
    
    # Cumulative HML performance (value vs growth)
    df['HML_cumret'] = (1 + df['hml'] / 100).cumprod()
    
    # Rolling 12-month HML performance
    df['HML_12m'] = df['hml'].rolling(12).sum()
    
    # Define valuation spread as:
    # Positive = cheap stocks outperforming (value regime)
    # Negative = expensive stocks outperforming (growth regime)
    df['V_spread'] = df['HML_12m']  # 12-month rolling HML
    
    # Standardize
    df['V_spread_z'] = (df['V_spread'] - df['V_spread'].mean()) / df['V_spread'].std()
    
    # Also use CAPE z-score as alternative valuation measure
    df['V_level'] = df['CAPE_zscore']
    
    return df

df_features = construct_valuation_spread(df_features)

print("✅ Valuation spread measures constructed")
print(f"\nV_spread (HML 12m) stats:")
print(df_features['V_spread'].describe())
df_features.tail(5)
✅ Valuation spread measures constructed

V_spread (HML 12m) stats:
count    420.000000
mean       0.018185
std        0.141636
min       -0.417500
25%       -0.078100
50%        0.016600
75%        0.095125
max        0.579800
Name: V_spread, dtype: float64
M2 FED_BS FFR Real_Rate T3M T10Y BAA AAA CPI CORE_CPI ... CPI_inflation log_CAPE CAPE_zscore Risk_Adj_Spread Net_Liq HML_cumret HML_12m V_spread V_spread_z V_level
2025-07-31 22028.7 6642578.0 4.33 0.838904 4.24 4.37 6.04 5.41 322.132 328.656 ... 2.731801 3.623541 1.242869 0.239234 6271856.555 1.005492 -0.0558 -0.0558 -0.522355 1.242869
2025-08-31 22108.3 6603384.0 4.33 1.523005 4.05 4.23 6.03 5.42 323.364 329.793 ... 2.939220 3.633631 1.306865 NaN 6013308.102 1.005936 -0.0006 -0.0006 -0.132625 1.306865
2025-09-30 22212.4 6608395.0 4.09 1.149868 3.86 4.16 5.83 5.22 324.368 330.542 ... 3.022700 3.652734 1.461587 0.429975 5803489.929 1.005831 0.0166 0.0166 -0.011187 1.461587
2025-10-31 22298.0 6587034.0 3.86 0.973381 3.73 4.11 5.80 5.22 NaN NaN ... 2.789925 3.668932 1.575675 1.433486 5628992.198 1.005520 -0.0229 -0.0229 -0.290070 1.575675
2025-11-30 22322.4 6552419.0 3.89 1.001027 3.73 4.02 5.80 5.18 325.031 331.068 ... 2.711969 3.666634 1.490940 NaN 5649017.439 1.005898 0.0132 0.0132 -0.035192 1.490940

5 rows × 45 columns

5. Liquidity Index Construction via Sparse PCA

5.1 Methodology

We construct a composite liquidity index using Sparse PCA (Zou et al. 2006, Witten et al. 2009) which:

  1. Maximizes explained variance like standard PCA
  2. Enforces sparsity via L1 penalty for interpretability
  3. Selects only the most relevant liquidity indicators

Liquidity proxies: - Fed balance sheet growth (QE intensity) - Real interest rates (opportunity cost) - Term spread (yield curve steepness) - Credit spreads (credit conditions) - VIX (risk appetite)

All variables are sign-adjusted so higher values = easier liquidity conditions.

df_features.tail(1)
M2 FED_BS FFR Real_Rate T3M T10Y BAA AAA CPI CORE_CPI ... Credit_Spread CPI_inflation log_CAPE CAPE_zscore Risk_Adj_Spread HML_cumret HML_12m V_spread V_spread_z V_level
2025-11-30 22322.4 6552419.0 3.89 1.001027 3.73 4.02 5.8 5.18 325.031 331.068 ... 0.62 2.711969 3.666634 1.49094 NaN 1.005898 0.0132 0.0132 -0.035192 1.49094

1 rows × 33 columns

def build_liquidity_index(df, method='sparse', n_components=1, alpha=1.0):
    """
    Build liquidity index using PCA or Sparse PCA
    
    Parameters:
    -----------
    df : DataFrame
        Input data
    method : str
        'pca' or 'sparse'
    n_components : int
        Number of components
    alpha : float
        Sparsity penalty (only for sparse PCA)
    """
    # Select liquidity variables
    liq_vars = [
        'M2_growth',
        'FED_BS_growth',
        'Term_Spread',
        'Real_Rate',
        'Credit_Spread',
        'VIX'
    ]
    
    # Extract data and drop NaN
    X = df[liq_vars].copy()
    X = X.dropna()
    
    # Sign adjustments: higher = easier liquidity
    X['Real_Rate'] = -X['Real_Rate']  # Lower real rates = easier
    X['Credit_Spread'] = -X['Credit_Spread']  # Tighter spreads = easier
    X['VIX'] = -X['VIX']  # Lower VIX = easier
    
    # Standardize
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)
    X_scaled_df = pd.DataFrame(X_scaled, index=X.index, columns=X.columns)
    
    # Fit PCA/Sparse PCA
    if method == 'sparse':
        model = SparsePCA(n_components=n_components, alpha=alpha, random_state=42)
    else:
        model = PCA(n_components=n_components)
    
    # Fit and transform
    L = model.fit_transform(X_scaled)
    
    # Create series
    L_series = pd.Series(L[:, 0], index=X.index, name='L')
    
    # Get loadings
    loadings = pd.DataFrame(
        model.components_.T,
        index=X.columns,
        columns=['PC1']
    )
    
    # Calculate explained variance
    if method == 'pca':
        explained_var = model.explained_variance_ratio_[0]
    else:
        # For sparse PCA, compute manually
        var_explained = np.var(L[:, 0]) / np.sum(np.var(X_scaled, axis=0))
        explained_var = var_explained
    
    print(f"\n{'='*70}")
    print(f"Liquidity Index Construction ({method.upper()})")
    print(f"{'='*70}")
    print(f"\nVariance explained: {explained_var:.2%}")
    print(f"\nLoadings:")
    print(loadings.sort_values('PC1', ascending=False))
    print(f"\nLiquidity index range: {L_series.min():.2f} to {L_series.max():.2f}")
    print(f"Mean: {L_series.mean():.2f}, Std: {L_series.std():.2f}")
    
    return L_series, loadings, model, X_scaled_df

# Build liquidity index
L_index, loadings, pca_model, scaled_features = build_liquidity_index(
    df_features,
    method='sparse',
    alpha=0.5  # Moderate sparsity
)

======================================================================
Liquidity Index Construction (SPARSE)
======================================================================

Variance explained: 44.98%

Loadings:
                    PC1
FED_BS_growth  0.525122
M2_growth      0.379190
Term_Spread    0.370752
Real_Rate      0.244985
Credit_Spread -0.414407
VIX           -0.459623

Liquidity index range: -3.15 to 7.11
Mean: 0.00, Std: 1.65

Liquidity Index Augmentation: Two-layer liquidity: Flow vs Stock / Excess

Add level / excess variables

Examples (all can be pulled from FRED):

Log level vs pre-2008 trend

Let \(m_t = \log M2_t\). Fit a linear trend on a pre-QE baseline (say 1985–2007):

\(m_t \approx a + bt \quad (t \le 2007)\)

Then define excess money stock:

\(EM_t = m_t - (a + bt)\)

After 2008, \(EM_t\) becomes increasingly positive if M2 grows above its historical trend.

Similarly for the Fed Balance Sheet, let \(b_t = \log(\text{FedBal}_t)\):

\(EB_t = b_t - (\alpha + \beta t)\) (pre-2007 trend)


Excess liquidity versus real economy

Use real GDP series (e.g., GDPC1) and define:

\(EL_t = \log M2_t - \log GDP_t\)

Or measure multi-year excess growth:

\(EL_t(3y) = \log M2_t - \log M2_{t-36} - (\log GDP_t - \log GDP_{t-36})\)


ZIRP / negative real-rate indicator

\(D_t^{ZIRP} = 1(r_t^{real} < 0)\)

From 2009–2015, this should be mostly 1.


Build an augmented feature vector

Instead of only:

\(x_t = \{ \Delta \log M2_t,\; \Delta \log FedBal_t,\; TS_t,\; r_t^{real},\; CS_t \}\)

Let’s expand to:

\(x_t^{aug} = \{ \Delta \log M2_t,\; \Delta \log FedBal_t,\; TS_t,\; r_t^{real},\; CS_t,\; EM_t,\; EB_t,\; EL_t,\; EL_t(3y),\; D_t^{ZIRP} \}\)

Then standardize and perform PCA / SparsePCA on this.

  • PC1 (or a combination) loading heavily on \(EM\), \(EB\), \(EL\), \(D_t^{ZIRP}\) → structural easy-money regime
  • PC2 (or your current PC1) loading on short-horizon changes → flow / shock liquidity
df_features.columns
Index(['M2', 'FED_BS', 'FFR', 'Real_Rate', 'T3M', 'T10Y', 'BAA', 'AAA', 'CPI',
       'CORE_CPI', 'GDP', 'IORB', 'RRP', 'TGA', 'Bank_Reserves', 'CI_Loans',
       'Consumer_Credit', 'Mortgage_Rate', 'Prime_Rate', 'STLFSI',
       'Dollar_Index', 'Household_Debt_Service', 'mkt_excess', 'smb', 'hml',
       'rmw', 'cma', 'rf', 'P', 'CAPE', 'VIX', 'M2_growth', 'FED_BS_growth',
       'Term_Spread', 'Credit_Spread', 'CPI_inflation', 'log_CAPE',
       'CAPE_zscore', 'Risk_Adj_Spread', 'Net_Liq', 'HML_cumret', 'HML_12m',
       'V_spread', 'V_spread_z', 'V_level'],
      dtype='object')
def augment_liquidity_index(df_features):
    df = df_features.copy()

    # Flow proxies
    df["dlog_M2"]      = np.log(df["M2"]).diff()
    df["dlog_FED_BS"]  = np.log(df["FED_BS"]).diff()

    # Levels
    df["log_M2"]      = np.log(df["M2"])
    df["log_FED_BS"] = np.log(df["FED_BS"])

    # Use data up to 2007-12-31 to fit trends
    pre = df.loc[: "2007-12-31"].copy()

    # --- Trend for M2 ---
    pre_m2 = pre["log_M2"].dropna()
    t_m2   = np.arange(len(pre_m2))
    coefs_M2 = np.polyfit(t_m2, pre_m2.values, deg=1)

    t_full = np.arange(len(df))
    trend_M2 = np.polyval(coefs_M2, t_full)
    df["EM"] = df["log_M2"] - trend_M2

    # --- Trend for Fed balance sheet ---
    pre_fb = pre["log_FED_BS"].dropna()
    t_fb   = np.arange(len(pre_fb))
    coefs_FB = np.polyfit(t_fb, pre_fb.values, deg=1)

    trend_FB = np.polyval(coefs_FB, t_full)
    df["EB"] = df["log_FED_BS"] - trend_FB

    # 3. Excess liquidity vs GDP over 3y
    df["log_GDP"] = np.log(df["GDP"])
    df["EL_3y"] = (df["log_M2"] - df["log_M2"].shift(36)) - \
                  (df["log_GDP"] - df["log_GDP"].shift(36))

    return df

df_features_aug = augment_liquidity_index(df_features)
df_features_aug.tail(1)
M2 FED_BS FFR Real_Rate T3M T10Y BAA AAA CPI CORE_CPI ... V_spread_z V_level dlog_M2 dlog_FED_BS log_M2 log_FED_BS EM EB log_GDP EL_3y
2025-11-30 22322.4 6552419.0 3.89 1.001027 3.73 4.02 5.8 5.18 325.031 331.068 ... -0.035192 1.49094 0.001094 -0.005269 10.013346 15.695345 0.19326 0.685233 NaN NaN

1 rows × 53 columns

def build_liquidity_index_aug(df, method='sparse', n_components=1, alpha=1.0):
    """
    Build liquidity index using PCA or Sparse PCA
    
    Parameters:
    -----------
    df : DataFrame
        Input data
    method : str
        'pca' or 'sparse'
    n_components : int
        Number of components
    alpha : float
        Sparsity penalty (only for sparse PCA)
    """
    # Select liquidity variables
    liq_vars = [
        'M2_growth',
        'FED_BS_growth',
        'EM',
        'EB',
        'EL_3y',
        'Term_Spread',
        'Real_Rate',
        'Credit_Spread',
        'VIX'
    ]
    
    # Extract data and drop NaN
    X = df[liq_vars].copy()
    X = X.dropna()
    
    # Sign adjustments: higher = easier liquidity
    X['Real_Rate'] = -X['Real_Rate']  # Lower real rates = easier
    X['Credit_Spread'] = -X['Credit_Spread']  # Tighter spreads = easier
    X['VIX'] = -X['VIX']  # Lower VIX = easier
    
    # Standardize
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)
    X_scaled_df = pd.DataFrame(X_scaled, index=X.index, columns=X.columns)
    
    # Fit PCA/Sparse PCA
    if method == 'sparse':
        model = SparsePCA(n_components=n_components, alpha=alpha, random_state=42)
    else:
        model = PCA(n_components=n_components)
    
    # Fit and transform
    L = model.fit_transform(X_scaled)
    
    # Create series
    L_series = pd.Series(L[:, 0], index=X.index, name='L')
    
    # Get loadings
    loadings = pd.DataFrame(
        model.components_.T,
        index=X.columns,
        columns=['PC1']
    )
    
    # Calculate explained variance
    if method == 'pca':
        explained_var = model.explained_variance_ratio_[0]
    else:
        # For sparse PCA, compute manually
        var_explained = np.var(L[:, 0]) / np.sum(np.var(X_scaled, axis=0))
        explained_var = var_explained
    
    print(f"\n{'='*70}")
    print(f"Liquidity Index Construction ({method.upper()})")
    print(f"{'='*70}")
    print(f"\nVariance explained: {explained_var:.2%}")
    print(f"\nLoadings:")
    print(loadings.sort_values('PC1', ascending=False))
    print(f"\nLiquidity index range: {L_series.min():.2f} to {L_series.max():.2f}")
    print(f"Mean: {L_series.mean():.2f}, Std: {L_series.std():.2f}")
    
    return L_series, loadings, model, X_scaled_df

# Build liquidity index
L_index, loadings, pca_model, scaled_features = build_liquidity_index_aug(
    df_features_aug,
    method='sparse',
    alpha=0.5  # Moderate sparsity
)

======================================================================
Liquidity Index Construction (SPARSE)
======================================================================

Variance explained: 37.40%

Loadings:
                    PC1
EL_3y          0.479785
FED_BS_growth  0.414107
M2_growth      0.386169
Term_Spread    0.365477
Real_Rate      0.348663
EB             0.135796
EM             0.091370
Credit_Spread -0.242676
VIX           -0.329204

Liquidity index range: -4.23 to 4.97
Mean: -0.00, Std: 1.84

Explained variance is too small = ~45%

def build_liquidity_index_aug_1(df, method='sparse', n_components=1, alpha=1.0):
    """
    Build liquidity index using PCA or Sparse PCA
    
    Parameters:
    -----------
    df : DataFrame
        Input data
    method : str
        'pca' or 'sparse'
    n_components : int
        Number of components
    alpha : float
        Sparsity penalty (only for sparse PCA)
    """
    # Select liquidity variables
    # liq_vars = [
    #     'Net_Liq',
    #     'M2_growth',
    #     'FED_BS_growth',
    #     'EM',
    #     'EB',
    #     'EL_3y',
    #     'Term_Spread',
    #     'Real_Rate',
    #     'Credit_Spread',
    #     'VIX',
    #     'Dollar_Index',
    #     'Household_Debt_Service',
    #     'STLFSI'
    # ]
    liq_vars = [
        'M2_growth',
        'FED_BS_growth',
        'Term_Spread',
        'Real_Rate',
        'Credit_Spread',
        'VIX',
        'EL_3y'
    ]
    
    # Extract data and drop NaN
    X = df[liq_vars].copy()
    X = X.dropna()
    
    # ========================================
    # CORRECT SIGN ADJUSTMENTS
    # Rule: Higher value = EASIER liquidity
    # ========================================
    
    # Variables to FLIP (higher originally means tighter)
    flip_vars = {
        'Real_Rate': 'Higher real rates = tighter → FLIP',
        'Credit_Spread': 'Wider spreads = tighter → FLIP',
        'VIX': 'Higher VIX = tighter → FLIP',
        'STLFSI': 'Higher stress = tighter → FLIP',
        'Dollar_Index': 'Stronger dollar = tighter → FLIP',
        'Household_Debt_Service': 'Higher burden = tighter → FLIP',
    }
    
    # Variables to KEEP (higher already means easier)
    keep_vars = {
        'Net_Liq': 'Higher net liquidity = easier → KEEP',
        'M2_growth': 'Higher growth = easier → KEEP',
        'FED_BS_growth': 'Higher growth = easier → KEEP',
        'EM': 'Higher excess M2 = easier → KEEP',
        'EB': 'Higher excess Fed BS = easier → KEEP',
        'EL_3y': 'Higher excess liquidity = easier → KEEP',
        'Term_Spread': 'Positive slope = easier → KEEP',
    }
    
    print("\n➖ Flipping (higher = tighter):")
    for var, reason in flip_vars.items():
        if var in X.columns:
            X[var] = -X[var]
            print(f"   {var:25s}: {reason}")
    
    print("\n✅ Keeping (higher = easier):")
    for var, reason in keep_vars.items():
        if var in X.columns:
            print(f"   {var:25s}: {reason}")
    
    # Standardize
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)
    X_scaled_df = pd.DataFrame(X_scaled, index=X.index, columns=X.columns)
    
    # Fit PCA/Sparse PCA
    if method == 'sparse':
        model = SparsePCA(n_components=n_components, alpha=alpha, random_state=42)
    else:
        model = PCA(n_components=n_components)
    
    # Fit and transform
    L = model.fit_transform(X_scaled)
    
    # Create series
    L_series = pd.Series(L[:, 0], index=X.index, name='L')
    
    # Get loadings
    loadings = pd.DataFrame(
        model.components_.T,
        index=X.columns,
        columns=['PC1']
    )
    
    # Calculate explained variance
    if method == 'pca':
        explained_var = model.explained_variance_ratio_[0]
    else:
        # For sparse PCA, compute manually
        var_explained = np.var(L[:, 0]) / np.sum(np.var(X_scaled, axis=0))
        explained_var = var_explained
    
    print(f"\n{'='*70}")
    print(f"Liquidity Index Construction ({method.upper()})")
    print(f"{'='*70}")
    print(f"\nVariance explained: {explained_var:.2%}")
    print(f"\nLoadings:")
    print(loadings.sort_values('PC1', ascending=False))
    print(f"\nLiquidity index range: {L_series.min():.2f} to {L_series.max():.2f}")
    print(f"Mean: {L_series.mean():.2f}, Std: {L_series.std():.2f}")
    
    return L_series, loadings, model, X_scaled_df

# Build liquidity index
L_index, loadings, pca_model, scaled_features = build_liquidity_index_aug_1(
    df_features_aug,
    method='sparse',
    alpha=0.5  # Moderate sparsity
)

➖ Flipping (higher = tighter):
   Real_Rate                : Higher real rates = tighter → FLIP
   Credit_Spread            : Wider spreads = tighter → FLIP
   VIX                      : Higher VIX = tighter → FLIP

✅ Keeping (higher = easier):
   M2_growth                : Higher growth = easier → KEEP
   FED_BS_growth            : Higher growth = easier → KEEP
   EL_3y                    : Higher excess liquidity = easier → KEEP
   Term_Spread              : Positive slope = easier → KEEP

======================================================================
Liquidity Index Construction (SPARSE)
======================================================================

Variance explained: 47.21%

Loadings:
                    PC1
EL_3y          0.472608
FED_BS_growth  0.428525
Term_Spread    0.385863
M2_growth      0.381727
Real_Rate      0.334564
Credit_Spread -0.269403
VIX           -0.337477

Liquidity index range: -4.46 to 5.37
Mean: -0.00, Std: 1.82

6. Regime Detection via Hidden Markov Models

6.1 HMM Specification

We employ a 3-state Gaussian HMM to identify: - State 0 (Tight): Low liquidity, restrictive monetary policy - State 1 (Neutral): Normal conditions - State 2 (High): Easy liquidity, expansionary policy

The model assumes: \[L_t | s_t = k \sim \mathcal{N}(\mu_k, \sigma_k^2)\]

With transition matrix capturing regime persistence.

def fit_hmm_regimes(L_series, n_states=3, n_iter=1000):
    """
    Fit Hidden Markov Model to identify liquidity regimes
    
    Parameters:
    -----------
    L_series : Series
        Liquidity index
    n_states : int
        Number of hidden states
    n_iter : int
        Maximum EM iterations
    """
    # Prepare data
    X = L_series.dropna().values.reshape(-1, 1)
    
    # Fit HMM
    model = hmm.GaussianHMM(
        n_components=n_states,
        covariance_type='full',
        n_iter=n_iter,
        random_state=42
    )
    
    model.fit(X)
    
    # Predict states (Viterbi algorithm)
    states = model.predict(X)
    
    # Get posterior probabilities
    posteriors = model.predict_proba(X)
    
    # Create results DataFrame
    results = pd.DataFrame({
        'L': L_series.dropna(),
        'state': states
    })
    
    for i in range(n_states):
        results[f'prob_state_{i}'] = posteriors[:, i]
    
    # Sort states by mean (0=Tight, 1=Neutral, 2=High)
    state_means = results.groupby('state')['L'].mean().sort_values()
    state_mapping = {old: new for new, old in enumerate(state_means.index)}
    results['state'] = results['state'].map(state_mapping)
    
    # Add labels
    state_labels = {0: 'Tight', 1: 'Neutral', 2: 'High'}
    results['state_label'] = results['state'].map(state_labels)
    
    # Print diagnostics
    print(f"\n{'='*70}")
    print(f"HMM Regime Detection ({n_states} states)")
    print(f"{'='*70}")
    print(f"\nModel converged: {model.monitor_.converged}")
    print(f"Log-likelihood: {model.score(X):.2f}")
    
    print(f"\nState statistics:")
    for state in range(n_states):
        mask = results['state'] == state
        n_obs = mask.sum()
        pct = n_obs / len(results) * 100
        mean_L = results.loc[mask, 'L'].mean()
        std_L = results.loc[mask, 'L'].std()
        label = state_labels[state]
        print(f"  State {state} ({label}): {n_obs} obs ({pct:.1f}%), "
              f"L = {mean_L:+.2f} ± {std_L:.2f}")
    
    print(f"\nTransition Matrix:")
    trans_matrix = pd.DataFrame(
        model.transmat_,
        index=[f"From_{state_labels[i]}" for i in range(n_states)],
        columns=[f"To_{state_labels[i]}" for i in range(n_states)]
    )
    print(trans_matrix.round(3))
    
    # Calculate regime persistence (diagonal elements)
    persistence = np.diag(model.transmat_)
    print(f"\nRegime persistence (prob. of staying):")
    for i, label in state_labels.items():
        print(f"  {label}: {persistence[i]:.1%}")
    
    return results, model, state_labels

# Fit HMM
regimes_df, hmm_model, state_labels = fit_hmm_regimes(L_index, n_states=3)

======================================================================
HMM Regime Detection (3 states)
======================================================================

Model converged: True
Log-likelihood: -301.76

State statistics:
  State 0 (Tight): 80 obs (43.2%), L = -1.62 ± 0.99
  State 1 (Neutral): 53 obs (28.6%), L = +1.21 ± 1.24
  State 2 (High): 52 obs (28.1%), L = +1.26 ± 1.29

Transition Matrix:
              To_Tight  To_Neutral  To_High
From_Tight       0.001       0.999    0.000
From_Neutral     0.935       0.009    0.056
From_High        0.022       0.004    0.974

Regime persistence (prob. of staying):
  Tight: 0.1%
  Neutral: 0.9%
  High: 97.4%

Validation

def map_to_known_liquidity_regimes():
    """
    Define known Fed policy regimes based on historical record
    """
    
    # Format: (start_date, end_date, regime_type, description)
    known_regimes = [
        # =====================================
        # Pre-Crisis Era (1990-2007)
        # =====================================
        ('1990-01', '1992-12', 'Neutral', 'Post-1990 Recession, Gradual Easing'),
        ('1993-01', '1994-01', 'High', 'Accommodative Policy, Economic Expansion'),
        ('1994-02', '1995-12', 'Tight', 'Fed Tightening Cycle (3.0% → 6.0%)'),
        ('1996-01', '1998-12', 'Neutral', 'Goldilocks Economy'),
        ('1999-01', '2000-05', 'Tight', 'Dot-com Bubble, Y2K Tightening'),
        ('2000-06', '2003-06', 'High', 'Dot-com Crash Response, Aggressive Easing'),
        ('2003-07', '2004-05', 'High', 'Ultra-Low Rates (1%), Housing Boom'),
        ('2004-06', '2006-06', 'Tight', 'Fed Hiking Cycle (1% → 5.25%)'),
        ('2006-07', '2007-07', 'Neutral', 'Pause at Peak, Pre-Crisis'),
        
        # =====================================
        # Financial Crisis & QE Era (2007-2015)
        # =====================================
        ('2007-08', '2008-08', 'Neutral', 'Initial Crisis Response, Rate Cuts Begin'),
        ('2008-09', '2010-03', 'High', 'QE1: Emergency Response ($1.7T)'),
        ('2010-04', '2010-10', 'Neutral', 'QE1 Taper, Brief Pause'),
        ('2010-11', '2011-06', 'High', 'QE2: $600B Purchase Program'),
        ('2011-07', '2012-08', 'Neutral', 'Operation Twist, Pre-QE3'),
        ('2012-09', '2014-10', 'High', 'QE3: $85B/month, then taper'),
        ('2014-11', '2015-11', 'Neutral', 'Post-QE3, ZIRP Maintained'),
        
        # =====================================
        # Normalization Attempt (2015-2019)
        # =====================================
        ('2015-12', '2018-12', 'Tight', 'Fed Hiking + QT (0.25% → 2.50%)'),
        ('2019-01', '2019-07', 'Neutral', 'Pause, Economic Slowdown'),
        ('2019-08', '2019-12', 'High', 'Insurance Cuts + Repo Crisis Response'),
        
        # =====================================
        # COVID Era (2020-2021)
        # =====================================
        ('2020-01', '2020-02', 'Neutral', 'Pre-COVID'),
        ('2020-03', '2021-10', 'High', 'COVID QE: Unlimited purchases, ZIRP'),
        ('2021-11', '2022-02', 'High', 'Taper Announcement but still easy'),
        
        # =====================================
        # Modern Tightening (2022-2025)
        # =====================================
        ('2022-03', '2023-06', 'Tight', 'Aggressive Hikes (0% → 5.25%)'),
        ('2023-07', '2024-08', 'Tight', 'Higher for Longer, QT Active'),
        ('2024-09', '2025-01', 'Neutral', 'Easing Cycle Begins'),
    ]
    
    # Convert to DataFrame
    regime_df = pd.DataFrame(known_regimes, columns=['start', 'end', 'regime', 'description'])
    regime_df['start'] = pd.to_datetime(regime_df['start'])
    regime_df['end'] = pd.to_datetime(regime_df['end'])
    
    return regime_df


def validate_hmm_against_known_regimes(regimes_df, known_regimes_df):
    """
    Compare HMM-detected regimes to known Fed policy periods
    """
    
    print(f"\n{'='*70}")
    print("REGIME VALIDATION: HMM vs Known Fed Policy")
    print(f"{'='*70}")
    
    # Expand known regimes to monthly observations
    known_monthly = []
    for _, row in known_regimes_df.iterrows():
        date_range = pd.date_range(row['start'], row['end'], freq='M')
        for date in date_range:
            known_monthly.append({
                'date': date,
                'known_regime': row['regime'],
                'period_description': row['description']
            })
    
    known_df = pd.DataFrame(known_monthly).set_index('date')
    
    # Merge with HMM results
    comparison = regimes_df.join(known_df, how='inner')
    comparison = comparison.dropna(subset=['known_regime', 'state_label'])
    
    print(f"\n📊 Overlap period: {comparison.index[0].strftime('%Y-%m')} to {comparison.index[-1].strftime('%Y-%m')}")
    print(f"   Total months: {len(comparison)}")
    
    # =====================================
    # Confusion Matrix
    # =====================================
    
    print(f"\n{'='*70}")
    print("Confusion Matrix: HMM vs Known Regimes")
    print(f"{'='*70}")
    
    from sklearn.metrics import confusion_matrix, accuracy_score, classification_report
    
    # Create confusion matrix
    cm = confusion_matrix(comparison['known_regime'], comparison['state_label'], 
                         labels=['Tight', 'Neutral', 'High'])
    
    cm_df = pd.DataFrame(cm, 
                        index=['Known_Tight', 'Known_Neutral', 'Known_High'],
                        columns=['HMM_Tight', 'HMM_Neutral', 'HMM_High'])
    
    print("\n", cm_df)
    
    # Calculate percentages
    cm_pct = cm_df.div(cm_df.sum(axis=1), axis=0) * 100
    print("\nPercentages (row-wise):")
    print(cm_pct.round(1))
    
    # Overall accuracy
    accuracy = accuracy_score(comparison['known_regime'], comparison['state_label'])
    print(f"\n📈 Overall Accuracy: {accuracy:.1%}")
    
    # =====================================
    # Regime-by-Regime Analysis
    # =====================================
    
    print(f"\n{'='*70}")
    print("Detailed Regime Concordance")
    print(f"{'='*70}")
    
    for known_regime in ['Tight', 'Neutral', 'High']:
        mask = comparison['known_regime'] == known_regime
        subset = comparison[mask]
        
        print(f"\n{known_regime.upper()} REGIME ({mask.sum()} months):")
        
        # How HMM classified these periods
        hmm_dist = subset['state_label'].value_counts(normalize=True) * 100
        for hmm_regime, pct in hmm_dist.sort_index().items():
            status = "✅" if hmm_regime == known_regime else "⚠️"
            print(f"   {status} Classified as {hmm_regime:8s}: {pct:5.1f}%")
    
    # =====================================
    # Key Historical Periods Check
    # =====================================
    
    print(f"\n{'='*70}")
    print("Key Historical Periods Validation")
    print(f"{'='*70}")
    
    key_periods = [
        ('2008-09', '2009-12', 'High', 'QE1 (Financial Crisis)'),
        ('2010-11', '2011-06', 'High', 'QE2'),
        ('2012-09', '2014-10', 'High', 'QE3'),
        ('2015-12', '2018-12', 'Tight', 'Fed Hiking + QT'),
        ('2020-03', '2021-06', 'High', 'COVID QE (Unlimited)'),
        ('2022-03', '2023-12', 'Tight', 'Aggressive Rate Hikes'),
    ]
    
    for start, end, expected, description in key_periods:
        try:
            mask = (comparison.index >= start) & (comparison.index <= end)
            subset = comparison[mask]
            
            if len(subset) == 0:
                continue
            
            # Modal HMM regime
            hmm_mode = subset['state_label'].mode()[0]
            hmm_pct = (subset['state_label'] == hmm_mode).sum() / len(subset) * 100
            
            # Check if matches expectation
            match = "✅" if hmm_mode == expected else "❌"
            
            print(f"\n{match} {description}")
            print(f"   Period: {start} to {end}")
            print(f"   Expected: {expected}")
            print(f"   HMM Mode: {hmm_mode} ({hmm_pct:.0f}% of period)")
            
            # Show distribution
            dist = subset['state_label'].value_counts()
            print(f"   Distribution: ", end="")
            for regime, count in dist.items():
                print(f"{regime}={count}, ", end="")
            print()
            
        except Exception as e:
            print(f"\n⚠️  {description}: Data not available")
    
    # =====================================
    # Mismatches Analysis
    # =====================================
    
    print(f"\n{'='*70}")
    print("Significant Mismatches (Known ≠ HMM)")
    print(f"{'='*70}")
    
    mismatches = comparison[comparison['known_regime'] != comparison['state_label']]
    
    if len(mismatches) > 0:
        # Find consecutive mismatch periods
        mismatch_periods = []
        current_start = None
        
        for i, (date, row) in enumerate(mismatches.iterrows()):
            if current_start is None:
                current_start = date
                current_known = row['known_regime']
                current_hmm = row['state_label']
            
            # Check if next row continues the mismatch pattern
            is_last = i == len(mismatches) - 1
            if is_last or (not is_last and mismatches.index[i+1] != date + pd.DateOffset(months=1)):
                mismatch_periods.append({
                    'start': current_start,
                    'end': date,
                    'known': current_known,
                    'hmm': current_hmm,
                    'months': (date.year - current_start.year) * 12 + date.month - current_start.month + 1
                })
                current_start = None
        
        # Show significant mismatches (>= 3 months)
        significant = [p for p in mismatch_periods if p['months'] >= 3]
        
        if significant:
            print(f"\nFound {len(significant)} significant mismatch periods (≥3 months):\n")
            for p in significant:
                print(f"   {p['start'].strftime('%Y-%m')} to {p['end'].strftime('%Y-%m')} ({p['months']} months)")
                print(f"      Known: {p['known']:8s} | HMM: {p['hmm']:8s}")
                
                # Get description for this period
                desc_mask = (known_regimes_df['start'] <= p['start']) & (known_regimes_df['end'] >= p['end'])
                if desc_mask.any():
                    desc = known_regimes_df[desc_mask]['description'].iloc[0]
                    print(f"      Context: {desc}")
                print()
        else:
            print("\n✅ No significant consecutive mismatches (all discrepancies < 3 months)")
    
    else:
        print("\n🎉 Perfect match! HMM exactly matches known regimes.")
    
    # =====================================
    # Summary Statistics
    # =====================================
    
    print(f"\n{'='*70}")
    print("SUMMARY STATISTICS")
    print(f"{'='*70}")
    
    print(f"\n📊 Agreement Rates by Regime:")
    for regime in ['Tight', 'Neutral', 'High']:
        mask = comparison['known_regime'] == regime
        if mask.sum() > 0:
            agreement = (comparison[mask]['state_label'] == regime).sum() / mask.sum() * 100
            print(f"   {regime:8s}: {agreement:5.1f}% correct")
    
    print(f"\n📊 Overall Performance:")
    print(f"   Total Accuracy:        {accuracy:.1%}")
    print(f"   Mismatch Rate:         {(1-accuracy):.1%}")
    print(f"   Consecutive Mismatches: {len([p for p in mismatch_periods if p['months'] >= 3])}")
    
    # Interpretation
    print(f"\n{'='*70}")
    print("INTERPRETATION")
    print(f"{'='*70}")
    
    if accuracy >= 0.80:
        print("\n✅ EXCELLENT: HMM successfully captures Fed policy regimes")
        print("   The model shows strong alignment with known monetary policy history.")
    elif accuracy >= 0.65:
        print("\n✔️  GOOD: HMM generally aligns with Fed policy regimes")
        print("   Some discrepancies exist but the model captures major policy shifts.")
    elif accuracy >= 0.50:
        print("\n⚠️  MODERATE: HMM shows partial alignment with Fed policy")
        print("   Consider adjusting liquidity index composition or HMM parameters.")
    else:
        print("\n❌ POOR: HMM does not align well with Fed policy regimes")
        print("   Major revisions needed to liquidity index or regime detection.")
    
    return comparison, cm_df, accuracy


# Run validation
known_regimes_df = map_to_known_liquidity_regimes()
comparison_df, confusion_matrix, accuracy = validate_hmm_against_known_regimes(
    regimes_df, 
    known_regimes_df
)

======================================================================
REGIME VALIDATION: HMM vs Known Fed Policy
======================================================================

📊 Overlap period: 2003-12 to 2024-12
   Total months: 166

======================================================================
Confusion Matrix: HMM vs Known Regimes
======================================================================

                HMM_Tight  HMM_Neutral  HMM_High
Known_Tight           44           11         9
Known_Neutral         21           12        13
Known_High             2           28        26

Percentages (row-wise):
               HMM_Tight  HMM_Neutral  HMM_High
Known_Tight         68.8         17.2      14.1
Known_Neutral       45.7         26.1      28.3
Known_High           3.6         50.0      46.4

📈 Overall Accuracy: 49.4%

======================================================================
Detailed Regime Concordance
======================================================================

TIGHT REGIME (64 months):
   ⚠️ Classified as High    :  14.1%
   ⚠️ Classified as Neutral :  17.2%
   ✅ Classified as Tight   :  68.8%

NEUTRAL REGIME (46 months):
   ⚠️ Classified as High    :  28.3%
   ✅ Classified as Neutral :  26.1%
   ⚠️ Classified as Tight   :  45.7%

HIGH REGIME (56 months):
   ✅ Classified as High    :  46.4%
   ⚠️ Classified as Neutral :  50.0%
   ⚠️ Classified as Tight   :   3.6%

======================================================================
Key Historical Periods Validation
======================================================================

✅ QE1 (Financial Crisis)
   Period: 2008-09 to 2009-12
   Expected: High
   HMM Mode: High (50% of period)
   Distribution: Neutral=5, High=5, 

✅ QE2
   Period: 2010-11 to 2011-06
   Expected: High
   HMM Mode: High (50% of period)
   Distribution: Neutral=3, High=3, 

✅ QE3
   Period: 2012-09 to 2014-10
   Expected: High
   HMM Mode: High (50% of period)
   Distribution: High=9, Neutral=9, 

✅ Fed Hiking + QT
   Period: 2015-12 to 2018-12
   Expected: Tight
   HMM Mode: Tight (58% of period)
   Distribution: Tight=15, Neutral=6, High=5, 

✅ COVID QE (Unlimited)
   Period: 2020-03 to 2021-06
   Expected: High
   HMM Mode: High (50% of period)
   Distribution: High=5, Neutral=5, 

✅ Aggressive Rate Hikes
   Period: 2022-03 to 2023-12
   Expected: Tight
   HMM Mode: Tight (53% of period)
   Distribution: Tight=8, Neutral=4, High=3, 

======================================================================
Significant Mismatches (Known ≠ HMM)
======================================================================

Found 1 significant mismatch periods (≥3 months):

   2007-12 to 2008-02 (3 months)
      Known: Neutral  | HMM: Tight   
      Context: Initial Crisis Response, Rate Cuts Begin


======================================================================
SUMMARY STATISTICS
======================================================================

📊 Agreement Rates by Regime:
   Tight   :  68.8% correct
   Neutral :  26.1% correct
   High    :  46.4% correct

📊 Overall Performance:
   Total Accuracy:        49.4%
   Mismatch Rate:         50.6%
   Consecutive Mismatches: 1

======================================================================
INTERPRETATION
======================================================================

❌ POOR: HMM does not align well with Fed policy regimes
   Major revisions needed to liquidity index or regime detection.
def build_advanced_regime_classification_system(df, L_series):
    """
    State-of-the-art regime classification with proper index alignment
    """
    
    import numpy as np
    import pandas as pd
    from sklearn.ensemble import IsolationForest, RandomForestClassifier
    from xgboost import XGBClassifier
    from sklearn.model_selection import TimeSeriesSplit, cross_val_score
    from sklearn.metrics import classification_report, confusion_matrix, roc_auc_score
    import shap
    import matplotlib.pyplot as plt
    import seaborn as sns
    
    print(f"\n{'='*80}")
    print("ADVANCED REGIME CLASSIFICATION SYSTEM")
    print(f"{'='*80}")
    
    # ==========================================
    # STEP 0: DATA ALIGNMENT
    # ==========================================
    
    print(f"\n{'='*80}")
    print("STEP 0: DATA ALIGNMENT")
    print(f"{'='*80}")
    
    # Find common index
    common_idx = df.index.intersection(L_series.index)
    
    if len(common_idx) == 0:
        raise ValueError("No overlapping dates between df and L_series!")
    
    print(f"\n📊 Data Coverage:")
    print(f"   df range:       {df.index[0]} to {df.index[-1]} ({len(df)} obs)")
    print(f"   L_series range: {L_series.index[0]} to {L_series.index[-1]} ({len(L_series)} obs)")
    print(f"   Overlap range:  {common_idx[0]} to {common_idx[-1]} ({len(common_idx)} obs)")
    
    # Use only overlapping data
    df_aligned = df.loc[common_idx].copy()
    L_aligned = L_series.loc[common_idx].copy()
    
    # Add L to dataframe
    df_aligned['L'] = L_aligned
    
    # Check for required columns
    required_cols = ['mkt_excess']
    missing_cols = [c for c in required_cols if c not in df_aligned.columns]
    if missing_cols:
        raise ValueError(f"Missing required columns: {missing_cols}")
    
    print(f"\n✅ Data aligned: {len(df_aligned)} observations")
    
    # ==========================================
    # STEP 1: Crisis Detection (Isolation Forest)
    # ==========================================
    
    print(f"\n{'='*80}")
    print("STEP 1: CRISIS DETECTION")
    print(f"{'='*80}")
    
    # Create crisis features
    df_crisis = df_aligned.copy()
    
    # Volatility features
    df_crisis['ret_vol_6m'] = df_crisis['mkt_excess'].rolling(6, min_periods=3).std()
    df_crisis['ret_vol_12m'] = df_crisis['mkt_excess'].rolling(12, min_periods=6).std()
    df_crisis['L_vol_3m'] = df_crisis['L'].rolling(3, min_periods=2).std()
    df_crisis['L_vol_6m'] = df_crisis['L'].rolling(6, min_periods=3).std()
    
    # Return features
    df_crisis['abs_ret'] = df_crisis['mkt_excess'].abs()
    df_crisis['cum_ret_3m'] = df_crisis['mkt_excess'].rolling(3, min_periods=2).sum()
    df_crisis['cum_ret_6m'] = df_crisis['mkt_excess'].rolling(6, min_periods=3).sum()
    
    # Drawdown
    cumret = (1 + df_crisis['mkt_excess']/100).cumprod()
    running_max = cumret.expanding().max()
    df_crisis['drawdown'] = (cumret / running_max - 1) * 100
    df_crisis['max_dd_6m'] = df_crisis['drawdown'].rolling(6, min_periods=3).min()
    
    # VIX if available
    crisis_features = ['ret_vol_6m', 'ret_vol_12m', 'L_vol_3m', 'L_vol_6m', 
                       'abs_ret', 'cum_ret_3m', 'cum_ret_6m', 'max_dd_6m']
    
    if 'VIX' in df_crisis.columns:
        df_crisis['VIX_level'] = df_crisis['VIX']
        df_crisis['VIX_zscore'] = (df_crisis['VIX'] - df_crisis['VIX'].rolling(60, min_periods=30).mean()) / \
                                   df_crisis['VIX'].rolling(60, min_periods=30).std()
        crisis_features.extend(['VIX_level', 'VIX_zscore'])
    
    # Check data availability
    print(f"\n📊 Crisis Features:")
    for feat in crisis_features:
        if feat in df_crisis.columns:
            n_valid = df_crisis[feat].notna().sum()
            print(f"   {feat:20s}: {n_valid}/{len(df_crisis)} valid")
    
    # Extract crisis features and drop NaN
    crisis_data = df_crisis[crisis_features].copy()
    
    # Fill forward initial NaN from rolling windows (up to 12 months)
    crisis_data = crisis_data.fillna(method='bfill', limit=12)
    crisis_data = crisis_data.dropna()
    
    print(f"\n📊 Crisis detection dataset: {len(crisis_data)} observations")
    
    if len(crisis_data) < 20:
        print(f"\n⚠️  WARNING: Not enough data for crisis detection ({len(crisis_data)} obs)")
        print(f"   Skipping crisis detection, will use all data for regime classification")
        
        df_crisis['is_crisis'] = False
        df_noncrisis = df_crisis.copy()
        
    else:
        # Fit Isolation Forest
        iso_forest = IsolationForest(
            contamination=0.05,
            random_state=42,
            n_estimators=100,
            max_samples='auto'
        )
        
        outlier_labels = iso_forest.fit_predict(crisis_data.values)
        
        # Crisis = outlier + negative return
        is_crisis_candidate = (outlier_labels == -1)
        is_negative_return = (df_crisis.loc[crisis_data.index, 'mkt_excess'] < -2)
        
        df_crisis['is_crisis'] = False
        df_crisis.loc[crisis_data.index, 'is_crisis'] = is_crisis_candidate & is_negative_return
        
        n_crisis = df_crisis['is_crisis'].sum()
        crisis_pct = n_crisis / len(df_crisis) * 100
        
        print(f"\n✅ Crisis Detection Complete")
        print(f"   Crisis periods:     {n_crisis} months ({crisis_pct:.1f}%)")
        print(f"   Non-crisis periods: {len(df_crisis) - n_crisis} months ({100-crisis_pct:.1f}%)")
        
        # Show crisis periods
        if n_crisis > 0:
            crisis_periods = df_crisis[df_crisis['is_crisis']].index
            print(f"\n📅 Detected Crisis Months:")
            for i, date in enumerate(crisis_periods[:10], 1):
                ret = df_crisis.loc[date, 'mkt_excess']
                print(f"   {i:2d}. {date.strftime('%Y-%m')}: Return = {ret:+.2f}%")
            if len(crisis_periods) > 10:
                print(f"   ... and {len(crisis_periods) - 10} more")
        
        # Filter to non-crisis
        df_noncrisis = df_crisis[~df_crisis['is_crisis']].copy()
    
    # ==========================================
    # STEP 2: Regime Classification Setup
    # ==========================================
    
    print(f"\n{'='*80}")
    print("STEP 2: REGIME CLASSIFICATION (NON-CRISIS PERIODS)")
    print(f"{'='*80}")
    
    # Create regime labels based on liquidity quantiles
    L_noncrisis = df_noncrisis['L']
    
    if len(L_noncrisis) < 30:
        raise ValueError(f"Not enough non-crisis data: {len(L_noncrisis)} observations")
    
    q33 = L_noncrisis.quantile(0.33)
    q67 = L_noncrisis.quantile(0.67)
    
    print(f"\n📊 Liquidity Quantiles:")
    print(f"   33rd percentile (q33): {q33:+.2f}")
    print(f"   67th percentile (q67): {q67:+.2f}")
    
    def assign_regime(L_val):
        if L_val < q33:
            return 0  # Tight
        elif L_val < q67:
            return 1  # Neutral
        else:
            return 2  # High
    
    df_noncrisis['regime'] = L_noncrisis.apply(assign_regime)
    df_noncrisis['regime_label'] = df_noncrisis['regime'].map({0: 'Tight', 1: 'Neutral', 2: 'High'})
    
    print(f"\n📊 Non-Crisis Regime Distribution:")
    for regime, label in [(0, 'Tight'), (1, 'Neutral'), (2, 'High')]:
        mask = df_noncrisis['regime'] == regime
        count = mask.sum()
        pct = count / len(df_noncrisis) * 100
        mean_L = df_noncrisis[mask]['L'].mean()
        print(f"   {label:8s}: {count:4d} months ({pct:5.1f}%), Mean L = {mean_L:+.2f}")
    
    # ==========================================
    # STEP 3: Feature Engineering
    # ==========================================
    
    print(f"\n{'='*80}")
    print("STEP 3: FEATURE ENGINEERING")
    print(f"{'='*80}")
    
    # Build feature set
    feature_cols = []
    
    # Core liquidity components
    liq_features = ['dlog_M2', 'dlog_FED_BS', 'EM', 'EB', 'EL_3y']
    if 'Net_Liq' in df_noncrisis.columns:
        liq_features.append('Net_Liq')
    
    for feat in liq_features:
        if feat in df_noncrisis.columns:
            feature_cols.append(feat)
    
    # Rates and spreads
    rate_features = ['Term_Spread', 'Credit_Spread']
    if 'Real_Rate' in df_noncrisis.columns:
        rate_features.append('Real_Rate')
    
    for feat in rate_features:
        if feat in df_noncrisis.columns:
            feature_cols.append(feat)
    
    # Market indicators
    if 'VIX' in df_noncrisis.columns:
        feature_cols.append('VIX')
    
    # Growth rates
    if 'M2_growth' in df_noncrisis.columns:
        feature_cols.append('M2_growth')
    if 'FED_BS_growth' in df_noncrisis.columns:
        feature_cols.append('FED_BS_growth')
    
    # Technical features on L
    df_noncrisis['L_ma_3'] = df_noncrisis['L'].rolling(3, min_periods=1).mean()
    df_noncrisis['L_ma_6'] = df_noncrisis['L'].rolling(6, min_periods=3).mean()
    df_noncrisis['L_momentum'] = df_noncrisis['L'].diff(3)
    
    feature_cols.extend(['L_ma_3', 'L_ma_6', 'L_momentum'])
    
    # Valuation if available
    if 'V_spread' in df_noncrisis.columns:
        feature_cols.append('V_spread')
    
    print(f"\n✅ Selected {len(feature_cols)} features:")
    for i, feat in enumerate(feature_cols, 1):
        if feat in df_noncrisis.columns:
            n_valid = df_noncrisis[feat].notna().sum()
            print(f"   {i:2d}. {feat:20s}: {n_valid}/{len(df_noncrisis)} valid")
    
    # Prepare feature matrix
    X = df_noncrisis[feature_cols].copy()
    y = df_noncrisis['regime'].copy()
    
    # Forward fill initial NaN from rolling windows
    X = X.fillna(method='ffill', limit=6)
    
    # Drop remaining NaN
    valid_idx = X.notna().all(axis=1) & y.notna()
    X = X[valid_idx]
    y = y[valid_idx]
    
    print(f"\n📊 Final ML dataset: {len(X)} samples × {len(feature_cols)} features")
    
    if len(X) < 30:
        raise ValueError(f"Not enough valid data for ML: {len(X)} samples")
    
    # ==========================================
    # STEP 4: Train ML Models
    # ==========================================
    
    print(f"\n{'='*80}")
    print("STEP 4: MODEL TRAINING")
    print(f"{'='*80}")
    
    # Time series split
    tscv = TimeSeriesSplit(n_splits=min(5, len(X) // 20))
    
    # Models
    models = {
        'Random Forest': RandomForestClassifier(
            n_estimators=100,
            max_depth=8,
            min_samples_split=10,
            min_samples_leaf=5,
            class_weight='balanced',
            random_state=42,
            n_jobs=-1
        ),
        'XGBoost': XGBClassifier(
            n_estimators=100,
            max_depth=5,
            learning_rate=0.1,
            subsample=0.8,
            colsample_bytree=0.8,
            random_state=42,
            eval_metric='mlogloss',
            use_label_encoder=False,
            n_jobs=-1
        )
    }
    
    results = {}
    
    for name, model in models.items():
        print(f"\n{'='*70}")
        print(f"Training {name}")
        print(f"{'='*70}")
        
        # Cross-validation
        try:
            cv_scores = cross_val_score(model, X, y, cv=tscv, scoring='accuracy', n_jobs=-1)
            print(f"\n📊 CV Accuracy: {cv_scores.mean():.3f}{cv_scores.std():.3f})")
        except Exception as e:
            print(f"\n⚠️  CV failed: {e}")
            cv_scores = np.array([])
        
        # Train full model
        model.fit(X, y)
        
        # Predictions
        y_pred = model.predict(X)
        y_pred_proba = model.predict_proba(X)
        
        # Metrics
        accuracy = (y_pred == y).mean()
        
        try:
            roc_auc = roc_auc_score(y, y_pred_proba, multi_class='ovr', average='weighted')
        except:
            roc_auc = np.nan
        
        print(f"\n📈 Performance:")
        print(f"   Accuracy: {accuracy:.3f}")
        if not np.isnan(roc_auc):
            print(f"   ROC-AUC:  {roc_auc:.3f}")
        
        # Classification report
        print(f"\n📋 Classification Report:")
        print(classification_report(y, y_pred, target_names=['Tight', 'Neutral', 'High'], zero_division=0))
        
        # Confusion matrix
        cm = confusion_matrix(y, y_pred)
        cm_pct = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis] * 100
        
        print(f"\n🎯 Confusion Matrix (%):")
        cm_df = pd.DataFrame(
            cm_pct,
            index=['True_Tight', 'True_Neutral', 'True_High'],
            columns=['Pred_Tight', 'Pred_Neutral', 'Pred_High']
        )
        print(cm_df.round(1))
        
        results[name] = {
            'model': model,
            'cv_scores': cv_scores,
            'accuracy': accuracy,
            'roc_auc': roc_auc,
            'y_pred': y_pred,
            'y_pred_proba': y_pred_proba
        }
    
    # ==========================================
    # STEP 5: Feature Importance
    # ==========================================
    
    print(f"\n{'='*80}")
    print("STEP 5: FEATURE IMPORTANCE")
    print(f"{'='*80}")
    
    best_model = results['XGBoost']['model']
    
    feature_importance = pd.DataFrame({
        'feature': X.columns,
        'importance': best_model.feature_importances_
    }).sort_values('importance', ascending=False)
    
    print(f"\n📊 Top 10 Features:")
    for i, row in feature_importance.head(10).iterrows():
        print(f"   {row['feature']:20s}: {row['importance']:.4f}")
    
    # ==========================================
    # STEP 6: Final Labels
    # ==========================================
    
    print(f"\n{'='*80}")
    print("STEP 6: FINAL REGIME ASSIGNMENT")
    print(f"{'='*80}")
    
    df_final = df_crisis.copy()
    df_final['regime'] = np.nan
    df_final['regime_label'] = 'Unknown'
    
    # Crisis
    if 'is_crisis' in df_final.columns:
        df_final.loc[df_final['is_crisis'], 'regime'] = 3
        df_final.loc[df_final['is_crisis'], 'regime_label'] = 'Crisis'
    
    # ML predictions
    ml_idx = X.index
    df_final.loc[ml_idx, 'regime'] = results['XGBoost']['y_pred']
    df_final.loc[ml_idx, 'regime_label'] = pd.Series(
        results['XGBoost']['y_pred'],
        index=ml_idx
    ).map({0: 'Tight', 1: 'Neutral', 2: 'High'})
    
    print(f"\n📊 Final Distribution:")
    dist = df_final['regime_label'].value_counts()
    for regime, count in dist.items():
        pct = count / len(df_final) * 100
        print(f"   {regime:8s}: {count:4d} ({pct:5.1f}%)")
    
    print(f"\n{'='*80}")
    print("✅ COMPLETE")
    print(f"{'='*80}")
    
    return df_final, results, feature_importance


# Run with proper inputs
df_regimes_advanced, ml_results, feature_importance = build_advanced_regime_classification_system(
    df_features_aug,
    L_index
)

================================================================================
ADVANCED REGIME CLASSIFICATION SYSTEM
================================================================================

================================================================================
STEP 0: DATA ALIGNMENT
================================================================================

📊 Data Coverage:
   df range:       1990-01-31 00:00:00 to 2025-11-30 00:00:00 (431 obs)
   L_series range: 2003-12-31 00:00:00 to 2025-09-30 00:00:00 (185 obs)
   Overlap range:  2003-12-31 00:00:00 to 2025-09-30 00:00:00 (185 obs)

✅ Data aligned: 185 observations

================================================================================
STEP 1: CRISIS DETECTION
================================================================================

📊 Crisis Features:
   ret_vol_6m          : 183/185 valid
   ret_vol_12m         : 180/185 valid
   L_vol_3m            : 184/185 valid
   L_vol_6m            : 183/185 valid
   abs_ret             : 185/185 valid
   cum_ret_3m          : 184/185 valid
   cum_ret_6m          : 183/185 valid
   max_dd_6m           : 183/185 valid
   VIX_level           : 185/185 valid
   VIX_zscore          : 156/185 valid

📊 Crisis detection dataset: 168 observations

✅ Crisis Detection Complete
   Crisis periods:     0 months (0.0%)
   Non-crisis periods: 185 months (100.0%)

================================================================================
STEP 2: REGIME CLASSIFICATION (NON-CRISIS PERIODS)
================================================================================

📊 Liquidity Quantiles:
   33rd percentile (q33): -0.80
   67th percentile (q67): +0.69

📊 Non-Crisis Regime Distribution:
   Tight   :   61 months ( 33.0%), Mean L = -1.95
   Neutral :   63 months ( 34.1%), Mean L = +0.00
   High    :   61 months ( 33.0%), Mean L = +1.95

================================================================================
STEP 3: FEATURE ENGINEERING
================================================================================

✅ Selected 16 features:
    1. dlog_M2             : 185/185 valid
    2. dlog_FED_BS         : 185/185 valid
    3. EM                  : 185/185 valid
    4. EB                  : 185/185 valid
    5. EL_3y               : 185/185 valid
    6. Net_Liq             : 118/185 valid
    7. Term_Spread         : 185/185 valid
    8. Credit_Spread       : 185/185 valid
    9. Real_Rate           : 185/185 valid
   10. VIX                 : 185/185 valid
   11. M2_growth           : 185/185 valid
   12. FED_BS_growth       : 185/185 valid
   13. L_ma_3              : 185/185 valid
   14. L_ma_6              : 183/185 valid
   15. L_momentum          : 182/185 valid
   16. V_spread            : 185/185 valid

📊 Final ML dataset: 155 samples × 16 features

================================================================================
STEP 4: MODEL TRAINING
================================================================================

======================================================================
Training Random Forest
======================================================================

📊 CV Accuracy: 0.800 (±0.104)

📈 Performance:
   Accuracy: 0.935
   ROC-AUC:  0.995

📋 Classification Report:
              precision    recall  f1-score   support

       Tight       0.95      0.95      0.95        42
     Neutral       0.89      0.93      0.91        54
        High       0.96      0.93      0.95        59

    accuracy                           0.94       155
   macro avg       0.94      0.94      0.94       155
weighted avg       0.94      0.94      0.94       155


🎯 Confusion Matrix (%):
              Pred_Tight  Pred_Neutral  Pred_High
True_Tight          95.2           4.8        0.0
True_Neutral         3.7          92.6        3.7
True_High            0.0           6.8       93.2

======================================================================
Training XGBoost
======================================================================

📊 CV Accuracy: 0.744 (±0.151)

📈 Performance:
   Accuracy: 1.000
   ROC-AUC:  1.000

📋 Classification Report:
              precision    recall  f1-score   support

       Tight       1.00      1.00      1.00        42
     Neutral       1.00      1.00      1.00        54
        High       1.00      1.00      1.00        59

    accuracy                           1.00       155
   macro avg       1.00      1.00      1.00       155
weighted avg       1.00      1.00      1.00       155


🎯 Confusion Matrix (%):
              Pred_Tight  Pred_Neutral  Pred_High
True_Tight         100.0           0.0        0.0
True_Neutral         0.0         100.0        0.0
True_High            0.0           0.0      100.0

================================================================================
STEP 5: FEATURE IMPORTANCE
================================================================================

📊 Top 10 Features:
   L_ma_6              : 0.3007
   L_ma_3              : 0.2085
   Term_Spread         : 0.1153
   FED_BS_growth       : 0.0552
   EL_3y               : 0.0481
   EM                  : 0.0375
   L_momentum          : 0.0371
   Real_Rate           : 0.0302
   Credit_Spread       : 0.0288
   VIX                 : 0.0272

================================================================================
STEP 6: FINAL REGIME ASSIGNMENT
================================================================================

📊 Final Distribution:
   High    :   59 ( 31.9%)
   Neutral :   54 ( 29.2%)
   Tight   :   42 ( 22.7%)
   Unknown :   30 ( 16.2%)

================================================================================
✅ COMPLETE
================================================================================
df_regimes_advanced.tail(2)
M2 FED_BS FFR Real_Rate T3M T10Y BAA AAA CPI CORE_CPI ... abs_ret cum_ret_3m cum_ret_6m drawdown max_dd_6m VIX_level VIX_zscore is_crisis regime regime_label
2025-07-31 22028.7 6642578.0 4.33 0.838904 4.24 4.37 6.04 5.41 322.132 328.656 ... 0.0198 0.0600 -0.0002 -0.041738 -0.110072 16.72 -0.536469 False 0.0 Tight
2025-09-30 22212.4 6608395.0 4.09 1.149868 3.86 4.16 5.83 5.22 324.368 330.542 ... 0.0339 0.1023 0.0057 -0.007852 -0.110072 16.28 -0.613043 False 0.0 Tight

2 rows × 68 columns

# factor_returns_by_regime = (
#     df_regimes_advanced.groupby("regime_label")[["mkt_excess", "hml", "rmw", "cma"]]
#     .mean() * 100
# )
# print(factor_returns_by_regime)
def diagnose_unknown_regimes(df_regimes_advanced, df_features_aug, L_index):
    """
    Understand why some periods are marked as Unknown
    """
    
    print(f"\n{'='*80}")
    print("DIAGNOSIS: Unknown Regime Periods")
    print(f"{'='*80}")
    
    # Find Unknown periods
    unknown_mask = df_regimes_advanced['regime_label'] == 'Unknown'
    unknown_dates = df_regimes_advanced[unknown_mask].index
    
    print(f"\n📊 Unknown Periods: {len(unknown_dates)} months ({len(unknown_dates)/len(df_regimes_advanced)*100:.1f}%)")
    
    if len(unknown_dates) > 0:
        print(f"\n📅 Unknown Date Range:")
        print(f"   First: {unknown_dates[0].strftime('%Y-%m')}")
        print(f"   Last:  {unknown_dates[-1].strftime('%Y-%m')}")
        
        # Check which features are missing
        print(f"\n🔍 Missing Data Analysis:")
        
        feature_cols = ['dlog_M2', 'dlog_FED_BS', 'EM', 'EB', 'EL_3y', 'Net_Liq',
                       'Term_Spread', 'Credit_Spread', 'Real_Rate', 'VIX',
                       'M2_growth', 'FED_BS_growth', 'V_spread']
        
        for feat in feature_cols:
            if feat in df_features_aug.columns:
                # Check missing in unknown periods
                unknown_data = df_features_aug.loc[unknown_dates, feat]
                n_missing = unknown_data.isna().sum()
                pct_missing = n_missing / len(unknown_dates) * 100
                
                if n_missing > 0:
                    print(f"   {feat:20s}: {n_missing:3d}/{len(unknown_dates)} missing ({pct_missing:5.1f}%)")
        
        # Check if Unknown periods cluster at start/end
        print(f"\n📍 Unknown Period Distribution:")
        
        # Group consecutive unknowns
        unknown_groups = []
        current_group = [unknown_dates[0]]
        
        for i in range(1, len(unknown_dates)):
            if (unknown_dates[i] - unknown_dates[i-1]).days < 60:  # Within 2 months
                current_group.append(unknown_dates[i])
            else:
                unknown_groups.append(current_group)
                current_group = [unknown_dates[i]]
        unknown_groups.append(current_group)
        
        print(f"   Found {len(unknown_groups)} clusters of unknown periods:")
        for i, group in enumerate(unknown_groups, 1):
            print(f"   Cluster {i}: {group[0].strftime('%Y-%m')} to {group[-1].strftime('%Y-%m')} ({len(group)} months)")

# Run diagnosis
diagnose_unknown_regimes(df_regimes_advanced, df_features_aug, L_index)

================================================================================
DIAGNOSIS: Unknown Regime Periods
================================================================================

📊 Unknown Periods: 30 months (16.2%)

📅 Unknown Date Range:
   First: 2003-12
   Last:  2012-10

🔍 Missing Data Analysis:
   Net_Liq             :  29/30 missing ( 96.7%)

📍 Unknown Period Distribution:
   Found 14 clusters of unknown periods:
   Cluster 1: 2003-12 to 2003-12 (1 months)
   Cluster 2: 2004-03 to 2004-04 (2 months)
   Cluster 3: 2004-12 to 2005-03 (4 months)
   Cluster 4: 2005-05 to 2005-06 (2 months)
   Cluster 5: 2005-08 to 2005-11 (4 months)
   Cluster 6: 2006-01 to 2006-03 (3 months)
   Cluster 7: 2006-05 to 2006-08 (4 months)
   Cluster 8: 2006-10 to 2006-11 (2 months)
   Cluster 9: 2007-01 to 2007-02 (2 months)
   Cluster 10: 2008-02 to 2008-02 (1 months)
   Cluster 11: 2009-11 to 2009-11 (1 months)
   Cluster 12: 2012-05 to 2012-05 (1 months)
   Cluster 13: 2012-07 to 2012-08 (2 months)
   Cluster 14: 2012-10 to 2012-10 (1 months)
df_features_aug.columns
Index(['M2', 'FED_BS', 'FFR', 'Real_Rate', 'T3M', 'T10Y', 'BAA', 'AAA', 'CPI',
       'CORE_CPI', 'GDP', 'IORB', 'RRP', 'TGA', 'Bank_Reserves', 'CI_Loans',
       'Consumer_Credit', 'Mortgage_Rate', 'Prime_Rate', 'STLFSI',
       'Dollar_Index', 'Household_Debt_Service', 'mkt_excess', 'smb', 'hml',
       'rmw', 'cma', 'rf', 'P', 'CAPE', 'VIX', 'M2_growth', 'FED_BS_growth',
       'Term_Spread', 'Credit_Spread', 'CPI_inflation', 'log_CAPE',
       'CAPE_zscore', 'Risk_Adj_Spread', 'Net_Liq', 'HML_cumret', 'HML_12m',
       'V_spread', 'V_spread_z', 'V_level', 'dlog_M2', 'dlog_FED_BS', 'log_M2',
       'log_FED_BS', 'EM', 'EB', 'log_GDP', 'EL_3y'],
      dtype='object')
def build_hml_regime_classification_with_crisis(df, L_series):
    """
    1. Detect crisis periods (Isolation Forest)
    2. Classify NON-CRISIS periods by HML regimes
    3. Predict HML regimes from liquidity features
    """
    
    import numpy as np
    import pandas as pd
    from sklearn.ensemble import IsolationForest, RandomForestClassifier
    from xgboost import XGBClassifier
    from sklearn.model_selection import TimeSeriesSplit, cross_val_score
    from sklearn.metrics import classification_report, confusion_matrix
    
    print(f"\n{'='*80}")
    print("HML REGIME CLASSIFICATION (WITH CRISIS DETECTION)")
    print(f"{'='*80}")
    
    # ==========================================
    # STEP 0: Data Alignment
    # ==========================================
    
    common_idx = df.index.intersection(L_series.index)
    df_aligned = df.loc[common_idx].copy()
    df_aligned['L'] = L_series.loc[common_idx]
    
    print(f"\n✅ Aligned: {len(df_aligned)} observations")
    
    # ==========================================
    # STEP 1: Crisis Detection (Isolation Forest)
    # ==========================================
    
    print(f"\n{'='*80}")
    print("STEP 1: CRISIS DETECTION")
    print(f"{'='*80}")
    
    df_crisis = df_aligned.copy()
    
    # Create crisis features (volatility, drawdowns, tail events)
    df_crisis['ret_vol_6m'] = df_crisis['mkt_excess'].rolling(6, min_periods=3).std()
    df_crisis['ret_vol_12m'] = df_crisis['mkt_excess'].rolling(12, min_periods=6).std()
    df_crisis['L_vol_3m'] = df_crisis['L'].rolling(3, min_periods=2).std()
    df_crisis['L_vol_6m'] = df_crisis['L'].rolling(6, min_periods=3).std()
    
    # Return features
    df_crisis['abs_ret'] = df_crisis['mkt_excess'].abs()
    df_crisis['cum_ret_3m'] = df_crisis['mkt_excess'].rolling(3, min_periods=2).sum()
    df_crisis['cum_ret_6m'] = df_crisis['mkt_excess'].rolling(6, min_periods=3).sum()
    
    # Drawdown
    cumret = (1 + df_crisis['mkt_excess']/100).cumprod()
    running_max = cumret.expanding().max()
    df_crisis['drawdown'] = (cumret / running_max - 1) * 100
    df_crisis['max_dd_6m'] = df_crisis['drawdown'].rolling(6, min_periods=3).min()
    
    # Crisis features list
    crisis_features = ['ret_vol_6m', 'ret_vol_12m', 'L_vol_3m', 'L_vol_6m', 
                       'abs_ret', 'cum_ret_3m', 'cum_ret_6m', 'max_dd_6m']
    
    # Add VIX if available
    if 'VIX' in df_crisis.columns:
        df_crisis['VIX_level'] = df_crisis['VIX']
        df_crisis['VIX_zscore'] = (df_crisis['VIX'] - df_crisis['VIX'].rolling(60, min_periods=30).mean()) / \
                                   df_crisis['VIX'].rolling(60, min_periods=30).std()
        crisis_features.extend(['VIX_level', 'VIX_zscore'])
    
    print(f"\n📊 Crisis Features Created: {len(crisis_features)}")
    
    # Extract and clean crisis data
    crisis_data = df_crisis[crisis_features].copy()
    crisis_data = crisis_data.fillna(method='bfill', limit=12)
    crisis_data = crisis_data.dropna()
    
    print(f"   Valid crisis detection samples: {len(crisis_data)}/{len(df_crisis)}")
    
    if len(crisis_data) >= 20:
        # Fit Isolation Forest
        iso_forest = IsolationForest(
            contamination=0.05,  # Expect ~5% crisis periods
            random_state=42,
            n_estimators=100,
            max_samples='auto'
        )
        
        outlier_labels = iso_forest.fit_predict(crisis_data.values)
        
        # Crisis = outlier + negative return (< -2%)
        is_crisis_candidate = (outlier_labels == -1)
        is_negative_return = (df_crisis.loc[crisis_data.index, 'mkt_excess'] < -2)
        
        df_crisis['is_crisis'] = False
        df_crisis.loc[crisis_data.index, 'is_crisis'] = is_crisis_candidate & is_negative_return
        
        n_crisis = df_crisis['is_crisis'].sum()
        crisis_pct = n_crisis / len(df_crisis) * 100
        
        print(f"\n✅ Crisis Detection Complete:")
        print(f"   Crisis periods:     {n_crisis} months ({crisis_pct:.1f}%)")
        print(f"   Non-crisis periods: {len(df_crisis) - n_crisis} months ({100-crisis_pct:.1f}%)")
        
        # Show detected crisis periods
        if n_crisis > 0:
            crisis_periods = df_crisis[df_crisis['is_crisis']].index
            print(f"\n📅 Detected Crisis Months:")
            for i, date in enumerate(crisis_periods[:10], 1):
                ret = df_crisis.loc[date, 'mkt_excess']
                L_val = df_crisis.loc[date, 'L']
                print(f"   {i:2d}. {date.strftime('%Y-%m')}: Return = {ret:+.2f}%, L = {L_val:+.2f}")
            if len(crisis_periods) > 10:
                print(f"   ... and {len(crisis_periods) - 10} more")
    else:
        print(f"\n⚠️  Insufficient data for crisis detection, marking all as non-crisis")
        df_crisis['is_crisis'] = False
    
    # Filter to non-crisis periods for HML regime analysis
    df_noncrisis = df_crisis[~df_crisis['is_crisis']].copy()
    
    print(f"\n📊 Non-Crisis Data: {len(df_noncrisis)} observations")
    
    # ==========================================
    # STEP 2: Define HML Regimes (Non-Crisis Only)
    # ==========================================
    
    print(f"\n{'='*80}")
    print("STEP 2: DEFINE HML REGIMES (NON-CRISIS PERIODS)")
    print(f"{'='*80}")
    
    # Check for HML data
    if 'hml' not in df_noncrisis.columns and 'HML' not in df_noncrisis.columns:
        raise ValueError("No HML column found in dataframe")
    
    hml_col = 'hml' if 'hml' in df_noncrisis.columns else 'HML'
    
    # Use 12-month rolling HML for smoother regimes
    if 'HML_12m' in df_noncrisis.columns:
        hml_metric = df_noncrisis['HML_12m']
        metric_col = 'HML_12m'
        metric_name = "12-month rolling HML"
    else:
        df_noncrisis['HML_12m'] = df_noncrisis[hml_col].rolling(12, min_periods=6).sum()
        hml_metric = df_noncrisis['HML_12m']
        metric_col = 'HML_12m'
        metric_name = "12-month rolling HML (created)"
    
    print(f"\n📊 Using: {metric_name}")
    print(f"   Range: {hml_metric.min():.2f} to {hml_metric.max():.2f}")
    print(f"   Mean:  {hml_metric.mean():.2f}")
    print(f"   Std:   {hml_metric.std():.2f}")
    
    # Define regime thresholds (tertiles on NON-CRISIS data)
    q33 = hml_metric.quantile(0.33)
    q67 = hml_metric.quantile(0.67)
    
    print(f"\n📊 HML Regime Thresholds (Non-Crisis):")
    print(f"   33rd percentile: {q33:.2f}")
    print(f"   67th percentile: {q67:.2f}")
    
    # Create regime labels
    df_noncrisis['hml_regime'] = pd.cut(
        hml_metric,
        bins=[-np.inf, q33, q67, np.inf],
        labels=[0, 1, 2]
    ).astype(int)
    
    df_noncrisis['hml_regime_label'] = df_noncrisis['hml_regime'].map({
        0: 'Growth-Friendly',
        1: 'Neutral', 
        2: 'Value-Friendly'
    })
    
    # Remove NaN regimes
    df_regimes = df_noncrisis.dropna(subset=['hml_regime'])
    
    print(f"\n📊 HML Regime Distribution (Non-Crisis):")
    for regime, label in [(0, 'Growth-Friendly'), (1, 'Neutral'), (2, 'Value-Friendly')]:
        mask = df_regimes['hml_regime'] == regime
        count = mask.sum()
        pct = count / len(df_regimes) * 100
        mean_hml = df_regimes.loc[mask, metric_col].mean()
        mean_L = df_regimes.loc[mask, 'L'].mean()
        print(f"   {label:16s}: {count:4d} ({pct:5.1f}%)")
        print(f"      Mean HML: {mean_hml:+.2f}, Mean L: {mean_L:+.2f}")
    
    # ==========================================
    # STEP 3: Feature Engineering (Liquidity)
    # ==========================================
    
    print(f"\n{'='*80}")
    print("STEP 3: LIQUIDITY FEATURES")
    print(f"{'='*80}")
    
    feature_cols = []
    
    # Core liquidity components
    liq_features = ['dlog_M2', 'dlog_FED_BS', 'EM', 'EB', 'EL_3y']
    feature_cols.extend([f for f in liq_features if f in df_regimes.columns])
    
    # Rates and spreads
    rate_features = ['Term_Spread', 'Credit_Spread', 'Real_Rate']
    feature_cols.extend([f for f in rate_features if f in df_regimes.columns])
    
    # Market indicators
    if 'VIX' in df_regimes.columns:
        feature_cols.append('VIX')
    
    # Growth rates
    if 'M2_growth' in df_regimes.columns:
        feature_cols.append('M2_growth')
    if 'FED_BS_growth' in df_regimes.columns:
        feature_cols.append('FED_BS_growth')
    
    # LIQUIDITY INDEX (L)
    feature_cols.append('L')
    
    # Technical on L
    df_regimes['L_ma_3'] = df_regimes['L'].rolling(3, min_periods=1).mean()
    df_regimes['L_ma_6'] = df_regimes['L'].rolling(6, min_periods=3).mean()
    df_regimes['L_momentum'] = df_regimes['L'].diff(3)
    df_regimes['L_accel'] = df_regimes['L_momentum'].diff(3)
    
    feature_cols.extend(['L_ma_3', 'L_ma_6', 'L_momentum', 'L_accel'])
    
    # Lagged HML (control for momentum)
    df_regimes['HML_lag3'] = df_regimes[hml_col].shift(3)
    df_regimes['HML_lag6'] = df_regimes[hml_col].shift(6)
    feature_cols.extend(['HML_lag3', 'HML_lag6'])
    
    print(f"\n✅ Selected {len(feature_cols)} features")
    
    # ==========================================
    # STEP 4: Prepare ML Dataset
    # ==========================================
    
    X_raw = df_regimes[feature_cols].copy()
    y = df_regimes['hml_regime'].copy()
    
    # Imputation
    X_imputed = X_raw.fillna(method='ffill', limit=6).fillna(method='bfill', limit=6)
    
    for col in X_imputed.columns:
        if X_imputed[col].isna().any():
            X_imputed[col] = X_imputed[col].fillna(X_imputed[col].median())
    
    # Final dataset
    valid_idx = X_imputed.notna().all(axis=1) & y.notna()
    X = X_imputed[valid_idx]
    y = y[valid_idx]
    
    print(f"\n📊 ML Dataset: {len(X)} samples × {len(feature_cols)} features")
    
    if len(X) < 50:
        raise ValueError(f"Insufficient data: {len(X)} samples")
    
    # ==========================================
    # STEP 5: Train Models with CV
    # ==========================================
    
    print(f"\n{'='*80}")
    print("STEP 4: MODEL TRAINING & CROSS-VALIDATION")
    print(f"{'='*80}")
    
    tscv = TimeSeriesSplit(n_splits=5)
    
    models = {
        'Random Forest': RandomForestClassifier(
            n_estimators=200,
            max_depth=8,
            min_samples_split=10,
            min_samples_leaf=5,
            class_weight='balanced',
            random_state=42,
            n_jobs=-1
        ),
        'XGBoost': XGBClassifier(
            n_estimators=200,
            max_depth=6,
            learning_rate=0.1,
            subsample=0.8,
            colsample_bytree=0.8,
            random_state=42,
            eval_metric='mlogloss',
            use_label_encoder=False,
            n_jobs=-1
        )
    }
    
    results = {}
    
    for name, model in models.items():
        print(f"\n{'='*70}")
        print(f"{name}")
        print(f"{'='*70}")
        
        # Cross-validation
        cv_scores = cross_val_score(model, X, y, cv=tscv, scoring='accuracy', n_jobs=-1)
        
        print(f"\n📊 Cross-Validation (Time Series):")
        print(f"   Mean Accuracy: {cv_scores.mean():.3f}{cv_scores.std():.3f})")
        print(f"   Fold Scores:   {[f'{s:.3f}' for s in cv_scores]}")
        
        baseline = 1/3
        improvement = cv_scores.mean() - baseline
        print(f"   vs Baseline:   {improvement:+.1%}")
        
        # Train full model
        model.fit(X, y)
        
        # In-sample
        y_pred = model.predict(X)
        train_acc = (y_pred == y).mean()
        
        print(f"\n📈 In-Sample Accuracy: {train_acc:.3f}")
        
        # Classification report
        print(f"\n📋 Classification Report:")
        print(classification_report(
            y, y_pred,
            target_names=['Growth-Friendly', 'Neutral', 'Value-Friendly'],
            zero_division=0
        ))
        
        # Confusion matrix
        cm = confusion_matrix(y, y_pred)
        cm_pct = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis] * 100
        
        print(f"\n🎯 Confusion Matrix (%):")
        cm_df = pd.DataFrame(
            cm_pct,
            index=['True_Growth', 'True_Neutral', 'True_Value'],
            columns=['Pred_Growth', 'Pred_Neutral', 'Pred_Value']
        )
        print(cm_df.round(1))
        
        results[name] = {
            'model': model,
            'cv_scores': cv_scores,
            'cv_mean': cv_scores.mean(),
            'cv_std': cv_scores.std(),
            'train_acc': train_acc,
            'y_pred': y_pred
        }
    
    # ==========================================
    # STEP 6: Feature Importance
    # ==========================================
    
    print(f"\n{'='*80}")
    print("STEP 5: FEATURE IMPORTANCE")
    print(f"{'='*80}")
    
    best_model_name = max(results, key=lambda k: results[k]['cv_mean'])
    best_model = results[best_model_name]['model']
    
    print(f"\n🏆 Best Model: {best_model_name} (CV: {results[best_model_name]['cv_mean']:.3f})")
    
    feature_importance = pd.DataFrame({
        'feature': X.columns,
        'importance': best_model.feature_importances_
    }).sort_values('importance', ascending=False)
    
    print(f"\n📊 Top 15 Features:")
    for i, row in feature_importance.head(15).iterrows():
        is_liquidity = any(liq in row['feature'] for liq in ['L', 'M2', 'FED', 'EM', 'EB', 'EL'])
        marker = "💧" if is_liquidity else "  "
        print(f"   {marker} {row['feature']:20s}: {row['importance']:.4f}")
    
    # Liquidity importance
    liquidity_features = [f for f in feature_importance['feature'] 
                         if any(liq in f for liq in ['L', 'M2', 'FED', 'EM', 'EB', 'EL'])]
    liquidity_importance = feature_importance[feature_importance['feature'].isin(liquidity_features)]['importance'].sum()
    
    print(f"\n💧 Total Liquidity Importance: {liquidity_importance:.1%}")
    
    # ==========================================
    # STEP 7: Hypothesis Test
    # ==========================================
    
    print(f"\n{'='*80}")
    print("HYPOTHESIS TEST: LIQUIDITY → HML REGIMES")
    print(f"{'='*80}")
    
    best_cv_acc = results[best_model_name]['cv_mean']
    baseline = 1/3
    improvement = best_cv_acc - baseline
    
    print(f"\n📊 Results:")
    print(f"   Baseline:      {baseline:.1%}")
    print(f"   CV Accuracy:   {best_cv_acc:.1%}")
    print(f"   Improvement:   {improvement:+.1%}")
    print(f"   Lift:          {best_cv_acc/baseline:.2f}x")
    
    if best_cv_acc >= 0.50:
        print(f"\n   ✅ STRONG SUPPORT for hypothesis")
    elif best_cv_acc >= 0.42:
        print(f"\n   ✔️  MODERATE SUPPORT for hypothesis")
    else:
        print(f"\n   ❌ WEAK SUPPORT for hypothesis")
    
    # ==========================================
    # STEP 8: Final Assignment
    # ==========================================
    
    df_final = df_crisis.copy()
    df_final['final_regime'] = np.nan
    df_final['final_regime_label'] = 'Unknown'
    
    # Crisis
    df_final.loc[df_final['is_crisis'], 'final_regime'] = 3
    df_final.loc[df_final['is_crisis'], 'final_regime_label'] = 'Crisis'
    
    # HML regimes (predicted)
    df_final.loc[X.index, 'final_regime'] = best_model.predict(X)
    df_final.loc[X.index, 'final_regime_label'] = pd.Series(
        best_model.predict(X),
        index=X.index
    ).map({0: 'Growth-Friendly', 1: 'Neutral', 2: 'Value-Friendly'})
    
    # Also keep actual HML regime
    df_final.loc[df_regimes.index, 'actual_hml_regime'] = df_regimes['hml_regime']
    df_final.loc[df_regimes.index, 'actual_hml_regime_label'] = df_regimes['hml_regime_label']
    
    print(f"\n📊 Final Distribution:")
    dist = df_final['final_regime_label'].value_counts()
    for regime in ['Crisis', 'Growth-Friendly', 'Neutral', 'Value-Friendly', 'Unknown']:
        if regime in dist.index:
            count = dist[regime]
            pct = count / len(df_final) * 100
            print(f"   {regime:16s}: {count:4d} ({pct:5.1f}%)")
    
    print(f"\n{'='*80}")
    print("✅ COMPLETE")
    print(f"{'='*80}")
    
    return df_final, results, feature_importance


# Run with crisis detection
df_hml_regimes, hml_results, hml_feat_imp = build_hml_regime_classification_with_crisis(
    df_features_aug,
    L_index
)

================================================================================
HML REGIME CLASSIFICATION (WITH CRISIS DETECTION)
================================================================================

✅ Aligned: 185 observations

================================================================================
STEP 1: CRISIS DETECTION
================================================================================

📊 Crisis Features Created: 10
   Valid crisis detection samples: 168/185

✅ Crisis Detection Complete:
   Crisis periods:     0 months (0.0%)
   Non-crisis periods: 185 months (100.0%)

📊 Non-Crisis Data: 185 observations

================================================================================
STEP 2: DEFINE HML REGIMES (NON-CRISIS PERIODS)
================================================================================

📊 Using: 12-month rolling HML
   Range: -0.42 to 0.31
   Mean:  -0.00
   Std:   0.13

📊 HML Regime Thresholds (Non-Crisis):
   33rd percentile: -0.07
   67th percentile: 0.06

📊 HML Regime Distribution (Non-Crisis):
   Growth-Friendly :   61 ( 33.0%)
      Mean HML: -0.14, Mean L: +0.16
   Neutral         :   63 ( 34.1%)
      Mean HML: -0.00, Mean L: -0.27
   Value-Friendly  :   61 ( 33.0%)
      Mean HML: +0.13, Mean L: +0.12

================================================================================
STEP 3: LIQUIDITY FEATURES
================================================================================

✅ Selected 18 features

📊 ML Dataset: 185 samples × 18 features

================================================================================
STEP 4: MODEL TRAINING & CROSS-VALIDATION
================================================================================

======================================================================
Random Forest
======================================================================

📊 Cross-Validation (Time Series):
   Mean Accuracy: 0.293 (±0.068)
   Fold Scores:   ['0.300', '0.300', '0.167', '0.367', '0.333']
   vs Baseline:   -4.0%

📈 In-Sample Accuracy: 0.924

📋 Classification Report:
                 precision    recall  f1-score   support

Growth-Friendly       0.94      0.97      0.95        61
        Neutral       0.92      0.89      0.90        63
 Value-Friendly       0.92      0.92      0.92        61

       accuracy                           0.92       185
      macro avg       0.92      0.92      0.92       185
   weighted avg       0.92      0.92      0.92       185


🎯 Confusion Matrix (%):
              Pred_Growth  Pred_Neutral  Pred_Value
True_Growth          96.7           3.3         0.0
True_Neutral          3.2          88.9         7.9
True_Value            3.3           4.9        91.8

======================================================================
XGBoost
======================================================================

📊 Cross-Validation (Time Series):
   Mean Accuracy: 0.293 (±0.108)
   Fold Scores:   ['0.433', '0.167', '0.167', '0.333', '0.367']
   vs Baseline:   -4.0%

📈 In-Sample Accuracy: 1.000

📋 Classification Report:
                 precision    recall  f1-score   support

Growth-Friendly       1.00      1.00      1.00        61
        Neutral       1.00      1.00      1.00        63
 Value-Friendly       1.00      1.00      1.00        61

       accuracy                           1.00       185
      macro avg       1.00      1.00      1.00       185
   weighted avg       1.00      1.00      1.00       185


🎯 Confusion Matrix (%):
              Pred_Growth  Pred_Neutral  Pred_Value
True_Growth         100.0           0.0         0.0
True_Neutral          0.0         100.0         0.0
True_Value            0.0           0.0       100.0

================================================================================
STEP 5: FEATURE IMPORTANCE
================================================================================

🏆 Best Model: Random Forest (CV: 0.293)

📊 Top 15 Features:
   💧 EB                  : 0.1444
   💧 EM                  : 0.1342
      Credit_Spread       : 0.0739
      VIX                 : 0.0566
   💧 FED_BS_growth       : 0.0543
   💧 L_momentum          : 0.0534
      Term_Spread         : 0.0515
   💧 EL_3y               : 0.0507
   💧 dlog_FED_BS         : 0.0505
   💧 L_ma_6              : 0.0471
   💧 L_accel             : 0.0464
   💧 HML_lag6            : 0.0436
   💧 dlog_M2             : 0.0347
   💧 L_ma_3              : 0.0329
   💧 HML_lag3            : 0.0327

💧 Total Liquidity Importance: 78.5%

================================================================================
HYPOTHESIS TEST: LIQUIDITY → HML REGIMES
================================================================================

📊 Results:
   Baseline:      33.3%
   CV Accuracy:   29.3%
   Improvement:   -4.0%
   Lift:          0.88x

   ❌ WEAK SUPPORT for hypothesis

📊 Final Distribution:
   Growth-Friendly :   63 ( 34.1%)
   Neutral         :   61 ( 33.0%)
   Value-Friendly  :   61 ( 33.0%)

================================================================================
✅ COMPLETE
================================================================================
factor_returns_by_regime = (
    df_regimes_fixed.groupby("regime_label")[["mkt_excess", "hml", "rmw", "cma"]]
    .mean() * 100
)
print(factor_returns_by_regime)
              mkt_excess       hml       rmw       cma
regime_label                                          
High            0.692131  0.374426  0.289180  0.525246
Neutral         0.826190  0.099524  0.430159 -0.144762
Tight           1.039672 -0.323115  0.007049 -0.536230
def build_hml_regime_classification_with_crisis(df, L_series):
    """
    1. Detect crisis periods (Isolation Forest)
    2. Classify NON-CRISIS periods by HML regimes
    3. Predict HML regimes from liquidity features
    """
    
    import numpy as np
    import pandas as pd
    from sklearn.ensemble import IsolationForest, RandomForestClassifier
    from xgboost import XGBClassifier
    from sklearn.model_selection import TimeSeriesSplit, cross_val_score
    from sklearn.metrics import classification_report, confusion_matrix
    
    print(f"\n{'='*80}")
    print("HML REGIME CLASSIFICATION (WITH CRISIS DETECTION)")
    print(f"{'='*80}")
    
    # ==========================================
    # STEP 0: Data Alignment
    # ==========================================
    
    common_idx = df.index.intersection(L_series.index)
    df_aligned = df.loc[common_idx].copy()
    df_aligned['L'] = L_series.loc[common_idx]
    
    print(f"\n✅ Aligned: {len(df_aligned)} observations")
    
    # ==========================================
    # STEP 1: Crisis Detection (Isolation Forest)
    # ==========================================
    
    print(f"\n{'='*80}")
    print("STEP 1: CRISIS DETECTION")
    print(f"{'='*80}")
    
    df_crisis = df_aligned.copy()
    
    # Create crisis features (volatility, drawdowns, tail events)
    df_crisis['ret_vol_6m'] = df_crisis['mkt_excess'].rolling(6, min_periods=3).std()
    df_crisis['ret_vol_12m'] = df_crisis['mkt_excess'].rolling(12, min_periods=6).std()
    df_crisis['L_vol_3m'] = df_crisis['L'].rolling(3, min_periods=2).std()
    df_crisis['L_vol_6m'] = df_crisis['L'].rolling(6, min_periods=3).std()
    
    # Return features
    df_crisis['abs_ret'] = df_crisis['mkt_excess'].abs()
    df_crisis['cum_ret_3m'] = df_crisis['mkt_excess'].rolling(3, min_periods=2).sum()
    df_crisis['cum_ret_6m'] = df_crisis['mkt_excess'].rolling(6, min_periods=3).sum()
    
    # Drawdown
    cumret = (1 + df_crisis['mkt_excess']/100).cumprod()
    running_max = cumret.expanding().max()
    df_crisis['drawdown'] = (cumret / running_max - 1) * 100
    df_crisis['max_dd_6m'] = df_crisis['drawdown'].rolling(6, min_periods=3).min()
    
    # Crisis features list
    crisis_features = ['ret_vol_6m', 'ret_vol_12m', 'L_vol_3m', 'L_vol_6m', 
                       'abs_ret', 'cum_ret_3m', 'cum_ret_6m', 'max_dd_6m']
    
    # Add VIX if available
    if 'VIX' in df_crisis.columns:
        df_crisis['VIX_level'] = df_crisis['VIX']
        df_crisis['VIX_zscore'] = (df_crisis['VIX'] - df_crisis['VIX'].rolling(60, min_periods=30).mean()) / \
                                   df_crisis['VIX'].rolling(60, min_periods=30).std()
        crisis_features.extend(['VIX_level', 'VIX_zscore'])
    
    print(f"\n📊 Crisis Features Created: {len(crisis_features)}")
    
    # Extract and clean crisis data
    crisis_data = df_crisis[crisis_features].copy()
    crisis_data = crisis_data.fillna(method='bfill', limit=12)
    crisis_data = crisis_data.dropna()
    
    print(f"   Valid crisis detection samples: {len(crisis_data)}/{len(df_crisis)}")
    
    if len(crisis_data) >= 20:
        # Fit Isolation Forest
        iso_forest = IsolationForest(
            contamination=0.05,  # Expect ~5% crisis periods
            random_state=42,
            n_estimators=100,
            max_samples='auto'
        )
        
        outlier_labels = iso_forest.fit_predict(crisis_data.values)
        
        # Crisis = outlier + negative return (< -2%)
        is_crisis_candidate = (outlier_labels == -1)
        is_negative_return = (df_crisis.loc[crisis_data.index, 'mkt_excess'] < -2)
        
        df_crisis['is_crisis'] = False
        df_crisis.loc[crisis_data.index, 'is_crisis'] = is_crisis_candidate
        
        n_crisis = df_crisis['is_crisis'].sum()
        crisis_pct = n_crisis / len(df_crisis) * 100
        
        print(f"\n✅ Crisis Detection Complete:")
        print(f"   Crisis periods:     {n_crisis} months ({crisis_pct:.1f}%)")
        print(f"   Non-crisis periods: {len(df_crisis) - n_crisis} months ({100-crisis_pct:.1f}%)")
        
        # Show detected crisis periods
        if n_crisis > 0:
            crisis_periods = df_crisis[df_crisis['is_crisis']].index
            print(f"\n📅 Detected Crisis Months:")
            for i, date in enumerate(crisis_periods[:10], 1):
                ret = df_crisis.loc[date, 'mkt_excess']
                L_val = df_crisis.loc[date, 'L']
                print(f"   {i:2d}. {date.strftime('%Y-%m')}: Return = {ret:+.2f}%, L = {L_val:+.2f}")
            if len(crisis_periods) > 10:
                print(f"   ... and {len(crisis_periods) - 10} more")
    else:
        print(f"\n⚠️  Insufficient data for crisis detection, marking all as non-crisis")
        df_crisis['is_crisis'] = False
    
    # Filter to non-crisis periods for HML regime analysis
    df_noncrisis = df_crisis[~df_crisis['is_crisis']].copy()
    
    print(f"\n📊 Non-Crisis Data: {len(df_noncrisis)} observations")
    
    # ==========================================
    # STEP 2: Define HML Regimes (Non-Crisis Only)
    # ==========================================
    
    print(f"\n{'='*80}")
    print("STEP 2: DEFINE HML REGIMES (NON-CRISIS PERIODS)")
    print(f"{'='*80}")
    
    # Check for HML data
    if 'hml' not in df_noncrisis.columns and 'HML' not in df_noncrisis.columns:
        raise ValueError("No HML column found in dataframe")
    
    hml_col = 'hml' if 'hml' in df_noncrisis.columns else 'HML'
    
    # Use 12-month rolling HML for smoother regimes
    # if 'HML_12m' in df_noncrisis.columns:
    #     hml_metric = df_noncrisis['HML_12m']
    #     metric_col = 'HML_12m'
    #     metric_name = "12-month rolling HML"
    # else:
    #     df_noncrisis['HML_12m'] = df_noncrisis[hml_col].rolling(12, min_periods=6).sum()
    #     hml_metric = df_noncrisis['HML_12m']
    #     metric_col = 'HML_12m'
    #     metric_name = "12-month rolling HML (created)"
    df_noncrisis['HML_forward_6m'] = df_noncrisis[hml_col].shift(-6).rolling(
        6, min_periods=3
    ).sum()
    hml_metric = df_noncrisis['HML_forward_6m']
    metric_col = 'HML_forward_6m'
    metric_name = "6-month FORWARD HML"
    
    print(f"\n📊 Using: {metric_name}")
    print(f"   Range: {hml_metric.min():.2f} to {hml_metric.max():.2f}")
    print(f"   Mean:  {hml_metric.mean():.2f}")
    print(f"   Std:   {hml_metric.std():.2f}")

    # Remove NaN BEFORE creating regimes
    df_noncrisis_valid = df_noncrisis.dropna(subset=['HML_forward_6m']).copy()
    hml_metric_valid = df_noncrisis_valid['HML_forward_6m']
    
    print(f"\n📊 After removing NaN: {len(df_noncrisis_valid)} observations")
    
    # Define regime thresholds (tertiles on NON-CRISIS data)
    q33 = hml_metric.quantile(0.33)
    q67 = hml_metric.quantile(0.67)
    
    print(f"\n📊 HML Regime Thresholds (Non-Crisis):")
    print(f"   33rd percentile: {q33:.2f}")
    print(f"   67th percentile: {q67:.2f}")
    
    # Create regime labels
    df_noncrisis['hml_regime'] = pd.cut(
        hml_metric_valid,
        bins=[-np.inf, q33, q67, np.inf],
        labels=[0, 1, 2]
    ).astype(int)
    
    df_noncrisis['hml_regime_label'] = df_noncrisis['hml_regime'].map({
        0: 'Growth-Friendly',
        1: 'Neutral', 
        2: 'Value-Friendly'
    })
    
    # Remove NaN regimes
    df_regimes = df_noncrisis.dropna(subset=['hml_regime'])
    
    print(f"\n📊 HML Regime Distribution (Non-Crisis):")
    for regime, label in [(0, 'Growth-Friendly'), (1, 'Neutral'), (2, 'Value-Friendly')]:
        mask = df_regimes['hml_regime'] == regime
        count = mask.sum()
        pct = count / len(df_regimes) * 100
        mean_hml = df_regimes.loc[mask, metric_col].mean()
        mean_L = df_regimes.loc[mask, 'L'].mean()
        print(f"   {label:16s}: {count:4d} ({pct:5.1f}%)")
        print(f"      Mean HML: {mean_hml:+.2f}, Mean L: {mean_L:+.2f}")
    
    # ==========================================
    # STEP 3: Feature Engineering (Liquidity)
    # ==========================================
    
    print(f"\n{'='*80}")
    print("STEP 3: LIQUIDITY FEATURES")
    print(f"{'='*80}")
    
    feature_cols = []
    
    # Core liquidity components
    liq_features = ['dlog_M2', 'dlog_FED_BS', 'EM', 'EB', 'EL_3y']
    feature_cols.extend([f for f in liq_features if f in df_regimes.columns])
    
    # Rates and spreads
    rate_features = ['Term_Spread', 'Credit_Spread', 'Real_Rate']
    feature_cols.extend([f for f in rate_features if f in df_regimes.columns])
    
    # Market indicators
    if 'VIX' in df_regimes.columns:
        feature_cols.append('VIX')
    
    # Growth rates
    if 'M2_growth' in df_regimes.columns:
        feature_cols.append('M2_growth')
    if 'FED_BS_growth' in df_regimes.columns:
        feature_cols.append('FED_BS_growth')
    
    # LIQUIDITY INDEX (L)
    feature_cols.append('L')
    
    # Technical on L
    df_regimes['L_ma_3'] = df_regimes['L'].rolling(3, min_periods=1).mean()
    df_regimes['L_ma_6'] = df_regimes['L'].rolling(6, min_periods=3).mean()
    df_regimes['L_momentum'] = df_regimes['L'].diff(3)
    df_regimes['L_accel'] = df_regimes['L_momentum'].diff(3)
    
    feature_cols.extend(['L_ma_3', 'L_ma_6', 'L_momentum', 'L_accel'])
    
    # Lagged HML (control for momentum)
    df_regimes['HML_lag3'] = df_regimes[hml_col].shift(3)
    df_regimes['HML_lag6'] = df_regimes[hml_col].shift(6)
    feature_cols.extend(['HML_lag3', 'HML_lag6'])
    
    print(f"\n✅ Selected {len(feature_cols)} features")
    
    # ==========================================
    # STEP 4: Prepare ML Dataset
    # ==========================================
    
    X_raw = df_regimes[feature_cols].copy()
    y = df_regimes['hml_regime'].copy()
    
    # Imputation
    X_imputed = X_raw.fillna(method='ffill', limit=6).fillna(method='bfill', limit=6)
    
    for col in X_imputed.columns:
        if X_imputed[col].isna().any():
            X_imputed[col] = X_imputed[col].fillna(X_imputed[col].median())
    
    # Final dataset
    valid_idx = X_imputed.notna().all(axis=1) & y.notna()
    X = X_imputed[valid_idx]
    y = y[valid_idx]
    
    print(f"\n📊 ML Dataset: {len(X)} samples × {len(feature_cols)} features")
    
    if len(X) < 50:
        raise ValueError(f"Insufficient data: {len(X)} samples")
    
    # ==========================================
    # STEP 5: Train Models with CV
    # ==========================================
    
    print(f"\n{'='*80}")
    print("STEP 4: MODEL TRAINING & CROSS-VALIDATION")
    print(f"{'='*80}")
    
    tscv = TimeSeriesSplit(n_splits=5)
    
    models = {
        'Random Forest': RandomForestClassifier(
            n_estimators=200,
            max_depth=8,
            min_samples_split=10,
            min_samples_leaf=5,
            class_weight='balanced',
            random_state=42,
            n_jobs=-1
        ),
        'XGBoost': XGBClassifier(
            n_estimators=200,
            max_depth=6,
            learning_rate=0.1,
            subsample=0.8,
            colsample_bytree=0.8,
            random_state=42,
            eval_metric='mlogloss',
            use_label_encoder=False,
            n_jobs=-1
        )
    }
    
    results = {}
    
    for name, model in models.items():
        print(f"\n{'='*70}")
        print(f"{name}")
        print(f"{'='*70}")
        
        # Cross-validation
        cv_scores = cross_val_score(model, X, y, cv=tscv, scoring='accuracy', n_jobs=-1)
        
        print(f"\n📊 Cross-Validation (Time Series):")
        print(f"   Mean Accuracy: {cv_scores.mean():.3f}{cv_scores.std():.3f})")
        print(f"   Fold Scores:   {[f'{s:.3f}' for s in cv_scores]}")
        
        baseline = 1/3
        improvement = cv_scores.mean() - baseline
        print(f"   vs Baseline:   {improvement:+.1%}")
        
        # Train full model
        model.fit(X, y)
        
        # In-sample
        y_pred = model.predict(X)
        train_acc = (y_pred == y).mean()
        
        print(f"\n📈 In-Sample Accuracy: {train_acc:.3f}")
        
        # Classification report
        print(f"\n📋 Classification Report:")
        print(classification_report(
            y, y_pred,
            target_names=['Growth-Friendly', 'Neutral', 'Value-Friendly'],
            zero_division=0
        ))
        
        # Confusion matrix
        cm = confusion_matrix(y, y_pred)
        cm_pct = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis] * 100
        
        print(f"\n🎯 Confusion Matrix (%):")
        cm_df = pd.DataFrame(
            cm_pct,
            index=['True_Growth', 'True_Neutral', 'True_Value'],
            columns=['Pred_Growth', 'Pred_Neutral', 'Pred_Value']
        )
        print(cm_df.round(1))
        
        results[name] = {
            'model': model,
            'cv_scores': cv_scores,
            'cv_mean': cv_scores.mean(),
            'cv_std': cv_scores.std(),
            'train_acc': train_acc,
            'y_pred': y_pred
        }
    
    # ==========================================
    # STEP 6: Feature Importance
    # ==========================================
    
    print(f"\n{'='*80}")
    print("STEP 5: FEATURE IMPORTANCE")
    print(f"{'='*80}")
    
    best_model_name = max(results, key=lambda k: results[k]['cv_mean'])
    best_model = results[best_model_name]['model']
    
    print(f"\n🏆 Best Model: {best_model_name} (CV: {results[best_model_name]['cv_mean']:.3f})")
    
    feature_importance = pd.DataFrame({
        'feature': X.columns,
        'importance': best_model.feature_importances_
    }).sort_values('importance', ascending=False)
    
    print(f"\n📊 Top 15 Features:")
    for i, row in feature_importance.head(15).iterrows():
        is_liquidity = any(liq in row['feature'] for liq in ['L', 'M2', 'FED', 'EM', 'EB', 'EL'])
        marker = "💧" if is_liquidity else "  "
        print(f"   {marker} {row['feature']:20s}: {row['importance']:.4f}")
    
    # Liquidity importance
    liquidity_features = [f for f in feature_importance['feature'] 
                         if any(liq in f for liq in ['L', 'M2', 'FED', 'EM', 'EB', 'EL'])]
    liquidity_importance = feature_importance[feature_importance['feature'].isin(liquidity_features)]['importance'].sum()
    
    print(f"\n💧 Total Liquidity Importance: {liquidity_importance:.1%}")
    
    # ==========================================
    # STEP 7: Hypothesis Test
    # ==========================================
    
    print(f"\n{'='*80}")
    print("HYPOTHESIS TEST: LIQUIDITY → HML REGIMES")
    print(f"{'='*80}")
    
    best_cv_acc = results[best_model_name]['cv_mean']
    baseline = 1/3
    improvement = best_cv_acc - baseline
    
    print(f"\n📊 Results:")
    print(f"   Baseline:      {baseline:.1%}")
    print(f"   CV Accuracy:   {best_cv_acc:.1%}")
    print(f"   Improvement:   {improvement:+.1%}")
    print(f"   Lift:          {best_cv_acc/baseline:.2f}x")
    
    if best_cv_acc >= 0.50:
        print(f"\n   ✅ STRONG SUPPORT for hypothesis")
    elif best_cv_acc >= 0.42:
        print(f"\n   ✔️  MODERATE SUPPORT for hypothesis")
    else:
        print(f"\n   ❌ WEAK SUPPORT for hypothesis")
    
    # ==========================================
    # STEP 8: Final Assignment
    # ==========================================
    
    df_final = df_crisis.copy()
    df_final['final_regime'] = np.nan
    df_final['final_regime_label'] = 'Unknown'
    
    # Crisis
    df_final.loc[df_final['is_crisis'], 'final_regime'] = 3
    df_final.loc[df_final['is_crisis'], 'final_regime_label'] = 'Crisis'
    
    # HML regimes (predicted)
    df_final.loc[X.index, 'final_regime'] = best_model.predict(X)
    df_final.loc[X.index, 'final_regime_label'] = pd.Series(
        best_model.predict(X),
        index=X.index
    ).map({0: 'Growth-Friendly', 1: 'Neutral', 2: 'Value-Friendly'})
    
    # Also keep actual HML regime
    df_final.loc[df_regimes.index, 'actual_hml_regime'] = df_regimes['hml_regime']
    df_final.loc[df_regimes.index, 'actual_hml_regime_label'] = df_regimes['hml_regime_label']
    
    print(f"\n📊 Final Distribution:")
    dist = df_final['final_regime_label'].value_counts()
    for regime in ['Crisis', 'Growth-Friendly', 'Neutral', 'Value-Friendly', 'Unknown']:
        if regime in dist.index:
            count = dist[regime]
            pct = count / len(df_final) * 100
            print(f"   {regime:16s}: {count:4d} ({pct:5.1f}%)")
    
    print(f"\n{'='*80}")
    print("✅ COMPLETE")
    print(f"{'='*80}")
    
    return df_final, results, feature_importance


# Run with crisis detection
df_hml_regimes, hml_results, hml_feat_imp = build_hml_regime_classification_with_crisis(
    df_features_aug,
    L_index
)

================================================================================
HML REGIME CLASSIFICATION (WITH CRISIS DETECTION)
================================================================================

✅ Aligned: 185 observations

================================================================================
STEP 1: CRISIS DETECTION
================================================================================

📊 Crisis Features Created: 10
   Valid crisis detection samples: 168/185

✅ Crisis Detection Complete:
   Crisis periods:     9 months (4.9%)
   Non-crisis periods: 176 months (95.1%)

📅 Detected Crisis Months:
    1. 2008-09: Return = -0.09%, L = +2.15
    2. 2008-10: Return = -0.17%, L = +5.37
    3. 2008-12: Return = +0.02%, L = +4.80
    4. 2009-03: Return = +0.09%, L = +4.66
    5. 2009-04: Return = +0.10%, L = +4.80
    6. 2009-06: Return = +0.00%, L = +3.73
    7. 2009-07: Return = +0.08%, L = +3.40
    8. 2020-03: Return = -0.13%, L = +1.94
    9. 2020-04: Return = +0.14%, L = +2.11

📊 Non-Crisis Data: 176 observations

================================================================================
STEP 2: DEFINE HML REGIMES (NON-CRISIS PERIODS)
================================================================================

📊 Using: 6-month FORWARD HML
   Range: -0.20 to 0.25
   Mean:  0.00
   Std:   0.07

📊 After removing NaN: 171 observations

📊 HML Regime Thresholds (Non-Crisis):
   33rd percentile: -0.03
   67th percentile: 0.03

📊 HML Regime Distribution (Non-Crisis):
   Growth-Friendly :   57 ( 33.3%)
      Mean HML: -0.07, Mean L: -0.54
   Neutral         :   58 ( 33.9%)
      Mean HML: +0.00, Mean L: +0.03
   Value-Friendly  :   56 ( 32.7%)
      Mean HML: +0.08, Mean L: +0.03

================================================================================
STEP 3: LIQUIDITY FEATURES
================================================================================

✅ Selected 18 features

📊 ML Dataset: 171 samples × 18 features

================================================================================
STEP 4: MODEL TRAINING & CROSS-VALIDATION
================================================================================

======================================================================
Random Forest
======================================================================

📊 Cross-Validation (Time Series):
   Mean Accuracy: 0.307 (±0.107)
   Fold Scores:   ['0.286', '0.214', '0.179', '0.464', '0.393']
   vs Baseline:   -2.6%

📈 In-Sample Accuracy: 0.895

📋 Classification Report:
                 precision    recall  f1-score   support

Growth-Friendly       0.89      0.95      0.92        57
        Neutral       0.84      0.84      0.84        58
 Value-Friendly       0.96      0.89      0.93        56

       accuracy                           0.89       171
      macro avg       0.90      0.90      0.90       171
   weighted avg       0.90      0.89      0.89       171


🎯 Confusion Matrix (%):
              Pred_Growth  Pred_Neutral  Pred_Value
True_Growth          94.7           5.3         0.0
True_Neutral         12.1          84.5         3.4
True_Value            0.0          10.7        89.3

======================================================================
XGBoost
======================================================================

📊 Cross-Validation (Time Series):
   Mean Accuracy: 0.214 (±0.106)
   Fold Scores:   ['0.179', '0.071', '0.179', '0.250', '0.393']
   vs Baseline:   -11.9%

📈 In-Sample Accuracy: 1.000

📋 Classification Report:
                 precision    recall  f1-score   support

Growth-Friendly       1.00      1.00      1.00        57
        Neutral       1.00      1.00      1.00        58
 Value-Friendly       1.00      1.00      1.00        56

       accuracy                           1.00       171
      macro avg       1.00      1.00      1.00       171
   weighted avg       1.00      1.00      1.00       171


🎯 Confusion Matrix (%):
              Pred_Growth  Pred_Neutral  Pred_Value
True_Growth         100.0           0.0         0.0
True_Neutral          0.0         100.0         0.0
True_Value            0.0           0.0       100.0

================================================================================
STEP 5: FEATURE IMPORTANCE
================================================================================

🏆 Best Model: Random Forest (CV: 0.307)

📊 Top 15 Features:
   💧 FED_BS_growth       : 0.1034
      Term_Spread         : 0.0878
   💧 M2_growth           : 0.0815
   💧 EM                  : 0.0766
   💧 EB                  : 0.0755
   💧 HML_lag6            : 0.0583
   💧 L                   : 0.0580
      Credit_Spread       : 0.0547
   💧 EL_3y               : 0.0530
   💧 L_ma_3              : 0.0501
   💧 L_ma_6              : 0.0485
   💧 dlog_M2             : 0.0465
   💧 L_momentum          : 0.0414
      Real_Rate           : 0.0399
   💧 dlog_FED_BS         : 0.0364

💧 Total Liquidity Importance: 78.5%

================================================================================
HYPOTHESIS TEST: LIQUIDITY → HML REGIMES
================================================================================

📊 Results:
   Baseline:      33.3%
   CV Accuracy:   30.7%
   Improvement:   -2.6%
   Lift:          0.92x

   ❌ WEAK SUPPORT for hypothesis

📊 Final Distribution:
   Crisis          :    9 (  4.9%)
   Growth-Friendly :   61 ( 33.0%)
   Neutral         :   58 ( 31.4%)
   Value-Friendly  :   52 ( 28.1%)
   Unknown         :    5 (  2.7%)

================================================================================
✅ COMPLETE
================================================================================