אימות גורמים בתיק מנוהל באמצעות LASSO לאחר בחירה כפולה

# Factor Zoo Validation via Post-Double-Selection LASSO


# Factor Zoo Validation via Post-Double-Selection LASSO

**Chapter 15: Causal Machine Learning**

This notebook uses post-double-selection LASSO to ask whether managed-portfolio
factor returns add incremental time-series explanatory power for SPY after ten
principal-component controls. It is a compact factor-spanning application of the
Belloni-Chernozhukov-Hansen selection logic discussed in this chapter.

## Learning objectives

- distinguish a naive single-factor association from a conditional factor loading;
- implement both LASSO selections with time-ordered cross-validation;
- carry the selected union into an intercept-inclusive OLS regression with HAC inference;
- interpret what the exercise does and does not establish about the factor zoo.

**Book references**: Chapter 14, Section 14.1 (Making the case for latent factors), and
Chapter 15, Section 15.4 (Isolating factor effects with DML).

**Prerequisites**: [`01_pca_equity_sectors`](../14_latent_factors/01_pca_equity_sectors.ipynb)
and [`03_econml_dml`](03_econml_dml.ipynb).

## Setup

The ETF data are daily, so the cube-root Newey-West bandwidth provides a transparent
default for serial dependence. The bandwidth is computed after the 60-day warm-up.

```python
"""Factor-spanning validation with post-double-selection LASSO."""

import matplotlib.pyplot as plt
import numpy as np
import polars as pl
import statsmodels.api as sm
from sklearn.decomposition import PCA
from sklearn.linear_model import Lasso
from sklearn.model_selection import GridSearchCV, TimeSeriesSplit
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

from data import load_etfs
from utils.reproducibility import set_global_seeds
from utils.style import COLORS, FIGSIZE, add_message_title, show_with_alt, zero_line
```

```python
START_DATE = "2006-01-01"
END_DATE = "2024-12-31"
MIN_OBSERVATIONS = 252
N_PCA_FACTORS = 10
N_CV_SPLITS = 5
N_ALPHAS = 60
SIGNIFICANCE_LEVEL = 0.05
MAX_SYMBOLS = 0  # 0 = all eligible symbols
OUTCOME_SYMBOL = "SPY"
WARMUP = 60
SEED = 42
```

```python
set_global_seeds(SEED)
```

## 1. Build a balanced, disjoint ETF panel

SPY is excluded from the factor-building universe and retained only as the outcome.
The source contains ETFs with different inception dates. Replacing pre-inception
observations with zero would manufacture returns, so the analysis retains the
long-history symbols observed on every date in the requested sample.

This eligibility rule uses the current curated ETF list and is not a point-in-time
historical universe. The result is an in-sample teaching exercise, not a
survivorship-free backtest, and nothing here is held out.

```python
etf_data = load_etfs(start_date=START_DATE, end_date=END_DATE).sort(["symbol", "timestamp"])

duplicate_keys = etf_data.select(pl.struct("symbol", "timestamp").is_duplicated().sum()).item()
if duplicate_keys:
    raise ValueError(f"ETF input contains {duplicate_keys} duplicate symbol-timestamp keys")

etf_returns = etf_data.with_columns(
    pl.col("close").pct_change().over("symbol").alias("return")
).drop_nulls(subset=["return"])

n_source_dates = etf_returns["timestamp"].n_unique()
symbol_coverage = etf_returns.group_by("symbol").len().sort("symbol")
eligible_symbols = symbol_coverage.filter(
    (pl.col("len") == n_source_dates) & (pl.col("len") >= MIN_OBSERVATIONS)
)["symbol"].to_list()

if OUTCOME_SYMBOL not in eligible_symbols:
    raise ValueError(f"Outcome {OUTCOME_SYMBOL!r} lacks complete sample coverage")
if MAX_SYMBOLS > 0:
    eligible_symbols = eligible_symbols[:MAX_SYMBOLS]
    if OUTCOME_SYMBOL not in eligible_symbols:
        eligible_symbols.append(OUTCOME_SYMBOL)
```

Pivoting after the coverage filter produces a genuinely balanced panel. The null
and finite-value assertions make the no-imputation contract executable.

```python
return_wide = (
    etf_returns.filter(pl.col("symbol").is_in(eligible_symbols))
    .pivot(on="symbol", index="timestamp", values="return")
    .sort("timestamp")
)

symbols_all = [column for column in return_wide.columns if column != "timestamp"]
if len(symbols_all) <= N_PCA_FACTORS + 1:
    raise ValueError("The balanced universe is too small for the requested PCA basis")
if return_wide.select(pl.sum_horizontal(pl.exclude("timestamp").null_count())).item() != 0:
    raise ValueError("Balanced return panel contains missing values")

returns_all = return_wide.select(symbols_all).to_numpy().astype(np.float64)
if not np.isfinite(returns_all).all():
    raise ValueError("Balanced return panel contains non-finite values")

outcome_idx = symbols_all.index(OUTCOME_SYMBOL)
outcome_return = returns_all[:, outcome_idx]
zoo_mask = np.ones(returns_all.shape[1], dtype=bool)
zoo_mask[outcome_idx] = False
zoo_returns = returns_all[:, zoo_mask]
zoo_symbols = [symbol for symbol in symbols_all if symbol != OUTCOME_SYMBOL]
T, N = zoo_returns.shape

print(
    f"Analysis window: {return_wide['timestamp'].min()} to {return_wide['timestamp'].max()} "
    f"({T:,} trading days)"
)
print(f"Outcome: {OUTCOME_SYMBOL}; factor-building universe: {N} long-history ETFs")
print("Missing returns imputed: 0")
```

## 2. Extract principal-component controls

Standardization prevents high-volatility ETFs from dominating the covariance
structure. The final basis is fitted on the full post-warm-up inference sample
because these are explicitly in-sample controls; uncertainty is conditional on
this estimated basis. Selection below refits both transformations inside each fold.

```python
pca_scaler = StandardScaler()
zoo_returns_inference = zoo_returns[WARMUP:]
zoo_returns_scaled = pca_scaler.fit_transform(zoo_returns_inference)
pca = PCA(n_components=N_PCA_FACTORS, random_state=SEED)
factor_returns = pca.fit_transform(zoo_returns_scaled)

factor_names = [f"PC{i + 1}" for i in range(N_PCA_FACTORS)]
variance_table = pl.DataFrame(
    {
        "factor": factor_names,
        "variance_explained": pca.explained_variance_ratio_,
        "cumulative_variance": np.cumsum(pca.explained_variance_ratio_),
    }
).with_columns(pl.selectors.numeric().round(4))

print(
    f"Cumulative variance explained by {N_PCA_FACTORS} PCs: {pca.explained_variance_ratio_.sum():.1%}"
)
variance_table
```

## 3. Construct lagged managed-portfolio factors

Each portfolio uses information ending at $t-1$ to set long and short baskets,
then records their return at $t$. The four candidates differ only in the ranking
statistic and lookback window.

```python
momentum_20d = np.zeros(T)
for t in range(20, T):
    momentum = zoo_returns[t - 20 : t].sum(axis=0)
    top = momentum >= np.percentile(momentum, 80)
    bottom = momentum <= np.percentile(momentum, 20)
    momentum_20d[t] = zoo_returns[t, top].mean() - zoo_returns[t, bottom].mean()

momentum_60d = np.zeros(T)
for t in range(60, T):
    momentum = zoo_returns[t - 60 : t].sum(axis=0)
    top = momentum >= np.percentile(momentum, 80)
    bottom = momentum <= np.percentile(momentum, 20)
    momentum_60d[t] = zoo_returns[t, top].mean() - zoo_returns[t, bottom].mean()
```

Low-volatility and mean-reversion portfolios reverse the ranking direction: they
buy the low-volatility or recent-loser quintile and sell the opposite quintile.

```python
low_vol = np.zeros(T)
for t in range(60, T):
    volatility = zoo_returns[t - 60 : t].std(axis=0)
    low = volatility <= np.percentile(volatility, 20)
    high = volatility >= np.percentile(volatility, 80)
    low_vol[t] = zoo_returns[t, low].mean() - zoo_returns[t, high].mean()

mean_reversion = np.zeros(T)
for t in range(5, T):
    recent_return = zoo_returns[t - 5 : t].sum(axis=0)
    losers = recent_return <= np.percentile(recent_return, 20)
    winners = recent_return >= np.percentile(recent_return, 80)
    mean_reversion[t] = zoo_returns[t, losers].mean() - zoo_returns[t, winners].mean()
```

The common warm-up removes the initialized zeros before any inference. Annualized
Sharpe ratios here are descriptive summaries, not selection criteria.

```python
candidate_names = ["Mom_20d", "Mom_60d", "LowVol", "MeanRev"]
candidates = np.column_stack([momentum_20d, momentum_60d, low_vol, mean_reversion])[WARMUP:]
outcome_trimmed = outcome_return[WARMUP:]
controls_trimmed = factor_returns
HAC_LAGS = max(1, int(len(outcome_trimmed) ** (1 / 3)))

candidate_summary = pl.DataFrame(
    {
        "factor": candidate_names,
        "annualized_sharpe": [
            candidates[:, i].mean() / candidates[:, i].std(ddof=1) * np.sqrt(252)
            for i in range(len(candidate_names))
        ],
        "daily_volatility": [candidates[:, i].std(ddof=1) for i in range(len(candidate_names))],
    }
).with_columns(pl.selectors.numeric().round(4))

print(f"Inference sample: {len(outcome_trimmed):,} days; Newey-West bandwidth: {HAC_LAGS} lags")
candidate_summary
```

## 4. Match the naive and conditional estimands

A factor-mean test and a regression loading answer different questions. The
comparison below holds the target fixed: the naive slope comes from SPY on one
candidate, while the post-selection slope adds selected PCA controls. Both models
include an intercept and use the same Newey-West covariance estimator.

```python
def ols_hac_test(outcome: np.ndarray, design: np.ndarray, hac_lags: int) -> dict[str, float]:
    """Estimate the first slope in an intercept-inclusive OLS-HAC regression."""
    design_2d = design.reshape(-1, 1) if design.ndim == 1 else design
    model = sm.OLS(outcome, sm.add_constant(design_2d)).fit(
        cov_type="HAC", cov_kwds={"maxlags": hac_lags}
    )
    return {
        "coef": float(model.params[1]),
        "se": float(model.bse[1]),
        "t_stat": float(model.tvalues[1]),
        "p_value": float(model.pvalues[1]),
    }
```

```python
naive_results = []
for index, name in enumerate(candidate_names):
    result = ols_hac_test(outcome_trimmed, candidates[:, index], HAC_LAGS)
    naive_results.append({"factor": name, **result})

naive_table = pl.DataFrame(naive_results).with_columns(
    pl.col("coef", "se").round(4),
    pl.col("t_stat").round(2),
    pl.col("p_value").round(4),
)
naive_table
```

## 5. Select controls without leaking future folds

Each LASSO uses an expanding time-series split. Its pipeline refits the input
scaler, PCA basis, and control scaler on each training fold, so later observations
cannot change an earlier validation prediction. A fixed broad alpha grid avoids
using the full target series to calibrate the candidate penalties.

```python
def fit_lasso_selector(target: np.ndarray, pca_inputs: np.ndarray) -> dict:
    """Tune a fold-local PCA-LASSO pipeline and return its full-sample support."""
    alpha_grid = np.geomspace(1e-8, 1e-2, N_ALPHAS)
    pipeline = make_pipeline(
        StandardScaler(),
        PCA(n_components=N_PCA_FACTORS, random_state=SEED),
        StandardScaler(),
        Lasso(max_iter=20_000, random_state=SEED),
    )
    search = GridSearchCV(
        pipeline,
        {"lasso__alpha": alpha_grid},
        cv=TimeSeriesSplit(n_splits=N_CV_SPLITS),
        scoring="neg_mean_squared_error",
        # The grid is small and the design is a few thousand rows; a worker per core would
        # take the whole machine from every other notebook executing beside this one.
        n_jobs=1,
    )
    search.fit(pca_inputs, target)
    coefficients = search.best_estimator_.named_steps["lasso"].coef_
    return {
        "selected": np.flatnonzero(np.abs(coefficients) > 1e-10),
        "alpha": float(search.best_params_["lasso__alpha"]),
    }
```

The first selection finds PCA controls that explain SPY. The second finds PCA
controls related to the candidate. Their union protects the candidate slope from
controls that one predictive equation alone might omit.

```python
def double_selection_test(
    outcome: np.ndarray,
    candidate: np.ndarray,
    controls: np.ndarray,
    pca_inputs: np.ndarray,
    control_names: list[str],
    hac_lags: int,
) -> dict:
    """Run post-double-selection and HAC inference for one candidate slope."""
    outcome_selection = fit_lasso_selector(outcome, pca_inputs)
    candidate_selection = fit_lasso_selector(candidate, pca_inputs)
    selected_union = sorted(
        set(outcome_selection["selected"]) | set(candidate_selection["selected"])
    )

    final_design = candidate.reshape(-1, 1)
    if selected_union:
        final_design = np.column_stack([candidate, controls[:, selected_union]])
    inference = ols_hac_test(outcome, final_design, hac_lags)

    return {
        **inference,
        "n_outcome": len(outcome_selection["selected"]),
        "n_candidate": len(candidate_selection["selected"]),
        "n_union": len(selected_union),
        "outcome_alpha": outcome_selection["alpha"],
        "candidate_alpha": candidate_selection["alpha"],
        "selected_names": [control_names[i] for i in selected_union],
    }
```

## 6. Run post-double-selection for the four candidates

```python
post_results = []
for index, name in enumerate(candidate_names):
    result = double_selection_test(
        outcome_trimmed,
        candidates[:, index],
        controls_trimmed,
        zoo_returns_inference,
        factor_names,
        HAC_LAGS,
    )
    post_results.append({"factor": name, **result})

post_table = pl.DataFrame(post_results).select(
    "factor",
    pl.col("coef").round(4),
    pl.col("se").round(4),
    pl.col("t_stat").round(2),
    pl.col("p_value").round(4),
    "n_outcome",
    "n_candidate",
    "n_union",
)
post_table
```

## 7. Compare uncertainty and selection breadth

The left panel compares the same SPY loading before and after PCA conditioning, with
Newey-West intervals at the conventional two-sided level. The right panel shows what each
of the two LASSO steps selected, since the union that enters the final regression is only
as interesting as the two selections behind it: an outcome equation that keeps the whole
basis makes the union the whole basis whatever the candidate equation chose, and
post-double-selection then reduces to controlling for everything.

```python
factor_positions = np.arange(len(candidate_names))
naive_coef = np.array([result["coef"] for result in naive_results])
naive_se = np.array([result["se"] for result in naive_results])
post_coef = np.array([result["coef"] for result in post_results])
post_se = np.array([result["se"] for result in post_results])
outcome_selected = np.array([result["n_outcome"] for result in post_results])
candidate_selected = np.array([result["n_candidate"] for result in post_results])
figure_subtitle = (
    f"HAC estimates; {N}-ETF factor zoo; {return_wide['timestamp'].min()} "
    f"to {return_wide['timestamp'].max()}"
)
```

The coefficient panel uses position and marker shape as well as color, so the
comparison remains legible in grayscale. The selection bars start at zero.

```python
fig, axes = plt.subplots(1, 2, figsize=FIGSIZE["dual_h_tall"], sharey=True)

axes[0].errorbar(
    naive_coef,
    factor_positions - 0.12,
    xerr=1.96 * naive_se,
    fmt="o",
    color=COLORS["neutral"],
    capsize=3,
    label="Naive",
)
axes[0].errorbar(
    post_coef,
    factor_positions + 0.12,
    xerr=1.96 * post_se,
    fmt="s",
    color=COLORS["blue"],
    capsize=3,
    label="Post-double-selection",
)
zero_line(axes[0], axis="x")
axes[0].set_yticks(factor_positions, candidate_names)
axes[0].set_xlabel("SPY loading (slope)")
_ = axes[0].legend(loc="best")

axes[1].barh(
    factor_positions - 0.18,
    outcome_selected,
    height=0.34,
    color=COLORS["blue"],
    alpha=0.85,
    label="Selected for SPY",
)
axes[1].barh(
    factor_positions + 0.18,
    candidate_selected,
    height=0.34,
    color=COLORS["amber"],
    alpha=0.85,
    label="Selected for the candidate",
)
axes[1].set_xlim(0, N_PCA_FACTORS)
axes[1].set_xlabel("Selected PCA controls (count)")
# The outcome LASSO retains all ten components for every candidate, so every navy bar spans
# the full axis and no corner inside the panel is free. A key placed in one sits on a bar of
# its own colour and cannot be read; it goes above the panel instead.
axes[1].legend(loc="lower left", bbox_to_anchor=(0.0, 1.0), ncol=2, frameon=False, fontsize=8)
axes[1].invert_yaxis()

add_message_title(
    axes[0],
    "SPY loading before and after PCA conditioning",
    subtitle=figure_subtitle,
)
show_with_alt(
    fig,
    "Two panels sharing a vertical axis of candidate factor names. The left panel plots each "
    "factor's SPY loading twice, the naive estimate and the post-double-selection estimate at "
    "slightly offset heights with different marker shapes, each with a horizontal "
    "Newey-West interval and a vertical line at zero. The right panel is a grouped "
    "horizontal bar chart of how many of the ten PCA controls each LASSO selected for that "
    "factor, one bar for the SPY equation and one for the candidate equation, with a legend "
    "naming them; the regression uses their union.",
)
```

**Interpretation**: the naive slopes mix each managed factor's association with
SPY and its correlation with broad co-movement. Once the PCA basis enters, the
coefficient magnitudes contract toward zero and their HAC intervals include zero.
The outcome-selection LASSO retains the full ten-component basis, so this
low-dimensional example effectively becomes a conservative all-PC spanning test.

This is not a direct replication of Feng, Giglio, and Xiu (2020). Their target is
a cross-sectional SDF loading estimated from many test assets. Here the target is
a conditional time-series loading for one disjoint ETF, which isolates the
post-double-selection mechanics without claiming an SDF or causal estimand.

## Key takeaways

1. **Compare like with like**: both columns test the SPY loading on a candidate;
   the conditional version differs only by the selected PCA controls.
2. **Temporal validation matters inside selection**: expanding folds and
   fold-local scaling and PCA keep later observations out of earlier validation.
3. **Serial dependence changes uncertainty**: intercept-inclusive Newey-West
   inference replaces both IID mean tests and a manual HC1 calculation.
4. **Missing data are part of the estimand**: the analysis uses a balanced
   long-history panel instead of converting pre-inception observations to zeros.
5. **Scope remains limited**: the coefficients are in-sample associations,
   conditional on estimated PCs and a current curated ETF universe. Cross-sectional
   SDF inference and point-in-time asset-pricing validation require a richer design.

**Connection to Section 15.3**: post-double-selection protects one target
coefficient by taking the union of controls predictive of the outcome and of the
focal factor. The next step for a production study would add a point-in-time
universe, cross-sectional test assets, and inference designed for estimated SDF
loadings.
![notebook output](figures/p1_1.png)