finalised verify_volatility with methods and tests #26

NAMESPACE

@@ -13,6 +13,8 @@ S3method(forecast,PosteriorBSVARMIX)
 S3method(forecast,PosteriorBSVARMSH)
 S3method(forecast,PosteriorBSVARSV)
 S3method(verify_volatility,PosteriorBSVAR)
+S3method(verify_volatility,PosteriorBSVARMIX)
+S3method(verify_volatility,PosteriorBSVARMSH)
 S3method(verify_volatility,PosteriorBSVARSV)
 export(compute_conditional_sd)
 export(compute_fitted_values)

NEWS.md

@@ -1,5 +1,7 @@
 # bsvars 2.0.0.9000
 
+1. Included Bayesian procedures for verifying structural shocks' heteroskedastiicty equation-by-equation using Savage-Dickey density rations [#26](https://github.com/bsvars/bsvars/issues/26)
+
 The package is under intensive development, and more functionalities will be provided soon! To see the package [ROADMAP](https://github.com/bsvars/bsvars/milestone/3) towards the next version 2.1.0.
 
 Have a question, or suggestion, or wanna get in touch? Join the package [DISCUSSION](https://github.com/bsvars/bsvars/discussions) forum.

R/verify.R

@@ -49,7 +49,7 @@
 #' specification  = specify_bsvar_sv$new(us_fiscal_lsuw, p = 1)
 #' set.seed(123)
 #' 
-#' # run the burn-in
+#' # estimate the model
 #' posterior      = estimate(specification, 60, thin = 1)
 #' 
 #' # verify heteroskedasticity
@@ -74,6 +74,48 @@ verify_volatility <- function(posterior) {
 
 
 
+#' @inherit verify_volatility
+#' @method verify_volatility PosteriorBSVAR
+#' @inheritParams verify_volatility
+#'
+#' @description Displays information that the model is homoskedastic.
+#' 
+#' @return Nothing. Just displays a message: The model is homoskedastic.
+#'
+#' @examples
+#' # simple workflow
+#' ############################################################
+#' # upload data
+#' data(us_fiscal_lsuw)
+#' 
+#' # specify the model and set seed
+#' specification  = specify_bsvar$new(us_fiscal_lsuw, p = 1)
+#' set.seed(123)
+#' 
+#' # estimate the model
+#' posterior      = estimate(specification, 5, thin = 1)
+#' 
+#' # verify heteroskedasticity
+#' sddr           = verify_volatility(posterior)
+#' 
+#' # workflow with the pipe |>
+#' ############################################################
+#' set.seed(123)
+#' us_fiscal_lsuw |>
+#'   specify_bsvar$new(p = 1) |>
+#'   estimate(S = 5, thin = 1) |> 
+#'   verify_volatility() -> sddr
+#'   
+#' @export
+verify_volatility.PosteriorBSVAR <- function(posterior) {
+  message("The model is homoskedastic.")
+}
+
+
+
+
+
+
 #' @inherit verify_volatility
 #' @method verify_volatility PosteriorBSVARSV
 #' @inheritParams verify_volatility
@@ -102,7 +144,7 @@ verify_volatility <- function(posterior) {
 #' specification  = specify_bsvar_sv$new(us_fiscal_lsuw, p = 1)
 #' set.seed(123)
 #' 
-#' # run the burn-in
+#' # estimate the model
 #' posterior      = estimate(specification, 60, thin = 1)
 #' 
 #' # verify heteroskedasticity
@@ -130,7 +172,7 @@ verify_volatility.PosteriorBSVARSV <- function(posterior) {
     X               = posterior$last_draw$data_matrices$X
 
     # estimate the SDDR
-    sddr            = .Call(`_bsvars_verify_volatility_cpp`, just_posterior, prior, Y, X, TRUE)
+    sddr            = .Call(`_bsvars_verify_volatility_sv_cpp`, just_posterior, prior, Y, X, TRUE)
 
     class(sddr)     = "SDDR"
     return(sddr)
@@ -139,13 +181,24 @@ verify_volatility.PosteriorBSVARSV <- function(posterior) {
 
 
 
+
 #' @inherit verify_volatility
-#' @method verify_volatility PosteriorBSVAR
+#' @method verify_volatility PosteriorBSVARMIX
 #' @inheritParams verify_volatility
 #'
-#' @description Displays information that the model is homoskedastic.
+#' @description Computes the logarithm of Bayes factor for the homoskedasticity hypothesis 
+#' for each of the structural shocks via Savage-Dickey Density Ration (SDDR).
+#' The hypothesis of homoskedasticity is represented by restriction:
+#' \deqn{H_0: \sigma^2_{n.1} = ... = \sigma^2_{n.M} = 1}
+#' The logarithm of Bayes factor for this hypothesis can be computed using the SDDR 
+#' as the difference of logarithms of the marginal posterior distribution ordinate at the restriction 
+#' less the marginal prior distribution ordinate at the same point:
+#' \deqn{log p(\omega_n = 0 | data) - log p(\omega_n = 0)}
+#' Therefore, a negative value of the difference is the evidence against 
+#' homoskedasticity of the structural shock. The estimation of both elements of the difference requires 
+#' numerical integration.
 #' 
-#' @return Nothing. Just displays a message: The model is homoskedastic.
+#' @seealso \code{\link{specify_bsvar_mix}}, \code{\link{estimate}}
 #'
 #' @examples
 #' # simple workflow
@@ -154,11 +207,11 @@ verify_volatility.PosteriorBSVARSV <- function(posterior) {
 #' data(us_fiscal_lsuw)
 #' 
 #' # specify the model and set seed
-#' specification  = specify_bsvar$new(us_fiscal_lsuw, p = 1)
+#' specification  = specify_bsvar_mix$new(us_fiscal_lsuw, p = 1, M = 2)
 #' set.seed(123)
 #' 
-#' # run the burn-in
-#' posterior      = estimate(specification, 5, thin = 1)
+#' # estimate the model
+#' posterior      = estimate(specification, 60, thin = 1)
 #' 
 #' # verify heteroskedasticity
 #' sddr           = verify_volatility(posterior)
@@ -167,11 +220,84 @@ verify_volatility.PosteriorBSVARSV <- function(posterior) {
 #' ############################################################
 #' set.seed(123)
 #' us_fiscal_lsuw |>
-#'   specify_bsvar$new(p = 1) |>
-#'   estimate(S = 5, thin = 1) |> 
+#'   specify_bsvar_mix$new(p = 1, M = 2) |>
+#'   estimate(S = 60, thin = 1) |> 
 #'   verify_volatility() -> sddr
 #'   
 #' @export
-verify_volatility.PosteriorBSVAR <- function(posterior) {
-  message("The model is homoskedastic.")
-}
+verify_volatility.PosteriorBSVARMIX <- function(posterior) {
+
+  # get the inputs to estimation
+  just_posterior  = posterior$posterior
+  prior           = posterior$last_draw$prior$get_prior()
+  Y               = posterior$last_draw$data_matrices$Y
+  X               = posterior$last_draw$data_matrices$X
+
+  # estimate the SDDR
+  sddr            = .Call(`_bsvars_verify_volatility_msh_cpp`, just_posterior, prior, Y, X)
+
+  class(sddr)     = "SDDR"
+  return(sddr)
+}
+
+
+
+
+
+#' @inherit verify_volatility
+#' @method verify_volatility PosteriorBSVARMSH
+#' @inheritParams verify_volatility
+#'
+#' @description Computes the logarithm of Bayes factor for the homoskedasticity hypothesis 
+#' for each of the structural shocks via Savage-Dickey Density Ration (SDDR).
+#' The hypothesis of homoskedasticity is represented by restriction:
+#' \deqn{H_0: \sigma^2_{n.1} = ... = \sigma^2_{n.M} = 1}
+#' The logarithm of Bayes factor for this hypothesis can be computed using the SDDR 
+#' as the difference of logarithms of the marginal posterior distribution ordinate at the restriction 
+#' less the marginal prior distribution ordinate at the same point:
+#' \deqn{log p(\omega_n = 0 | data) - log p(\omega_n = 0)}
+#' Therefore, a negative value of the difference is the evidence against 
+#' homoskedasticity of the structural shock. The estimation of both elements of the difference requires 
+#' numerical integration.
+#' 
+#' @seealso \code{\link{specify_bsvar_msh}}, \code{\link{estimate}}
+#'
+#' @examples
+#' # simple workflow
+#' ############################################################
+#' # upload data
+#' data(us_fiscal_lsuw)
+#' 
+#' # specify the model and set seed
+#' specification  = specify_bsvar_msh$new(us_fiscal_lsuw, p = 1, M = 2)
+#' set.seed(123)
+#' 
+#' # estimate the model
+#' posterior      = estimate(specification, 60, thin = 1)
+#' 
+#' # verify heteroskedasticity
+#' sddr           = verify_volatility(posterior)
+#' 
+#' # workflow with the pipe |>
+#' ############################################################
+#' set.seed(123)
+#' us_fiscal_lsuw |>
+#'   specify_bsvar_msh$new(p = 1, M = 2) |>
+#'   estimate(S = 60, thin = 1) |> 
+#'   verify_volatility() -> sddr
+#'   
+#' @export
+verify_volatility.PosteriorBSVARMSH <- function(posterior) {
+
+  # get the inputs to estimation
+  just_posterior  = posterior$posterior
+  prior           = posterior$last_draw$prior$get_prior()
+  Y               = posterior$last_draw$data_matrices$Y
+  X               = posterior$last_draw$data_matrices$X
+
+  # estimate the SDDR
+  sddr            = .Call(`_bsvars_verify_volatility_msh_cpp`, just_posterior, prior, Y, X)
+
+  class(sddr)     = "SDDR"
+  return(sddr)
+}

inst/include/bsvars_RcppExports.h

@@ -898,17 +898,17 @@ namespace bsvars {
         return Rcpp::as<std::string >(rcpp_result_gen);
     }
 
-    inline Rcpp::List verify_volatility_cpp(const Rcpp::List& posterior, const Rcpp::List& prior, const arma::mat& Y, const arma::mat& X, const bool sample_s_ = true) {
-        typedef SEXP(*Ptr_verify_volatility_cpp)(SEXP,SEXP,SEXP,SEXP,SEXP);
-        static Ptr_verify_volatility_cpp p_verify_volatility_cpp = NULL;
-        if (p_verify_volatility_cpp == NULL) {
-            validateSignature("Rcpp::List(*verify_volatility_cpp)(const Rcpp::List&,const Rcpp::List&,const arma::mat&,const arma::mat&,const bool)");
-            p_verify_volatility_cpp = (Ptr_verify_volatility_cpp)R_GetCCallable("bsvars", "_bsvars_verify_volatility_cpp");
+    inline Rcpp::List verify_volatility_sv_cpp(const Rcpp::List& posterior, const Rcpp::List& prior, const arma::mat& Y, const arma::mat& X, const bool sample_s_ = true) {
+        typedef SEXP(*Ptr_verify_volatility_sv_cpp)(SEXP,SEXP,SEXP,SEXP,SEXP);
+        static Ptr_verify_volatility_sv_cpp p_verify_volatility_sv_cpp = NULL;
+        if (p_verify_volatility_sv_cpp == NULL) {
+            validateSignature("Rcpp::List(*verify_volatility_sv_cpp)(const Rcpp::List&,const Rcpp::List&,const arma::mat&,const arma::mat&,const bool)");
+            p_verify_volatility_sv_cpp = (Ptr_verify_volatility_sv_cpp)R_GetCCallable("bsvars", "_bsvars_verify_volatility_sv_cpp");
         }
         RObject rcpp_result_gen;
         {
             RNGScope RCPP_rngScope_gen;
-            rcpp_result_gen = p_verify_volatility_cpp(Shield<SEXP>(Rcpp::wrap(posterior)), Shield<SEXP>(Rcpp::wrap(prior)), Shield<SEXP>(Rcpp::wrap(Y)), Shield<SEXP>(Rcpp::wrap(X)), Shield<SEXP>(Rcpp::wrap(sample_s_)));
+            rcpp_result_gen = p_verify_volatility_sv_cpp(Shield<SEXP>(Rcpp::wrap(posterior)), Shield<SEXP>(Rcpp::wrap(prior)), Shield<SEXP>(Rcpp::wrap(Y)), Shield<SEXP>(Rcpp::wrap(X)), Shield<SEXP>(Rcpp::wrap(sample_s_)));
         }
         if (rcpp_result_gen.inherits("interrupted-error"))
             throw Rcpp::internal::InterruptedException();
@@ -919,11 +919,11 @@ namespace bsvars {
         return Rcpp::as<Rcpp::List >(rcpp_result_gen);
     }
 
-    inline double dig2dirichlet(const arma::vec& x, const arma::vec& a, const arma::vec& b, const bool logarithm = true) {
+    inline double dig2dirichlet(const arma::rowvec& x, const arma::rowvec& a, const arma::rowvec& b, const bool logarithm = true) {
         typedef SEXP(*Ptr_dig2dirichlet)(SEXP,SEXP,SEXP,SEXP);
         static Ptr_dig2dirichlet p_dig2dirichlet = NULL;
         if (p_dig2dirichlet == NULL) {
-            validateSignature("double(*dig2dirichlet)(const arma::vec&,const arma::vec&,const arma::vec&,const bool)");
+            validateSignature("double(*dig2dirichlet)(const arma::rowvec&,const arma::rowvec&,const arma::rowvec&,const bool)");
             p_dig2dirichlet = (Ptr_dig2dirichlet)R_GetCCallable("bsvars", "_bsvars_dig2dirichlet");
         }
         RObject rcpp_result_gen;
@@ -940,11 +940,11 @@ namespace bsvars {
         return Rcpp::as<double >(rcpp_result_gen);
     }
 
-    inline double ddirichlet(const arma::vec& x, const arma::vec& a, const bool logarithm = true) {
+    inline double ddirichlet(const arma::rowvec& x, const arma::rowvec& a, const bool logarithm = true) {
         typedef SEXP(*Ptr_ddirichlet)(SEXP,SEXP,SEXP);
         static Ptr_ddirichlet p_ddirichlet = NULL;
         if (p_ddirichlet == NULL) {
-            validateSignature("double(*ddirichlet)(const arma::vec&,const arma::vec&,const bool)");
+            validateSignature("double(*ddirichlet)(const arma::rowvec&,const arma::rowvec&,const bool)");
             p_ddirichlet = (Ptr_ddirichlet)R_GetCCallable("bsvars", "_bsvars_ddirichlet");
         }
         RObject rcpp_result_gen;
@@ -961,6 +961,27 @@ namespace bsvars {
         return Rcpp::as<double >(rcpp_result_gen);
     }
 
+    inline Rcpp::List verify_volatility_msh_cpp(const Rcpp::List& posterior, const Rcpp::List& prior, const arma::mat& Y, const arma::mat& X) {
+        typedef SEXP(*Ptr_verify_volatility_msh_cpp)(SEXP,SEXP,SEXP,SEXP);
+        static Ptr_verify_volatility_msh_cpp p_verify_volatility_msh_cpp = NULL;
+        if (p_verify_volatility_msh_cpp == NULL) {
+            validateSignature("Rcpp::List(*verify_volatility_msh_cpp)(const Rcpp::List&,const Rcpp::List&,const arma::mat&,const arma::mat&)");
+            p_verify_volatility_msh_cpp = (Ptr_verify_volatility_msh_cpp)R_GetCCallable("bsvars", "_bsvars_verify_volatility_msh_cpp");
+        }
+        RObject rcpp_result_gen;
+        {
+            RNGScope RCPP_rngScope_gen;
+            rcpp_result_gen = p_verify_volatility_msh_cpp(Shield<SEXP>(Rcpp::wrap(posterior)), Shield<SEXP>(Rcpp::wrap(prior)), Shield<SEXP>(Rcpp::wrap(Y)), Shield<SEXP>(Rcpp::wrap(X)));
+        }
+        if (rcpp_result_gen.inherits("interrupted-error"))
+            throw Rcpp::internal::InterruptedException();
+        if (Rcpp::internal::isLongjumpSentinel(rcpp_result_gen))
+            throw Rcpp::LongjumpException(rcpp_result_gen);
+        if (rcpp_result_gen.inherits("try-error"))
+            throw Rcpp::exception(Rcpp::as<std::string>(rcpp_result_gen).c_str());
+        return Rcpp::as<Rcpp::List >(rcpp_result_gen);
+    }
+
 }
 
 #endif // RCPP_bsvars_RCPPEXPORTS_H_GEN_

inst/tinytest/test_verify_volatility.R

@@ -0,0 +1,88 @@
+
+data(us_fiscal_lsuw)
+
+set.seed(1)
+suppressMessages(
+  specification_no1 <- specify_bsvar$new(us_fiscal_lsuw)
+)
+run_no1             <- estimate(specification_no1, 60, 1, show_progress = FALSE)
+
+expect_message(
+  verify_volatility(run_no1),
+  pattern = "*homoskedastic.",
+  info = "verify_volatility: the model to be homoskedastic"
+)
+
+#############################################
+set.seed(1)
+suppressMessages(
+  specification_no1 <- specify_bsvar_sv$new(us_fiscal_lsuw, centred_sv = TRUE)
+)
+run_no1             <- estimate(specification_no1, 60, 1, show_progress = FALSE)
+
+expect_message(
+  verify_volatility(run_no1),
+  pattern = "*known*",
+  info = "verify_volatility: for centred SV"
+)
+
+#############################################
+set.seed(1)
+suppressMessages(
+  specification_no1 <- specify_bsvar_sv$new(us_fiscal_lsuw)
+)
+run_no1             <- estimate(specification_no1, 60, 1, show_progress = FALSE)
+
+expect_true(
+  is.numeric(verify_volatility(run_no1)$logSDDR),
+  info = "verify_volatility: for non-centred SV; logSDDR is numeric"
+)
+
+#############################################
+set.seed(1)
+suppressMessages(
+  specification_no1 <- specify_bsvar_msh$new(us_fiscal_lsuw, M = 2)
+)
+run_no1             <- estimate(specification_no1, 60, 1, show_progress = FALSE)
+
+expect_true(
+  is.numeric(verify_volatility(run_no1)$logSDDR),
+  info = "verify_volatility: for MSH model; logSDDR is numeric"
+)
+
+#############################################
+set.seed(1)
+suppressMessages(
+  specification_no1 <- specify_bsvar_mix$new(us_fiscal_lsuw, M = 2)
+)
+run_no1             <- estimate(specification_no1, 60, 1, show_progress = FALSE)
+
+expect_true(
+  is.numeric(verify_volatility(run_no1)$logSDDR),
+  info = "verify_volatility: for MIX model; logSDDR is numeric"
+)
+
+#############################################
+set.seed(1)
+suppressMessages(
+  specification_no1 <- specify_bsvar_msh$new(us_fiscal_lsuw, finiteM = FALSE)
+)
+run_no1             <- estimate(specification_no1, 60, 1, show_progress = FALSE)
+
+expect_true(
+  is.numeric(verify_volatility(run_no1)$logSDDR),
+  info = "verify_volatility: for sparse MSH model; logSDDR is numeric"
+)
+
+#############################################
+set.seed(1)
+suppressMessages(
+  specification_no1 <- specify_bsvar_mix$new(us_fiscal_lsuw, finiteM = FALSE)
+)
+run_no1             <- estimate(specification_no1, 60, 1, show_progress = FALSE)
+
+expect_true(
+  is.numeric(verify_volatility(run_no1)$logSDDR),
+  info = "verify_volatility: for sparse  MIX model; logSDDR is numeric"
+)
+