add giglio xiu 2021

NAMESPACE

@@ -4,6 +4,7 @@ export(ChenFang2019BetaRankTest)
 export(FGXFactorsTest)
 export(FRP)
 export(GKRFactorScreening)
+export(GiglioXiu2021RiskPremia)
 export(HACcovariance)
 export(HJMisspecificationDistance)
 export(IterativeKleibergenPaap2006BetaRankTest)

R/RcppExports.R

@@ -13,6 +13,10 @@ FRPCpp <- function(returns, factors, misspecification_robust = TRUE, include_sta
     .Call(`_intrinsicFRP_FRPCpp`, returns, factors, misspecification_robust, include_standard_errors, hac_prewhite, target_level_gkr2014_screening)
 }
 
+GiglioXiu2021RiskPremiaCpp <- function(returns, factors, which_n_pca = 0L) {
+    .Call(`_intrinsicFRP_GiglioXiu2021RiskPremiaCpp`, returns, factors, which_n_pca)
+}
+
 GKRFactorScreeningCpp <- function(returns, factors, target_level = 0.05, hac_prewhite = FALSE) {
     .Call(`_intrinsicFRP_GKRFactorScreeningCpp`, returns, factors, target_level, hac_prewhite)
 }

R/frp.R

@@ -109,3 +109,65 @@ FRP = function(
   ))
 
 }
+
+#' @title Compute Factor Risk Premia using Giglio and Xiu (2021) Method
+#'
+#' @description Computes the factor risk premia using the 3-step procedure of
+#' Giglio and Xiu (2021) <doi:10.1086/714090>. The procedure consists in
+#' 1) extracting p PCA from returns, where p is a consistent estimator of the
+#' number of strong latent factors, 2) compute a cross-sectional regression of
+#' average returns on the estimated betas of the p latent factors, and 3) compute
+#' a time-series regression of observed factors on the p latent factors to obtain
+#' their mimicking portfolio risk premia.
+#'
+#' @param returns A `n_observations x n_returns` matrix of asset returns.
+#' @param factors A `n_observations x n_factors` matrix of risk factors.
+#' @param which_n_pca Method to determine the number of PCA components:
+#' `0` for Giglio Xiu (2021) <doi:10.1086/714090>,
+#' `-1` for Ahn Horenstein (2013) <doi:10.3982/ECTA8968>,
+#' or any positive integer for a specific number of PCs.
+#' @param check_arguments A boolean: `TRUE` for internal check of all function
+#' arguments; `FALSE` otherwise. Default is `TRUE`.
+#'
+#' @return A list containing:
+#' \describe{
+#'   \item{risk_premia}{A matrix of factor risk premia.}
+#'   \item{n_pca}{The number of principal components used.}
+#' }
+#' @examples
+#' \dontrun{
+#' returns <- matrix(rnorm(200), nrow=20, ncol=10)
+#' factors <- matrix(rnorm(60), nrow=20, ncol=3)
+#' which_n_pca <- 0
+#' result <- GiglioXiu2021RiskPremia(returns, factors, which_n_pca)
+#' print(result)
+#' }
+#' @export
+GiglioXiu2021RiskPremia = function(
+  returns,
+  factors,
+  which_n_pca = 0,
+  check_arguments = TRUE
+) {
+
+  # check function arguments
+  if (check_arguments) {
+
+    stopifnot("`returns` must contain numeric values" = is.numeric(returns))
+    stopifnot("`factors` must contain numeric values" = is.numeric(factors))
+    stopifnot("`returns` and `factors` must have the same number of rows" = nrow(returns) == nrow(factors))
+    stopifnot("`returns` must not contain missing values (NA/NaN)" = !anyNA(returns))
+    stopifnot("`factors` must not contain missing values (NA/NaN)" = !anyNA(factors))
+
+    stopifnot("`which_n_pca` must be numeric" = is.numeric(which_n_pca))
+
+  }
+
+  # Ensure the function is defined in the C++ file and compiled
+  return(.Call(`_intrinsicFRP_GiglioXiu2021RiskPremiaCpp`,
+    returns,
+    factors,
+    which_n_pca
+  ))
+
+}

src/RcppExports.cpp

@@ -54,6 +54,19 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// GiglioXiu2021RiskPremiaCpp
+Rcpp::List GiglioXiu2021RiskPremiaCpp(const arma::mat& returns, const arma::mat& factors, const int which_n_pca);
+RcppExport SEXP _intrinsicFRP_GiglioXiu2021RiskPremiaCpp(SEXP returnsSEXP, SEXP factorsSEXP, SEXP which_n_pcaSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const arma::mat& >::type returns(returnsSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type factors(factorsSEXP);
+    Rcpp::traits::input_parameter< const int >::type which_n_pca(which_n_pcaSEXP);
+    rcpp_result_gen = Rcpp::wrap(GiglioXiu2021RiskPremiaCpp(returns, factors, which_n_pca));
+    return rcpp_result_gen;
+END_RCPP
+}
 // GKRFactorScreeningCpp
 Rcpp::List GKRFactorScreeningCpp(const arma::mat& returns, const arma::mat& factors, const double target_level, const bool hac_prewhite);
 RcppExport SEXP _intrinsicFRP_GKRFactorScreeningCpp(SEXP returnsSEXP, SEXP factorsSEXP, SEXP target_levelSEXP, SEXP hac_prewhiteSEXP) {
@@ -194,6 +207,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_intrinsicFRP_FGXThreePassCovarianceCpp", (DL_FUNC) &_intrinsicFRP_FGXThreePassCovarianceCpp, 4},
     {"_intrinsicFRP_FGXThreePassCovarianceNoControlsCpp", (DL_FUNC) &_intrinsicFRP_FGXThreePassCovarianceNoControlsCpp, 3},
     {"_intrinsicFRP_FRPCpp", (DL_FUNC) &_intrinsicFRP_FRPCpp, 6},
+    {"_intrinsicFRP_GiglioXiu2021RiskPremiaCpp", (DL_FUNC) &_intrinsicFRP_GiglioXiu2021RiskPremiaCpp, 3},
     {"_intrinsicFRP_GKRFactorScreeningCpp", (DL_FUNC) &_intrinsicFRP_GKRFactorScreeningCpp, 4},
     {"_intrinsicFRP_HACCovarianceMatrixCpp", (DL_FUNC) &_intrinsicFRP_HACCovarianceMatrixCpp, 2},
     {"_intrinsicFRP_HACVarianceCpp", (DL_FUNC) &_intrinsicFRP_HACVarianceCpp, 2},

src/frp.cpp

@@ -4,6 +4,7 @@
 #include "hac_covariance.h"
 #include "gkr_factor_screening.h"
 #include "utils.h"
+#include "n_pca.h"
 
 /////////////////
 ///// FRPCpp ////
@@ -289,3 +290,78 @@ arma::vec StandardErrorsKRSFRPCpp(
     std::sqrt(static_cast<double>(returns.n_rows));
 
 }
+
+//////////////////////////////////////
+///// GiglioXiu2021RiskPremiaCpp ////
+
+Rcpp::List GiglioXiu2021RiskPremiaCpp(
+  const arma::mat& returns,
+  const arma::mat& factors,
+  const int which_n_pca
+) {
+
+  // parameters
+  const unsigned int n_assets = returns.n_cols;
+  const unsigned int n_observations = returns.n_rows;
+
+  // center returns
+  const arma::rowvec mean_returns = arma::mean(returns, 0);
+  const arma::mat centered_returns = returns.each_row() - mean_returns;
+  const arma::mat centered_factors = factors.each_row() - arma::mean(factors, 0);
+
+  // returns decomposition
+  arma::vec e_vals;
+  arma::mat U;
+  arma::mat V;
+  arma::svd(U, e_vals, V, centered_returns.t() / (n_assets * n_observations));
+  e_vals = arma::square(e_vals);
+
+
+  // compute the number of PCA
+  unsigned int n_pca;
+
+  if (which_n_pca == 0) {
+
+    // use the method in Giglio Xiu 2021
+    unsigned int n_max_pca = std::floor(2 * n_assets / 3);
+    n_pca = NPCA_GiglioXiu2021Cpp(e_vals, n_assets, n_observations, n_max_pca);
+
+  } else if (which_n_pca < 0) {
+
+    // use the method in Ahn Horenstein 2013
+    unsigned int n_max_pca = std::floor(2 * n_assets / 3);
+    n_pca = NPCA_AhnHorenstein2013Cpp(e_vals, n_max_pca)["er"];
+
+  } else {
+
+    // set the number of PCA to the supplied value
+    n_pca = which_n_pca;
+
+  }
+
+  // latent factors and corresponding beta
+  const arma::mat pca_factors = std::sqrt(n_observations) * V.head_cols(n_pca);
+  const arma::mat beta_hat = centered_returns.t() * pca_factors / n_observations;
+
+  // cross-sectional regression of average returns on estimated betas
+  const arma::vec gamma = SolveSympd(
+    beta_hat.t() * beta_hat,
+    beta_hat.t() * mean_returns.t()
+  );
+
+  // time-series regression of factors on latent factors
+  const arma::mat eta = SolveSympd(
+    pca_factors.t() * pca_factors,
+    pca_factors.t() * centered_factors
+  ).t();
+
+  // factor risk premia
+  const arma::vec risk_premia = eta * gamma;
+
+  // output list containing factor risk premia and the number of pca used
+  return Rcpp::List::create(
+    Rcpp::Named("risk_premia") = risk_premia,
+    Rcpp::Named("n_pca") = n_pca
+  );
+
+}

src/frp.h

@@ -133,4 +133,12 @@ arma::vec StandardErrorsKRSFRPCpp(
   const bool hac_prewhite = false
 );
 
+// Compute the factor risk premia estimates of Giglio Xiu 2021.
+// [[Rcpp::export]]
+Rcpp::List GiglioXiu2021RiskPremiaCpp(
+  const arma::mat& returns,
+  const arma::mat& factors,
+  const int which_n_pca = 0
+);
+
 #endif

src/n_pca.cpp

@@ -0,0 +1,78 @@
+// Author: Alberto Quaini
+
+#include "utils.h"
+
+////////////////////////////////////////////
+/////// NPCA_GiglioXiu2021Cpp //////////////
+
+unsigned int NPCA_GiglioXiu2021Cpp(
+  const arma::vec& e_vals,
+  const unsigned int n_assets,
+  const unsigned int n_observations,
+  unsigned int p_max
+) {
+
+  // minimum between n_assets and n_observations
+  const unsigned int min_NT = std::min(n_assets, n_observations);
+
+  // adjust p_max if it is too large
+  if (p_max > min_NT) {
+    p_max = min_NT - 1;
+  }
+
+  // compute the penalty term
+  const double scaling = 0.5 * arma::median(e_vals);
+  const double penalty = scaling * (std::log(n_assets) + std::log(n_observations)) *
+    (std::pow(n_assets, -0.5) + std::pow(n_observations, -0.5));
+
+  // compute the criterion for each index
+  arma::vec criterion(p_max);
+  for (unsigned int i = 0; i < p_max; ++i) {
+    criterion(i) = e_vals(i) + (i + 1) * penalty;
+  }
+
+  // return the index that minimizes the criterion
+  return criterion.index_min();
+
+}
+
+////////////////////////////////////////////////
+/////// NPCA_AhnHorenstein2013Cpp //////////////
+
+Rcpp::List NPCA_AhnHorenstein2013Cpp(
+  const arma::vec& evals,
+  unsigned int n_max = 0
+) {
+
+  // set the number of eigenvalues
+  const unsigned int n_evals = evals.n_elem;
+
+  // if n_max <= zero or >= n_evals, set it to n_evals - 1
+  if (n_max <= 0 || n_max >= n_evals) {
+    n_max = n_evals - 1;
+  }
+
+  const arma::uvec top = arma::regspace<arma::uvec>(0, n_max - 1);
+  const arma::uvec bottom = arma::regspace<arma::uvec>(1, n_max);
+
+  // calculate the eigenvalue ratios and the index where the ratio is maximized
+  const arma::vec evals_ratio = evals(top) / evals(bottom);
+  const unsigned int er = evals_ratio.index_max() + 1; // +1 to match R's 1-based index
+
+  // calculate the sum of the eigenvalues after each index
+  arma::vec evals_revcumsum = arma::reverse(arma::cumsum(arma::reverse(evals)));
+  evals_revcumsum.shed_row(n_evals - 1); // remove the last element
+  evals_revcumsum.insert_rows(n_evals - 1, arma::zeros(1)); // add zero at the end
+
+  // divide the eigenvalues by the sum
+  const arma::vec evals_tr = evals / evals_revcumsum;
+  const arma::vec growth_ratio = arma::log(
+    1 + evals_tr(top)) / arma::log(1 + evals_tr(bottom)
+  );
+  const unsigned int gr = growth_ratio.index_max() + 1; // +1 to match R's 1-based index
+
+  return Rcpp::List::create(
+    Rcpp::Named("er") = er,
+    Rcpp::Named("gr") = gr
+  );
+}

src/n_pca.h

@@ -0,0 +1,23 @@
+// Author: Alberto Quaini
+
+#ifndef N_PCA_H
+#define N_PCA_H
+
+#include <RcppArmadillo.h>
+
+// Function for internal use.
+// Implement the procedure to determine the number of PCA
+// that summarize the risk in returns of Giglio Xiu (2021).
+unsigned int NPCA_GiglioXiu2021Cpp(
+  const arma::vec& e_vals,
+  const unsigned int n_assets,
+  const unsigned int n_observations,
+  unsigned int p_max
+);
+
+Rcpp::List NPCA_AhnHorenstein2013Cpp(
+  const arma::vec& evals,
+  unsigned int n_max = 0
+);
+
+#endif