own definition of solve and inv

README.md

@@ -138,7 +138,7 @@ covariance matrix of the model's pricing errors $\varepsilon_t$.
 
 To cite `intrinsicFRP` in publications, please use:
 
-> Quaini, A. (2023). `intrinsicFRP`: An R Package for Factor Model Asset Pricing. `R` package version 1.0.0. URL: https://CRAN.R-project.org/package=intrinsicFRP.
+> Quaini, A. (2023). `intrinsicFRP`: An R Package for Factor Model Asset Pricing. `R` package version 2.0.0. URL: https://CRAN.R-project.org/package=intrinsicFRP.
 
 ## Installation
 

src/adaptive_weights.cpp

@@ -1,6 +1,7 @@
 // Author: Alberto Quaini
 
 #include "adaptive_weights.h"
+#include "utils.h"
 
 ///////////////////////////////
 ///// AdaptiveWeightsCpp //////
@@ -21,18 +22,16 @@ arma::vec AdaptiveWeightsCpp(
   case 'a' :
     // Type 'a': Compute weights based on the first-step intrinsic factor risk premia estimates
     return AdaptiveWeightsFromVectorCpp(
-      arma::cov(factors, returns) * arma::solve(
+      arma::cov(factors, returns) * SolveSympd(
         arma::cov(returns),
-        arma::mean(returns).t(),
-        arma::solve_opts::likely_sympd
+        arma::mean(returns).t()
     ));
 
   case 'b' :
     // Type 'b': Compute weights based on the matrix of factors regression coefficients on test asset excess returns
-    return AdaptiveWeightsFromMatrixCpp(arma::solve(
+    return AdaptiveWeightsFromMatrixCpp(SolveSympd(
       arma::cov(factors),
-      arma::cov(factors, returns),
-      arma::solve_opts::likely_sympd
+      arma::cov(factors, returns)
     ));
 
   default :

src/frp.cpp

@@ -3,6 +3,7 @@
 #include "frp.h"
 #include "hac_covariance.h"
 #include "gkr_factor_screening.h"
+#include "utils.h"
 
 /////////////////
 ///// FRPCpp ////
@@ -88,10 +89,9 @@ Rcpp::List ReturnFRPCpp(
     const arma::vec mean_returns = arma::mean(returns).t();
 
     // Solve for beta coefficients
-    const arma::mat beta = arma::solve(
+    const arma::mat beta = SolveSympd(
       arma::cov(factors),
-      covariance_factors_returns,
-      arma::solve_opts::likely_sympd
+      covariance_factors_returns
     ).t();
 
     // Compute risk premia based on whether misspecification is robust or not
@@ -154,10 +154,9 @@ Rcpp::List ReturnFRPCpp(
   } else {
 
     // Solve for beta coefficients when standard errors are not included
-    const arma::mat beta = arma::solve(
+    const arma::mat beta = SolveSympd(
       arma::cov(factors),
-      arma::cov(factors, returns),
-      arma::solve_opts::likely_sympd
+      arma::cov(factors, returns)
     ).t();
 
     // Return risk premia only in a list
@@ -187,10 +186,7 @@ arma::vec FMFRPCpp(
 ) {
 
   // Solve for risk premia using Fama-MacBeth method
-  return arma::solve(
-    beta.t() * beta,
-    beta.t(),
-    arma::solve_opts::likely_sympd) * mean_returns;
+  return SolveSympd(beta.t() * beta, beta.t()) * mean_returns;
 
 }
 
@@ -204,17 +200,15 @@ arma::vec KRSFRPCpp(
 ) {
 
   // Solve for beta coefficients with the weighting matrix
-  const arma::mat beta_t_wei_mat_inv = arma::solve(
+  const arma::mat beta_t_wei_mat_inv = SolveSympd(
     weighting_matrix,
-    beta,
-    arma::solve_opts::likely_sympd
+    beta
   ).t();
 
   // Compute risk premia using Kan-Robotti-Shanken method
-  return arma::solve(
+  return SolveSympd(
     beta_t_wei_mat_inv * beta,
-    beta_t_wei_mat_inv,
-    arma::solve_opts::likely_sympd
+    beta_t_wei_mat_inv
   ) * mean_returns;
 
 }
@@ -234,7 +228,7 @@ arma::vec StandardErrorsFRPCpp(
 ) {
 
   // Compute the matrix H as the inverse of beta' * beta
-  const arma::mat h_matrix  = arma::inv_sympd(beta.t() * beta);
+  const arma::mat h_matrix  = InvSympd(beta.t() * beta);
 
   // Compute matrix A used in the standard error calculation
   const arma::mat a_matrix = h_matrix * beta.t();
@@ -252,10 +246,9 @@ arma::vec StandardErrorsFRPCpp(
 
   // Compute variance of factors and solve for fac_centred_var_fac_inv
   const arma::mat variance_factors = arma::cov(factors);
-  const arma::mat fac_centred_var_fac_inv = arma::solve(
+  const arma::mat fac_centred_var_fac_inv = SolveSympd(
     variance_factors,
-    factors_centred.t(),
-    arma::solve_opts::likely_sympd
+    factors_centred.t()
   ).t();
 
   // Compute matrix Z used in standard error calculation
@@ -267,10 +260,9 @@ arma::vec StandardErrorsFRPCpp(
     // mean term
     (gamma.each_row() - (a_matrix * mean_returns).t()) -
     // beta term
-    phi_centred.each_col() % (factors_centred * arma::solve(
+    phi_centred.each_col() % (factors_centred * SolveSympd(
         variance_factors,
-        gamma_true,
-        arma::solve_opts::likely_sympd
+        gamma_true
     )) +
     // error term
     (fac_centred_var_fac_inv.each_col() % (
@@ -298,15 +290,14 @@ arma::vec StandardErrorsKRSFRPCpp(
 ) {
 
   // Compute the inverse of variance_returns and useful matrices that use it
-  const arma::mat var_ret_inv = arma::inv_sympd(variance_returns);
+  const arma::mat var_ret_inv = InvSympd(variance_returns);
   const arma::mat var_ret_inv_beta = var_ret_inv * beta;
   const arma::vec var_ret_inv_mean_ret = var_ret_inv * mean_returns;
 
   // Compute matrix A used in the standard error calculation
-  const arma::mat a_matrix = arma::solve(
+  const arma::mat a_matrix = SolveSympd(
     beta.t() * var_ret_inv_beta,
-    var_ret_inv_beta.t(),
-    arma::solve_opts::likely_sympd
+    var_ret_inv_beta.t()
   );
 
   // Center returns and factors
@@ -316,11 +307,10 @@ arma::vec StandardErrorsKRSFRPCpp(
   // Compute intermediate terms for standard error calculation
   const arma::vec var_ret_inv_err_krs = var_ret_inv_mean_ret -
     var_ret_inv_beta * krs_frp;
-  const arma::mat var_fac_inv = arma::inv_sympd(arma::cov(factors));
-  const arma::mat hkrs_var_fac_inv = arma::solve(
+  const arma::mat var_fac_inv = InvSympd(arma::cov(factors));
+  const arma::mat hkrs_var_fac_inv = SolveSympd(
     beta.t() * var_ret_inv_beta,
-    var_fac_inv,
-    arma::solve_opts::likely_sympd
+    var_fac_inv
   );
 
   // Compute terms used in the standard error calculation

src/gkr_factor_screening.cpp

@@ -3,6 +3,7 @@
 #include "gkr_factor_screening.h"
 #include "hac_covariance.h"
 #include "frp.h"
+#include "utils.h"
 
 ////////////////////////////////
 ///// GKRFactorScreeningCpp ////
@@ -145,17 +146,15 @@ arma::vec GKRSDFCoefficientsCpp(
 ) {
 
   // Solve for `V[R]^{-1} * Cov[R, F]
-  const arma::mat var_ret_inv_cov_ret_fac = arma::solve(
+  const arma::mat var_ret_inv_cov_ret_fac = SolveSympd(
     variance_returns,
-    covariance_returns_factors,
-    arma::solve_opts::likely_sympd
+    covariance_returns_factors
   );
 
   // Compute sdf coefficients using the GKR method
-  return arma::solve(
+  return SolveSympd(
     var_ret_inv_cov_ret_fac.t() * covariance_returns_factors,
-    var_ret_inv_cov_ret_fac.t(),
-    arma::solve_opts::likely_sympd
+    var_ret_inv_cov_ret_fac.t()
   ) * mean_returns;
 
 }
@@ -174,12 +173,12 @@ arma::vec StandardErrorsGKRSDFCoefficientsCpp(
 ) {
 
   // Compute the inverse of variance_returns and useful matrices that use it
-  const arma::mat var_ret_inv = arma::inv_sympd(variance_returns);
+  const arma::mat var_ret_inv = InvSympd(variance_returns);
   const arma::mat var_ret_inv_cov_ret_fac = var_ret_inv * covariance_returns_factors;
   const arma::vec var_ret_inv_mean_ret = var_ret_inv * mean_returns;
 
   // Compute matrix H used in the standard error calculation
-  const arma::mat h_matrix = arma::inv_sympd(
+  const arma::mat h_matrix = InvSympd(
     covariance_returns_factors.t() * var_ret_inv_cov_ret_fac
   );
 

src/hac_covariance.cpp

@@ -244,11 +244,20 @@ void HACRevertPrewhiteningCpp(
   const unsigned int n_variables = coefficients.n_rows;
 
   // Compute the inverse of the identity -  coefficients.t()
-  const arma::mat inv_temp = arma::inv(
-    arma::eye(n_variables, n_variables) - coefficients.t()
-  );
+  arma::mat inv_temp;
+  try {
+    // Try arma::inv_sympd
+    inv_temp = arma::inv(arma::eye(n_variables, n_variables) - coefficients.t(), arma::inv_opts::no_ugly);
+
+  } catch (const std::runtime_error&) {
+    // Fallback to generalized inverse
+    inv_temp = pinv(
+      arma::eye(n_variables, n_variables) - coefficients.t()
+    );
+
+  }
 
-  // Adjust the HAC covariance matrix using the inverse of the prewhitening transformation
+  // Adjust the HAC covariance matrix using the inverse of the pre-whitening transformation
   hac_covariance = inv_temp * hac_covariance * inv_temp.t();
 
 }

src/hj_misspecification_distance.cpp

@@ -2,6 +2,7 @@
 
 #include "hj_misspecification_distance.h"
 #include "hac_covariance.h"
+#include "utils.h"
 
 ///////////////////////////////////////
 //// HJMisspecificationDistanceCpp ////
@@ -38,26 +39,23 @@ Rcpp::List HJMisspecificationDistanceCpp(
 ) {
 
   // Compute the inverse variance of returns.
-  const arma::mat inv_var_ret = arma::inv_sympd(variance_returns);
+  const arma::mat inv_var_ret = InvSympd(variance_returns);
 
   // Calculate required vectors and matrices to compure the KRS SDF coefficients.
-  const arma::vec var_ret_inv_mean_ret = arma::solve(
+  const arma::vec var_ret_inv_mean_ret = SolveSympd(
     variance_returns,
-    mean_returns,
-    arma::solve_opts::likely_sympd
+    mean_returns
   );
   const arma::mat covariance_factors_returns = arma::cov(factors, returns);
-  const arma::mat var_ret_inv_cov_ret_fac = arma::solve(
+  const arma::mat var_ret_inv_cov_ret_fac = SolveSympd(
     variance_returns,
-    covariance_factors_returns.t(),
-    arma::solve_opts::likely_sympd
+    covariance_factors_returns.t()
   );
 
   // Compute the KRS SDF coefficients.
-  const arma::vec krs_sdf_coefficients = arma::solve(
+  const arma::vec krs_sdf_coefficients = SolveSympd(
     covariance_factors_returns * var_ret_inv_cov_ret_fac,
-    covariance_factors_returns,
-    arma::solve_opts::likely_sympd
+    covariance_factors_returns
   ) * var_ret_inv_mean_ret;
 
   // Calculate the squared Hansen-Jagannathan (HJ) distance.

src/identification_tests.cpp

@@ -1,6 +1,7 @@
 // Author: Alberto Quaini
 
 #include "identification_tests.h"
+#include "utils.h"
 
 ///////////////////////////////////////
 ///// ChenFang2019BetaRankTestCpp /////
@@ -144,15 +145,13 @@ arma::vec2 KleibergenPaap2006BetaRankTestStatisticAndPvalueCpp(
   const arma::mat sqrt_V22 = arma::sqrtmat_sympd(V22 * V22.t());
 
   // Compute the A_qperp and B_qperp matrices according to equation (12) in the KP 2006 paper
-  const arma::mat A_qperp = U.tail_cols(n_factors - q) * arma::solve(
+  const arma::mat A_qperp = U.tail_cols(n_factors - q) * SolveSympd(
     U22,
-    sqrt_U22,
-    arma::solve_opts::likely_sympd
+    sqrt_U22
   );
-  const arma::mat B_qperp = sqrt_V22 * arma::solve(
+  const arma::mat B_qperp = sqrt_V22 * SolveSympd(
     V22.t(),
-    V.tail_cols(n_returns - q).t(),
-    arma::solve_opts::likely_sympd
+    V.tail_cols(n_returns - q).t()
   );
 
   // Compute the Kronecker product of B_qperp and A_qperp transposed
@@ -167,10 +166,9 @@ arma::vec2 KleibergenPaap2006BetaRankTestStatisticAndPvalueCpp(
 
   // Calculate the test statistic, which does not require multiplication by n_observations
   // This is based on the quadratic form of lambda_q and Omega_q
-  const double statistic = arma::dot(lambda_q, arma::solve(
+  const double statistic = arma::dot(lambda_q, SolveSympd(
     Omega_q,
-    lambda_q,
-    arma::solve_opts::likely_sympd
+    lambda_q
   ));
 
   // Calculate the p-value using the Chi-squared distribution with degrees of freedom

src/tfrp.cpp

@@ -2,6 +2,7 @@
 
 #include "tfrp.h"
 #include "hac_covariance.h"
+#include "utils.h"
 
 //////////////////
 ///// TFRPCpp ////
@@ -61,10 +62,9 @@ arma::vec TFRPCpp(
 ) {
 
   // Return the risk premia computed using the formula Cov(F, R) * Var(R)^(-1) * E(R)
-  return covariance_factors_returns * arma::solve(
+  return covariance_factors_returns * SolveSympd(
     variance_returns,
-    mean_returns,
-    arma::solve_opts::likely_sympd
+    mean_returns
   );
 
 }
@@ -82,7 +82,7 @@ arma::vec StandardErrorsTFRPCpp(
 ) {
 
   // Compute the inverse of variance returns applied to two different matrices
-  const arma::mat var_ret_inv = arma::inv_sympd(variance_returns);
+  const arma::mat var_ret_inv = InvSympd(variance_returns);
   const arma::mat var_ret_inv_cov_ret_fac = var_ret_inv *
     covariance_factors_returns.t();