Add an opt-in override for the fewer-observations-than-assets check

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, ...=...)
 

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)
+})