Case study: sparse factors in S&P 500 returns

The data

The dataset snp500 is the market-deflated correlation matrix of daily log-returns for p = 423 S&P 500 constituents with complete price histories from January 2010 to December 2019 (n = 2,515 trading days).

library("msPCA")
data(snp500)

dim(snp500)
#> [1] 423 423
round(snp500[1:4, 1:4], 3)
#>           A   AAPL   ABT   ACGL
#> A     0.433  0.004 0.031 -0.027
#> AAPL  0.004  0.716 0.012 -0.048
#> ABT   0.031  0.012 0.614  0.010
#> ACGL -0.027 -0.048 0.010  0.631

The matrix was constructed from the S&P 500 daily update dataset on Kaggle, released under CC0 1.0. Typically, stock returns are dominated by a “market factor” that absorbs a disproportionate share of total variance. To expose cross-sectional structure (sector and style effects) rather than market-wide movements, we projected out the leading eigenvector \(v_1\) of the empirical correlation matrix \(\Sigma\), so \(\texttt{snp500} = P^\top \Sigma P\) with \(P = I - v_1 v_1^\top\). The matrix is rank \(p - 1=422\):

ev <- eigen(snp500, symmetric = TRUE, only.values = TRUE)$values
sum(ev > 1e-8)          # 422: the market direction has been removed
#> [1] 422

The full data processing script is provided in data-raw/snp500.R in the package repository, and ?snp500 documents the format.

Sparse factor extraction

We extract r = 4 sparse factors, each allowed to load on at most k stocks, varying k from 5 to 35 in steps of 5 and running the analysis under both constraint types. The code below produces the results for the orthogonality constraint (feasibilityConstraintType = 0); setting feasibilityConstraintType = 1 gives the zero pairwise correlation results.

ks_grid <- seq(5, 35, by = 5)

results <- lapply(ks_grid, function(k) {
  set.seed(42)
  res <- mspca(snp500, r = 4, ks = rep(k, 4), verbose = FALSE,
               maxIter = 100, feasibilityConstraintType = 0)
  data.frame(
    k      = k,
    fve    = fraction_variance_explained(snp500, res$x_best),
    orth   = feasibility_violation_off(snp500, res$x_best, 0),
    pwcorr = feasibility_violation_off(snp500, res$x_best, 1)
  )
})
results_df <- do.call(rbind, results)

We also run nsprcomp::nsprcomp() (Sigg 2019) at the same budgets as a reference. Note that nsprcomp() requires a data matrix rather than a covariance matrix, so it needs the deflated returns XR <- X %*% P rather than snp500. That matrix is 2,515 x 423 and is not shipped with the package; data-raw/snp500.R documents how to rebuild it from the raw prices. We include a pre-computed comparison instead here, so the numbers below can be inspected and re-plotted without a rerun:

res_grid <- read.csv(system.file("vignette-data", "snp_varyingk_results.csv",
                                 package = "msPCA"))

ks  <- sort(unique(res_grid$k))
by_constraint <- function(cn) {
  sub <- res_grid[res_grid$constraint == cn, ]
  sub[match(ks, sub$k), c("fve", "orth_violation")]
}

tab <- cbind(k = ks,
             by_constraint("orthogonality"),
             by_constraint("zero-correlation"),
             by_constraint("nsprcomp"))

knitr::kable(
  tab, digits = 4, row.names = FALSE,
  col.names = c("k", "FVE (msPCA - orth)", "orth viol (msPCA - orth)", "FVE (msPCA - zero-corr)",
                "orth viol (msPCA - zero-corr)", "FVE (nsprcomp)", "orth viol (nsprcomp)"),
  caption = paste("FVE and orthogonality violation across the sparsity grid.",
                  "Every violation column reports the orthogonality violation,",
                  "including for the fits run under the zero-correlation constraint.")
)

Comparing the two constraint types

The figures below report the fraction of variance explained (FVE), the orthogonality violation, and the uncorrelatedness violation as a function of k, for msPCA under each constraint type, against nsprcomp::nsprcomp() at the same sparsity budgets.

Left: fraction of variance explained vs. sparsity budget k. Center: orthogonality violation vs. k. Right: uncorrelatedness violation vs. k. Results are shown for msPCA with orthogonality constraints (blue, solid), msPCA with zero pairwise correlation constraints (green, dashed), and nsprcomp::nsprcomp() (orange, dotted). All methods use r = 4 components.

On this dataset, nsprcomp::nsprcomp() returns exactly orthogonal loading vectors for small-to-moderate sparsity budgets (k <= 20), reflecting the effectiveness of the deflation procedure when component supports can easily be disjoint. Beyond that orthogonality breaks down, and does so steeply: the violation is 0.015 at k = 25, 0.12 at k = 30 and 0.23 at k = 35.

By contrast, msPCA with orthogonality constraints holds the violation at or below 1e-4 — the default feasibility tolerance — at every budget from k = 10 upward, because the penalty on constraint violation is explicitly tightened throughout the algorithm. In terms of FVE the two methods are comparable, with a small edge for nsprcomp::nsprcomp() up to k = 25. At k = 30 and k = 35 — precisely the budgets where nsprcomp::nsprcomp() gives up orthogonality — msPCA overtakes it, reaching an FVE of 0.111 and 0.114 against 0.102 and 0.105 while keeping the violation below tolerance.

Neither nsprcomp::nsprcomp() nor orthogonality-constrained msPCA yields uncorrelated PCs here. To obtain uncorrelated PCs we run msPCA with pairwise correlation constraints instead, which yields PCs with near-zero pairwise correlation that are not mutually orthogonal. On this dataset, requiring zero pairwise correlation rather than orthogonality is possible only at the expense of a substantially lower FVE.

The two constraints correspond to different feasibility definitions and lead to meaningfully different factor compositions. msPCA lets the user choose and enforce whichever is relevant to their use case, with predictable behavior across the full range of sparsity levels.

Interpreting the sparse components

Fixing k = 10, each factor loads on 10 stocks out of 423, making it possible to associate each component with an economic theme. Sector labels below follow the Global Industry Classification Standard [GICS; MSCI and S&P Dow Jones Indices (2023)].

set.seed(42)
res_orth <- mspca(snp500, r = 4, ks = rep(10, 4), verbose = FALSE,
                  maxIter = 100, feasibilityConstraintType = 0)
set.seed(42)
res_corr <- mspca(snp500, r = 4, ks = rep(10, 4), verbose = FALSE,
                  maxIter = 100, feasibilityConstraintType = 1)

print(res_orth)
print(res_corr)

summary() on either fit reports the violations under the constraint that fit enforced, and labels them as such:

summary(res_orth)
summary(res_corr)

Loadings of the 4 PCs (sparsity k = 10) returned by msPCA with orthogonality (left) or zero-correlation (right) constraints.

Takeaway

Orthogonality and zero pairwise correlation are not interchangeable. On strongly correlated data they recover qualitatively different factor structures. On this stock return data, orthogonality returns disjoint supports and spends all four components inside two sectors: three separate utility clusters and one REIT basket. Zero pairwise correlation, on the other hand, returns overlapping supports and recovers four distinct themes instead —utilities, REITs, packaged food and household staples, and casino operators— each expressed as a contrast between positively and negatively loaded groups. We also observe a difference in the distribution of variance explained by each PC: under zero correlation the leading component carries most of the explained variance (2.28% against 0.45–0.64% for the others), whereas orthogonality spreads it more evenly (2.13% down to 1.21%), and total FVE is markedly lower with zero-correlation constraints at every sparsity budget.