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