From dd7b197cd6ab97b9c23b6bab87c5d7b51aabf533 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20St=C3=B6ckl?= Date: Sun, 13 Sep 2026 11:35:26 +0200 Subject: [PATCH] Add an opt-in override for the fewer-observations-than-assets check statistical.factor.model() refuses to fit whenever m < N. That is the one shape a statistical factor model is normally reached for -- a universe wider than its history -- but the restriction is sound guidance and stays the default. This adds a way past it for callers who can say why. statistical.factor.model(R, k = 1, ..., allow.fewer.obs = FALSE) With the default, nothing changes: the same check, the same message, the same results, and tests/testthat/test-stat-factor-model.R is untouched and still passes unchanged on this branch. With allow.fewer.obs = TRUE the fit proceeds and k is bounded by m - 2, which is where the arithmetic binds: extractCovariance() divides the residual sums of squares by m - k - 1, so k = m - 1 would return an infinite covariance rather than an error. The bound lives inside the m < N branch, so the m >= N path is not touched in any respect. The number of assets does not enter it, because on this path there are more assets than observations. The documentation carries what the option costs, as asked: why it is FALSE by default, what a short wide panel does to sample second moments, when lifting the restriction is defensible, and citations on both sides. For: Connor and Korajczyk (1986, 1988), whose asymptotic principal components are built for a cross-section that grows against a fixed time series; Bai and Ng (2002) for choosing k in that regime; Fan, Liao and Mincheva (2013) for the covariance estimator built on it. Against: Michaud (1989) on error maximization, Ledoit and Wolf (2004) and Fan, Fan and Lv (2008) on why regularization exists, and DeMiguel, Garlappi and Uppal (2009) on losing to equal weighting out of sample. The section says plainly that the covariance is full rank by construction, which is not the same as being well estimated. Tests are a new file. Ten tests pin both sides: that the default refuses with the message it always had, with the argument absent and with it FALSE; that the opt-in path fits, and is bounded, and that every admissible k returns a finite covariance; that a positional third argument still reaches prcomp(); and that the tall and square cases are identical with the option set either way. No existing test is modified. Co-Authored-By: Claude Opus 5 --- R/stat.factor.model.R | 94 ++++++++++++++++- man/statistical.factor.model.Rd | 78 +++++++++++++- tests/testthat/test-factor-model-wide.R | 130 ++++++++++++++++++++++++ 3 files changed, 299 insertions(+), 3 deletions(-) create mode 100644 tests/testthat/test-factor-model-wide.R diff --git a/R/stat.factor.model.R b/R/stat.factor.model.R index 504ba1a2..cbe018f1 100644 --- a/R/stat.factor.model.R +++ b/R/stat.factor.model.R @@ -22,9 +22,54 @@ #' loadings, factor realizations, and residuals are computed and returned #' given the number of factors used for the model. #' +#' By default the model is not fitted when there are fewer observations than +#' assets. The decomposition itself does not require it -- \code{prcomp} +#' returns \code{min(m, N)} components and the fit uses \code{k} of them -- +#' but a short, wide panel is where sample second moments are least +#' trustworthy. The largest sample eigenvalues are biased upward and the +#' smallest downward as \code{N/m} grows, so a factor structure extracted +#' there can be mostly noise, and any instability in the underlying series is +#' harder to detect in a short window than in a long one. This is the standard +#' guidance and it is the default. +#' +#' \code{allow.fewer.obs = TRUE} lifts the restriction. The case for lifting it +#' is the asymptotic principal components literature, which was built for this +#' shape: Connor and Korajczyk (1986, 1988) show that the factor realizations +#' are consistently estimated as the cross-section grows with the time series +#' held fixed, which is the opposite of the usual asymptotics and the reason +#' the method exists. Bai and Ng (2002) give criteria for choosing \code{k} in +#' that regime, and Fan, Liao and Mincheva (2013) build a covariance estimator +#' on top of it for exactly the case \code{N > m}. +#' +#' The case against lifting it is not that the fit fails but that it succeeds +#' quietly. The covariance returned is full rank by construction, which is not +#' the same as being well estimated: it will invert, and an optimizer will take +#' the inverse at face value. That is the error-maximizing behaviour Michaud +#' (1989) describes, and the reason shrinkage (Ledoit and Wolf 2004) and +#' factor-based regularization (Fan, Fan and Lv 2008) exist. DeMiguel, +#' Garlappi and Uppal (2009) is the reminder that a portfolio built on a badly +#' estimated covariance can lose to equal weighting out of sample. +#' +#' Set it to \code{TRUE} when the factors are wanted as a low-dimensional +#' summary rather than as an estimate of every pairwise correlation, when +#' \code{k} is small relative to the number of observations, and when the +#' result is judged out of sample rather than in it. Do not set it merely to +#' obtain an invertible covariance matrix from data that does not support one. +#' +#' With the restriction lifted, \code{k} is bounded by \code{m - 2}: +#' \code{extractCovariance} divides the residual sums of squares by +#' \code{m - k - 1}, so \code{k = m - 1} would return an infinite covariance +#' rather than an error. The number of assets does not enter, because on this +#' path there are more of them than observations. The bound applies only here; +#' behaviour when \code{m >= N} is unchanged in every respect. +#' #' @param R xts of asset returns #' @param k number of factors to use #' @param \dots additional arguments passed to \code{prcomp} +#' @param allow.fewer.obs logical. Fit the model even when there are fewer +#' observations than assets. \code{FALSE} by default; see Details for when +#' lifting the restriction is defensible and what it costs. Must be named, as +#' it follows \code{\dots}. #' @return #' #' \describe{ #' \item{factor_loadings}{ N x k matrix of factor loadings (i.e. betas)} @@ -34,8 +79,37 @@ #' } #' Where N is the number of assets, k is the number of factors, and m is the #' number of observations. +#' @references +#' Bai, J. and Ng, S. (2002). Determining the Number of Factors in Approximate +#' Factor Models. \emph{Econometrica}, 70(1), 191-221. +#' +#' Connor, G. and Korajczyk, R. A. (1986). Performance Measurement with the +#' Arbitrage Pricing Theory: A New Framework for Analysis. \emph{Journal of +#' Financial Economics}, 15(3), 373-394. +#' +#' Connor, G. and Korajczyk, R. A. (1988). Risk and Return in an Equilibrium +#' APT: Application of a New Test Methodology. \emph{Journal of Financial +#' Economics}, 21(2), 255-289. +#' +#' DeMiguel, V., Garlappi, L. and Uppal, R. (2009). Optimal Versus Naive +#' Diversification: How Inefficient is the 1/N Portfolio Strategy? +#' \emph{Review of Financial Studies}, 22(5), 1915-1953. +#' +#' Fan, J., Fan, Y. and Lv, J. (2008). High Dimensional Covariance Matrix +#' Estimation Using a Factor Model. \emph{Journal of Econometrics}, 147(1), +#' 186-197. +#' +#' Fan, J., Liao, Y. and Mincheva, M. (2013). Large Covariance Estimation by +#' Thresholding Principal Orthogonal Complements. \emph{Journal of the Royal +#' Statistical Society: Series B}, 75(4), 603-680. +#' +#' Ledoit, O. and Wolf, M. (2004). Honey, I Shrunk the Sample Covariance +#' Matrix. \emph{Journal of Portfolio Management}, 30(4), 110-119. +#' +#' Michaud, R. O. (1989). The Markowitz Optimization Enigma: Is Optimized +#' Optimal? \emph{Financial Analysts Journal}, 45(1), 31-42. #' @export -statistical.factor.model <- function(R, k=1, ...){ +statistical.factor.model <- function(R, k=1, ..., allow.fewer.obs=FALSE){ if(!is.xts(R)){ R <- try(as.xts(R)) if(inherits(R, "try-error")) stop("R must be an xts object or coercible to an xts object") @@ -45,13 +119,29 @@ statistical.factor.model <- function(R, k=1, ...){ N <- ncol(R) # checks for R - if(m < N) stop("fewer observations than assets") + if(m < N && !isTRUE(allow.fewer.obs)) stop("fewer observations than assets") x <- coredata(R) # Make sure k is an integer if(k <= 0) stop("k must be a positive integer") k <- as.integer(k) + if(m < N){ + # Only reachable with allow.fewer.obs = TRUE. prcomp() returns min(m, N) + # components, of which min(m - 1, N) carry variance once the data has been + # centred. The binding limit is one below that, because extractCovariance() + # divides the residual sums of squares by m - k - 1, so k = m - 1 would + # return an infinite covariance rather than an error. Nothing here touches + # the m >= N path. + max_k <- m - 2L + if(max_k < 1L) + stop("at least three observations are needed to fit a factor model") + if(k > max_k) + stop(sprintf(paste("k = %d requests more factors than the data supports:", + "%d observations and %d assets give at most %d"), + k, m, N, max_k)) + } + # Fit a statistical factor model using Principal Component Analysis (PCA) fit <- prcomp(x, ...=...) diff --git a/man/statistical.factor.model.Rd b/man/statistical.factor.model.Rd index 223f0bcb..86497f4f 100644 --- a/man/statistical.factor.model.Rd +++ b/man/statistical.factor.model.Rd @@ -4,7 +4,7 @@ \alias{statistical.factor.model} \title{Statistical Factor Model} \usage{ -statistical.factor.model(R, k = 1, ...) +statistical.factor.model(R, k = 1, ..., allow.fewer.obs = FALSE) } \arguments{ \item{R}{xts of asset returns} @@ -12,6 +12,11 @@ statistical.factor.model(R, k = 1, ...) \item{k}{number of factors to use} \item{\dots}{additional arguments passed to \code{prcomp}} + +\item{allow.fewer.obs}{logical. Fit the model even when there are fewer +observations than assets. \code{FALSE} by default; see Details for when +lifting the restriction is defensible and what it costs. Must be named, as +it follows \code{\dots}.} } \value{ #' \describe{ @@ -30,4 +35,75 @@ Fit a statistical factor model using Principal Component Analysis (PCA) The statistical factor model is fitted using \code{prcomp}. The factor loadings, factor realizations, and residuals are computed and returned given the number of factors used for the model. + +By default the model is not fitted when there are fewer observations than +assets. The decomposition itself does not require it -- \code{prcomp} +returns \code{min(m, N)} components and the fit uses \code{k} of them -- +but a short, wide panel is where sample second moments are least +trustworthy. The largest sample eigenvalues are biased upward and the +smallest downward as \code{N/m} grows, so a factor structure extracted +there can be mostly noise, and any instability in the underlying series is +harder to detect in a short window than in a long one. This is the standard +guidance and it is the default. + +\code{allow.fewer.obs = TRUE} lifts the restriction. The case for lifting it +is the asymptotic principal components literature, which was built for this +shape: Connor and Korajczyk (1986, 1988) show that the factor realizations +are consistently estimated as the cross-section grows with the time series +held fixed, which is the opposite of the usual asymptotics and the reason +the method exists. Bai and Ng (2002) give criteria for choosing \code{k} in +that regime, and Fan, Liao and Mincheva (2013) build a covariance estimator +on top of it for exactly the case \code{N > m}. + +The case against lifting it is not that the fit fails but that it succeeds +quietly. The covariance returned is full rank by construction, which is not +the same as being well estimated: it will invert, and an optimizer will take +the inverse at face value. That is the error-maximizing behaviour Michaud +(1989) describes, and the reason shrinkage (Ledoit and Wolf 2004) and +factor-based regularization (Fan, Fan and Lv 2008) exist. DeMiguel, +Garlappi and Uppal (2009) is the reminder that a portfolio built on a badly +estimated covariance can lose to equal weighting out of sample. + +Set it to \code{TRUE} when the factors are wanted as a low-dimensional +summary rather than as an estimate of every pairwise correlation, when +\code{k} is small relative to the number of observations, and when the +result is judged out of sample rather than in it. Do not set it merely to +obtain an invertible covariance matrix from data that does not support one. + +With the restriction lifted, \code{k} is bounded by \code{m - 2}: +\code{extractCovariance} divides the residual sums of squares by +\code{m - k - 1}, so \code{k = m - 1} would return an infinite covariance +rather than an error. The number of assets does not enter, because on this +path there are more of them than observations. The bound applies only here; +behaviour when \code{m >= N} is unchanged in every respect. +} +\references{ +Bai, J. and Ng, S. (2002). Determining the Number of Factors in Approximate +Factor Models. \emph{Econometrica}, 70(1), 191-221. + +Connor, G. and Korajczyk, R. A. (1986). Performance Measurement with the +Arbitrage Pricing Theory: A New Framework for Analysis. \emph{Journal of +Financial Economics}, 15(3), 373-394. + +Connor, G. and Korajczyk, R. A. (1988). Risk and Return in an Equilibrium +APT: Application of a New Test Methodology. \emph{Journal of Financial +Economics}, 21(2), 255-289. + +DeMiguel, V., Garlappi, L. and Uppal, R. (2009). Optimal Versus Naive +Diversification: How Inefficient is the 1/N Portfolio Strategy? +\emph{Review of Financial Studies}, 22(5), 1915-1953. + +Fan, J., Fan, Y. and Lv, J. (2008). High Dimensional Covariance Matrix +Estimation Using a Factor Model. \emph{Journal of Econometrics}, 147(1), +186-197. + +Fan, J., Liao, Y. and Mincheva, M. (2013). Large Covariance Estimation by +Thresholding Principal Orthogonal Complements. \emph{Journal of the Royal +Statistical Society: Series B}, 75(4), 603-680. + +Ledoit, O. and Wolf, M. (2004). Honey, I Shrunk the Sample Covariance +Matrix. \emph{Journal of Portfolio Management}, 30(4), 110-119. + +Michaud, R. O. (1989). The Markowitz Optimization Enigma: Is Optimized +Optimal? \emph{Financial Analysts Journal}, 45(1), 31-42. } diff --git a/tests/testthat/test-factor-model-wide.R b/tests/testthat/test-factor-model-wide.R new file mode 100644 index 00000000..b8bec5cf --- /dev/null +++ b/tests/testthat/test-factor-model-wide.R @@ -0,0 +1,130 @@ +############################################################################### +# tests/testthat/test-factor-model-wide.R +# +# statistical.factor.model() refuses a panel with fewer observations than +# assets. That restriction is the default and stays the default: a short, wide +# panel is where sample second moments are least trustworthy. +# +# allow.fewer.obs = TRUE lifts it, for the case the asymptotic principal +# components literature is built for. These tests pin both sides: that the +# default is unchanged in every respect, and that the opt-in path is bounded +# where the arithmetic requires it. +############################################################################### + +skip_on_cran() +library(PortfolioAnalytics) + +utils::data(edhec) +R <- edhec[, 1:8] # 8 assets, plenty of months +wide <- R[1:5, ] # 5 observations, 8 assets: m < N + +test_that("the default refuses a wide panel, with the message it always had", { + expect_error(statistical.factor.model(wide, k = 2), + "fewer observations than assets") + expect_error(statistical.factor.model(wide, k = 2, allow.fewer.obs = FALSE), + "fewer observations than assets") +}) + +test_that("allow.fewer.obs = TRUE fits the wide panel", { + expect_no_error(statistical.factor.model(wide, k = 2, allow.fewer.obs = TRUE)) +}) + +test_that("the wide fit produces a covariance of the right shape and rank", { + S <- extractCovariance( + statistical.factor.model(wide, k = 2, allow.fewer.obs = TRUE)) + expect_equal(dim(S), c(ncol(wide), ncol(wide))) + expect_equal(qr(S)$rank, ncol(wide)) + expect_gt(min(eigen(S, symmetric = TRUE, only.values = TRUE)$values), 0) +}) + +test_that("on the opt-in path k is bounded by m - 2", { + m <- nrow(wide) + # extractCovariance() divides by m - k - 1, so k = m - 1 would return an + # infinite covariance instead of an error. The asset count does not enter: + # on this path there are more assets than observations. + expect_no_error(statistical.factor.model(wide, k = m - 2L, + allow.fewer.obs = TRUE)) + expect_error(statistical.factor.model(wide, k = m - 1L, + allow.fewer.obs = TRUE), + "requests more factors than the data supports") + expect_error(statistical.factor.model(wide, k = 99L, allow.fewer.obs = TRUE), + "requests more factors than the data supports") +}) + +test_that("every k the bound admits yields a finite covariance", { + # The bound exists for this reason, so check the reason and not only the + # boundary. + for (k in seq_len(nrow(wide) - 2L)) { + S <- extractCovariance( + statistical.factor.model(wide, k = k, allow.fewer.obs = TRUE)) + expect_true(all(is.finite(S)), info = paste("k =", k)) + expect_gt(min(diag(S)), 0) + } +}) + +test_that("the message says what was asked for and what was available", { + msg <- tryCatch( + statistical.factor.model(wide, k = 99L, allow.fewer.obs = TRUE), + error = function(e) conditionMessage(e)) + expect_match(msg, "k = 99") + expect_match(msg, as.character(nrow(wide))) + expect_match(msg, as.character(ncol(wide))) +}) + +test_that("too few observations for any factor are refused", { + # m - 2 < 1 leaves no admissible k at all. + expect_error(statistical.factor.model(R[1, ], k = 1, allow.fewer.obs = TRUE), + "at least three observations") + expect_error(statistical.factor.model(R[1:2, ], k = 1, + allow.fewer.obs = TRUE), + "at least three observations") + expect_no_error(statistical.factor.model(R[1:3, ], k = 1, + allow.fewer.obs = TRUE)) +}) + +test_that("k <= 0 is still refused, and before the new check", { + expect_error(statistical.factor.model(wide, k = 0, allow.fewer.obs = TRUE), + "positive integer") + expect_error(statistical.factor.model(wide, k = -1, allow.fewer.obs = TRUE), + "positive integer") +}) + +test_that("the argument must be named, so nothing positional changes meaning", { + # allow.fewer.obs follows ..., so a third positional argument still reaches + # prcomp() as it always did. + expect_error(statistical.factor.model(wide, 2, TRUE), + "fewer observations than assets") +}) + +test_that("a tall panel is untouched, with the option set either way", { + tall <- R[1:60, ] + a <- extractCovariance(statistical.factor.model(tall, k = 3)) + b <- extractCovariance(statistical.factor.model(tall, k = 3, + allow.fewer.obs = TRUE)) + expect_equal(as.numeric(a), as.numeric(b)) + + # Reproduce the decomposition independently. + x <- zoo::coredata(tall) + pc <- prcomp(x) + betas <- pc$rotation[, 1:3, drop = FALSE] + f <- x %*% betas + resid <- x - f %*% t(betas) + resid <- sweep(resid, 2, colMeans(resid)) + # extractCovariance() divides the residual sums of squares by m - k - 1, + # not by m - 1, so var() is not the right denominator here. + stock_m2 <- colSums(resid^2) / (nrow(x) - 3L - 1L) + S_ref <- betas %*% stats::cov(f) %*% t(betas) + diag(stock_m2) + + expect_equal(as.numeric(a), as.numeric(S_ref), tolerance = 1e-8) +}) + +test_that("the square case is untouched", { + # m = N is not the wide path, so nothing here changes: a full basis + # reconstructs the data and the covariance comes back finite, as before. + square <- R[seq_len(ncol(R)), ] # m = N = 8 + S <- extractCovariance(statistical.factor.model(square, k = ncol(square))) + expect_true(all(is.finite(S))) + expect_equal(as.numeric(S), + as.numeric(stats::cov(zoo::coredata(square))), + tolerance = 1e-12) +})