Replicating He et al. (2023)
In sdim: An R Package for Supervised Dimension Reduction

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

This vignette reproduces Table 3 of He, Huang, Li, and Zhou (2023). The table compares four ways of summarising a large set of candidate factor proxies down to a few risk factors that price the 48 Fama-French value-weighted industry portfolios. Performance is measured by the total adjusted $R^2$ (%) of the pricing regressions.

The four methods are:

FF: use the Fama-French five factors plus momentum directly, taking the first $K$ of them.
PCA: extract the $K$ principal components that maximise the variance of the factor-proxy matrix $G$.
PLS: extract the $K$ components of $G$ that are most predictive of the test-asset returns.
RRA: the reduced-rank approach of He et al. (2023). It searches for $K$ linear combinations of the proxies in $G$ whose loadings best explain the returns of the basis assets. Formally, the loading matrix in the regression of returns on proxies is restricted to rank $K$, and that restriction is what shrinks $L \approx 70$ proxies down to a handful of usable factors.

The motivation is empirical: many candidate factors have been proposed in the literature, but most carry little incremental pricing information once one accounts for the others. RRA is designed to find the small linear subspace that retains the pricing-relevant content.

Setup

The bundled he2023_* datasets come from the authors' replication package. The factor proxies in he2023_factors end twelve months earlier than the portfolio panels, so we slice the rows to align them and convert percentages to decimals. Returns are taken in excess of the one-month Treasury bill rate RF:

library(sdim)

he2023_ff48 <- he2023_ff48vw[1:516, -1] / 100 - he2023_ff5$RF[127:642] / 100
G <- he2023_factors[1:516, -1] / 100

# First 6 columns of G are Fama-French 5 + momentum
f5 <- G[, 1:6]

Replication

We loop over the same factor counts as the paper. For each $K$:

the FF row uses the first $K$ of the six Fama-French/momentum factors directly (only defined for $K \le 6$);
the PCA, PLS, and RRA rows fit the corresponding *_est() function on the full proxy set $G$, then pass the extracted factors to eval_factors() to get the total adjusted $R^2$ on the test assets:

nfact   <- c(1, 3, 5, 6, 10)
methods <- c("FF", "PCA", "PLS", "RRA")

total_r2 <- matrix(NA, nrow = length(methods), ncol = length(nfact))
rownames(total_r2) <- methods
colnames(total_r2) <- paste(nfact, "factors")

for (j in seq_along(nfact)) {

  k <- nfact[j]

  if (k <= 6) {

    total_r2["FF", j] <- eval_factors(he2023_ff48, f5[, 1:k])["TotalR2"]

  }

  fit_pca <- pca_est(target = he2023_ff48, X = G, nfac = k)
  total_r2["PCA", j] <- eval_factors(he2023_ff48, fit_pca$factors)["TotalR2"]

  fit_pls <- pls_est(target = he2023_ff48, X = G, nfac = k)
  total_r2["PLS", j] <- eval_factors(he2023_ff48, fit_pls$factors)["TotalR2"]

  fit_rra <- rra_est(target = he2023_ff48, X = G, nfac = k)
  total_r2["RRA", j] <- eval_factors(he2023_ff48, fit_rra$factors)["TotalR2"]

}

Results

round(total_r2, 2)

RRA delivers the highest total adjusted $R^2$ at every factor count. This is the headline finding of He et al. (2023): once we look for factors that are constructed to price the basis assets — rather than factors that maximise own-variance (PCA) or predictive covariance with returns one column at a time (PLS) — a small number of linear combinations of the 70 proxies recovers nearly all the pricing information.