Merge branch 'develop'

R/RcppExports.R

@@ -1066,6 +1066,65 @@ sim_wishart <- function(mat_scale, shape) {
     .Call(`_bvhar_sim_wishart`, mat_scale, shape)
 }
 
+#' Quasi-density of GIG
+#' 
+#' @param x postivie support
+#' @param lambda Index of modified Bessel function of third kind.
+#' @param beta Square of the multiplication of the other two parameters.
+#' 
+#' @noRd
+dgig_quasi <- function(x, lambda, beta) {
+    .Call(`_bvhar_dgig_quasi`, x, lambda, beta)
+}
+
+#' AR-Mehod for non-concave part
+#' 
+#' @param num_sim Number to generate process
+#' @param lambda Index of modified Bessel function of third kind.
+#' @param beta Square of the multiplication of the other two parameters.
+#' 
+#' @noRd
+rgig_nonconcave <- function(num_sim, lambda, beta) {
+    .Call(`_bvhar_rgig_nonconcave`, num_sim, lambda, beta)
+}
+
+#' Ratio-of-Uniforms without Mode Shift
+#' 
+#' @param num_sim Number to generate process
+#' @param lambda Index of modified Bessel function of third kind.
+#' @param beta Square of the multiplication of the other two parameters.
+#' 
+#' @noRd
+rgig_without_mode <- function(num_sim, lambda, beta) {
+    .Call(`_bvhar_rgig_without_mode`, num_sim, lambda, beta)
+}
+
+#' Ratio-of-Uniforms with Mode Shift
+#' 
+#' @param num_sim Number to generate process
+#' @param lambda Index of modified Bessel function of third kind.
+#' @param beta Square of the multiplication of the other two parameters.
+#' 
+#' @noRd
+rgig_with_mode <- function(num_sim, lambda, beta) {
+    .Call(`_bvhar_rgig_with_mode`, num_sim, lambda, beta)
+}
+
+#' Generate Generalized Inverse Gaussian Distribution
+#' 
+#' This function samples GIG(lambda, psi, chi) random variates.
+#' 
+#' @param num_sim Number to generate process
+#' @param lambda Index of modified Bessel function of third kind.
+#' @param psi Second parameter of GIG
+#' @param chi Third parameter of GIG
+#' 
+#' @references Hörmann, W., Leydold, J. Generating generalized inverse Gaussian random variates. Stat Comput 24, 547–557 (2014).
+#' @noRd
+sim_gig <- function(num_sim, lambda, psi, chi) {
+    .Call(`_bvhar_sim_gig`, num_sim, lambda, psi, chi)
+}
+
 #' VAR(1) Representation Given VAR Coefficient Matrix
 #' 
 #' Compute the VAR(1) coefficient matrix form

R/bvar-sv.R

@@ -280,27 +280,14 @@ bvar_sv <- function(y,
   # Preprocess the results--------------------------------
   thin_id <- seq(from = 1, to = num_iter - num_burn, by = thinning)
   res$alpha_record <- res$alpha_record[thin_id,]
+  res$h_record <- res$h_record[thin_id,]
   res$a_record <- res$a_record[thin_id,]
   res$h0_record <- res$h0_record[thin_id,]
   res$sigh_record <- res$sigh_record[thin_id,]
-  res$h_record <- split.data.frame(
-    res$h_record,
-    gl(num_design, num_iter + 1)
-  ) %>%
-    lapply(
-      function(x) {
-        colnames(x) <- paste0("h[", seq_len(ncol(x)), "]")
-        x[thin_id + 1,] # since num_iter + 1 rows
-      }
-    )
-  res$h_record <- lapply(
-    seq_along(num_design),
-    function(i) {
-      colnames(res$h_record[[i]]) <- paste0(colnames(res$h_record[[i]]), i)
-      res$h_record[[i]]
-    }
+  colnames(res$h_record) <- paste0(
+    paste0("h[", seq_len(dim_data), "]"),
+    gl(num_design, dim_data)
   )
-  res$h_record <- do.call(cbind, res$h_record)
   res$h_record <- as_draws_df(res$h_record)
   res$coefficients <- matrix(colMeans(res$alpha_record), ncol = dim_data)
   mat_lower <- matrix(0L, nrow = dim_data, ncol = dim_data)

R/bvhar-sv.R

@@ -379,27 +379,14 @@ bvhar_sv <- function(y,
   names(res) <- gsub(pattern = "^alpha", replacement = "phi", x = names(res)) # alpha to phi
   thin_id <- seq(from = 1, to = num_iter - num_burn, by = thinning)
   res$phi_record <- res$phi_record[thin_id,]
+  res$h_record <- res$h_record[thin_id,]
   res$a_record <- res$a_record[thin_id,]
   res$h0_record <- res$h0_record[thin_id,]
   res$sigh_record <- res$sigh_record[thin_id,]
-  res$h_record <- split.data.frame(
-    res$h_record,
-    gl(num_design, num_iter + 1)
-  ) %>%
-    lapply(
-      function(x) {
-        colnames(x) <- paste0("h[", seq_len(ncol(x)), "]")
-        x[thin_id + 1,] # since num_iter + 1 rows
-      }
-    )
-  res$h_record <- lapply(
-    seq_along(num_design),
-    function(x) {
-      colnames(res$h_record[[x]]) <- paste0(colnames(res$h_record[[x]]), x)
-      res$h_record[[x]]
-    }
+  colnames(res$h_record) <- paste0(
+    paste0("h[", seq_len(dim_data), "]"),
+    gl(num_design, dim_data)
   )
-  res$h_record <- do.call(cbind, res$h_record)
   res$h_record <- as_draws_df(res$h_record)
   res$coefficients <- matrix(colMeans(res$phi_record), ncol = dim_data)
   mat_lower <- matrix(0L, nrow = dim_data, ncol = dim_data)

src/RcppExports.cpp

@@ -932,6 +932,72 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// dgig_quasi
+double dgig_quasi(double x, double lambda, double beta);
+RcppExport SEXP _bvhar_dgig_quasi(SEXP xSEXP, SEXP lambdaSEXP, SEXP betaSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< double >::type x(xSEXP);
+    Rcpp::traits::input_parameter< double >::type lambda(lambdaSEXP);
+    Rcpp::traits::input_parameter< double >::type beta(betaSEXP);
+    rcpp_result_gen = Rcpp::wrap(dgig_quasi(x, lambda, beta));
+    return rcpp_result_gen;
+END_RCPP
+}
+// rgig_nonconcave
+Eigen::VectorXd rgig_nonconcave(int num_sim, double lambda, double beta);
+RcppExport SEXP _bvhar_rgig_nonconcave(SEXP num_simSEXP, SEXP lambdaSEXP, SEXP betaSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< int >::type num_sim(num_simSEXP);
+    Rcpp::traits::input_parameter< double >::type lambda(lambdaSEXP);
+    Rcpp::traits::input_parameter< double >::type beta(betaSEXP);
+    rcpp_result_gen = Rcpp::wrap(rgig_nonconcave(num_sim, lambda, beta));
+    return rcpp_result_gen;
+END_RCPP
+}
+// rgig_without_mode
+Eigen::VectorXd rgig_without_mode(int num_sim, double lambda, double beta);
+RcppExport SEXP _bvhar_rgig_without_mode(SEXP num_simSEXP, SEXP lambdaSEXP, SEXP betaSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< int >::type num_sim(num_simSEXP);
+    Rcpp::traits::input_parameter< double >::type lambda(lambdaSEXP);
+    Rcpp::traits::input_parameter< double >::type beta(betaSEXP);
+    rcpp_result_gen = Rcpp::wrap(rgig_without_mode(num_sim, lambda, beta));
+    return rcpp_result_gen;
+END_RCPP
+}
+// rgig_with_mode
+Eigen::VectorXd rgig_with_mode(int num_sim, double lambda, double beta);
+RcppExport SEXP _bvhar_rgig_with_mode(SEXP num_simSEXP, SEXP lambdaSEXP, SEXP betaSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< int >::type num_sim(num_simSEXP);
+    Rcpp::traits::input_parameter< double >::type lambda(lambdaSEXP);
+    Rcpp::traits::input_parameter< double >::type beta(betaSEXP);
+    rcpp_result_gen = Rcpp::wrap(rgig_with_mode(num_sim, lambda, beta));
+    return rcpp_result_gen;
+END_RCPP
+}
+// sim_gig
+Eigen::VectorXd sim_gig(int num_sim, double lambda, double psi, double chi);
+RcppExport SEXP _bvhar_sim_gig(SEXP num_simSEXP, SEXP lambdaSEXP, SEXP psiSEXP, SEXP chiSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< int >::type num_sim(num_simSEXP);
+    Rcpp::traits::input_parameter< double >::type lambda(lambdaSEXP);
+    Rcpp::traits::input_parameter< double >::type psi(psiSEXP);
+    Rcpp::traits::input_parameter< double >::type chi(chiSEXP);
+    rcpp_result_gen = Rcpp::wrap(sim_gig(num_sim, lambda, psi, chi));
+    return rcpp_result_gen;
+END_RCPP
+}
 // compute_stablemat
 Eigen::MatrixXd compute_stablemat(Eigen::MatrixXd x);
 RcppExport SEXP _bvhar_compute_stablemat(SEXP xSEXP) {
@@ -1653,6 +1719,11 @@ static const R_CallMethodDef CallEntries[] = {
     {"_bvhar_sim_iw", (DL_FUNC) &_bvhar_sim_iw, 2},
     {"_bvhar_sim_mniw", (DL_FUNC) &_bvhar_sim_mniw, 5},
     {"_bvhar_sim_wishart", (DL_FUNC) &_bvhar_sim_wishart, 2},
+    {"_bvhar_dgig_quasi", (DL_FUNC) &_bvhar_dgig_quasi, 3},
+    {"_bvhar_rgig_nonconcave", (DL_FUNC) &_bvhar_rgig_nonconcave, 3},
+    {"_bvhar_rgig_without_mode", (DL_FUNC) &_bvhar_rgig_without_mode, 3},
+    {"_bvhar_rgig_with_mode", (DL_FUNC) &_bvhar_rgig_with_mode, 3},
+    {"_bvhar_sim_gig", (DL_FUNC) &_bvhar_sim_gig, 4},
     {"_bvhar_compute_stablemat", (DL_FUNC) &_bvhar_compute_stablemat, 1},
     {"_bvhar_compute_var_stablemat", (DL_FUNC) &_bvhar_compute_var_stablemat, 1},
     {"_bvhar_compute_vhar_stablemat", (DL_FUNC) &_bvhar_compute_vhar_stablemat, 1},

src/bvhardraw.h

@@ -2,8 +2,8 @@
 #define BVHARDRAW_H
 
 #include "bvharomp.h"
+#include "randsim.h"
 #include "bvharmisc.h"
-#include "bvharprob.h"
 
 Eigen::VectorXd build_ssvs_sd(Eigen::VectorXd spike_sd, Eigen::VectorXd slab_sd, Eigen::VectorXd mixture_dummy);
 

src/bvharmisc.h

@@ -3,24 +3,10 @@
 
 Eigen::MatrixXd scale_har(int dim, int week, int month, bool include_mean);
 
-Eigen::MatrixXd sim_mgaussian(int num_sim, Eigen::VectorXd mu, Eigen::MatrixXd sig);
-
-Eigen::MatrixXd sim_mgaussian_chol(int num_sim, Eigen::VectorXd mu, Eigen::MatrixXd sig);
-
-Eigen::MatrixXd sim_mstudent(int num_sim, double df, Eigen::VectorXd mu, Eigen::MatrixXd sig, int method);
-
 Eigen::MatrixXd VARcoeftoVMA(Eigen::MatrixXd var_coef, int var_lag, int lag_max);
 
 Eigen::MatrixXd VHARcoeftoVMA(Eigen::MatrixXd vhar_coef, Eigen::MatrixXd HARtrans_mat, int lag_max);
 
-Eigen::MatrixXd sim_matgaussian(Eigen::MatrixXd mat_mean, Eigen::MatrixXd mat_scale_u, Eigen::MatrixXd mat_scale_v);
-
-Eigen::MatrixXd sim_iw(Eigen::MatrixXd mat_scale, double shape);
-
-Rcpp::List sim_mniw(int num_sim, Eigen::MatrixXd mat_mean, Eigen::MatrixXd mat_scale_u, Eigen::MatrixXd mat_scale, double shape);
-
-Eigen::MatrixXd sim_wishart(Eigen::MatrixXd mat_scale, double shape);
-
 Eigen::MatrixXd kronecker_eigen(Eigen::MatrixXd x, Eigen::MatrixXd y);
 
 Eigen::VectorXd vectorize_eigen(Eigen::MatrixXd x);

src/estimate-hierarchical.cpp

@@ -1,6 +1,6 @@
 #include <RcppEigen.h>
 #include "bvharmisc.h"
-#include "bvharprob.h"
+#include "randsim.h"
 #include <progress.hpp>
 #include <progress_bar.hpp>
 

src/estimate-sv.cpp

@@ -107,7 +107,8 @@ Rcpp::List estimate_var_sv(int num_iter, int num_burn,
   Eigen::MatrixXd contem_coef_record = Eigen::MatrixXd::Zero(num_iter + 1, num_lowerchol); // a = a21, a31, a32, ..., ak1, ..., ak(k-1)
   Eigen::MatrixXd lvol_sig_record = Eigen::MatrixXd::Zero(num_iter + 1, dim); // sigma_h^2 = (sigma_(h1i)^2, ..., sigma_(hki)^2)
   Eigen::MatrixXd lvol_init_record = Eigen::MatrixXd::Zero(num_iter + 1, dim); // h0 = h10, ..., hk0
-  Eigen::MatrixXd lvol_record = Eigen::MatrixXd::Zero(num_design * (num_iter + 1), dim); // time-varying h = (h_1, ..., h_k) with h_j = (h_j1, ..., h_jn): h_ij in each dim-block
+  // Eigen::MatrixXd lvol_record = Eigen::MatrixXd::Zero(num_design * (num_iter + 1), dim); // time-varying h = (h_1, ..., h_k) with h_j = (h_j1, ..., h_jn): h_ij in each dim-block
+	Eigen::MatrixXd lvol_record = Eigen::MatrixXd::Zero(num_iter + 1, num_design * dim); // time-varying h = (h_1, ..., h_k) with h_j = (h_j1, ..., h_jn), row-binded
   // SSVS--------------
   Eigen::MatrixXd coef_dummy_record(num_iter + 1, num_alpha);
   Eigen::MatrixXd coef_weight_record(num_iter + 1, num_grp);
@@ -122,7 +123,8 @@ Rcpp::List estimate_var_sv(int num_iter, int num_burn,
   coef_record.row(0) = coefvec_ols;
   contem_coef_record.row(0) = Eigen::VectorXd::Zero(num_lowerchol); // initialize a as 0
   lvol_init_record.row(0) = (y - x * coef_mat).transpose().array().square().rowwise().mean().log(); // initialize h0 as mean of log((y - x alpha)^T (y - x alpha))
-  lvol_record.block(0, 0, num_design, dim) = lvol_init_record.row(0).replicate(num_design, 1);
+  // lvol_record.block(0, 0, num_design, dim) = lvol_init_record.row(0).replicate(num_design, 1);
+	Eigen::MatrixXd lvol_draw = lvol_init_record.row(0).replicate(num_design, 1); // h_j = (h_j1, ..., h_jn) for MCMC update
   lvol_sig_record.row(0) = .1 * Eigen::VectorXd::Ones(dim);
   // SSVS--------------
   coef_dummy_record.row(0) = Eigen::VectorXd::Ones(num_alpha);
@@ -136,7 +138,7 @@ Rcpp::List estimate_var_sv(int num_iter, int num_burn,
   Eigen::MatrixXd chol_lower = Eigen::MatrixXd::Zero(dim, dim); // L in Sig_t^(-1) = L D_t^(-1) LT
   Eigen::MatrixXd latent_innov(num_design, dim); // Z0 = Y0 - X0 A = (eps_p+1, eps_p+2, ..., eps_n+p)^T
   Eigen::MatrixXd ortho_latent(num_design, dim); // orthogonalized Z0
-  Eigen::MatrixXd lvol_draw = lvol_record.block(0, 0, num_design, dim);
+  // Eigen::MatrixXd lvol_draw = lvol_record.block(0, 0, num_design, dim);
   // Corrected triangular factorization-------
   Eigen::VectorXd prior_mean_j = Eigen::VectorXd::Zero(dim_design); // Prior mean vector of j-th column of A
   Eigen::MatrixXd prior_prec_j = Eigen::MatrixXd::Identity(dim_design, dim_design); // Prior precision of j-th column of A
@@ -305,7 +307,8 @@ Rcpp::List estimate_var_sv(int num_iter, int num_burn,
         ortho_latent.col(t)
       );
     }
-    lvol_record.block(num_design * i, 0, num_design, dim) = lvol_draw;
+    // lvol_record.block(num_design * i, 0, num_design, dim) = lvol_draw;
+		lvol_record.row(i) = vectorize_eigen(lvol_draw.transpose());
     // 3. a---------------------------------
     switch (prior_type) {
     case 2:
@@ -372,7 +375,7 @@ Rcpp::List estimate_var_sv(int num_iter, int num_burn,
   if (prior_type == 2) {
     return Rcpp::List::create(
       Rcpp::Named("alpha_record") = coef_record.bottomRows(num_iter - num_burn),
-      Rcpp::Named("h_record") = lvol_record,
+      Rcpp::Named("h_record") = lvol_record.bottomRows(num_iter - num_burn),
       Rcpp::Named("a_record") = contem_coef_record.bottomRows(num_iter - num_burn),
       Rcpp::Named("h0_record") = lvol_init_record.bottomRows(num_iter - num_burn),
       Rcpp::Named("sigh_record") = lvol_sig_record.bottomRows(num_iter - num_burn),
@@ -382,7 +385,7 @@ Rcpp::List estimate_var_sv(int num_iter, int num_burn,
     shrink_record.row(num_iter) = (Eigen::MatrixXd::Identity(num_coef, num_coef) + prior_alpha_prec).inverse().diagonal();
     return Rcpp::List::create(
       Rcpp::Named("alpha_record") = coef_record.bottomRows(num_iter - num_burn),
-      Rcpp::Named("h_record") = lvol_record,
+      Rcpp::Named("h_record") = lvol_record.bottomRows(num_iter - num_burn),
       Rcpp::Named("a_record") = contem_coef_record.bottomRows(num_iter - num_burn),
       Rcpp::Named("h0_record") = lvol_init_record.bottomRows(num_iter - num_burn),
       Rcpp::Named("sigh_record") = lvol_sig_record.bottomRows(num_iter - num_burn),
@@ -393,7 +396,7 @@ Rcpp::List estimate_var_sv(int num_iter, int num_burn,
   }
   return Rcpp::List::create(
     Rcpp::Named("alpha_record") = coef_record.bottomRows(num_iter - num_burn),
-    Rcpp::Named("h_record") = lvol_record,
+    Rcpp::Named("h_record") = lvol_record.bottomRows(num_iter - num_burn),
     Rcpp::Named("a_record") = contem_coef_record.bottomRows(num_iter - num_burn),
     Rcpp::Named("h0_record") = lvol_init_record.bottomRows(num_iter - num_burn),
     Rcpp::Named("sigh_record") = lvol_sig_record.bottomRows(num_iter - num_burn)

src/forecast-bvar.cpp

@@ -1,5 +1,6 @@
 #include <RcppEigen.h>
 #include "bvharmisc.h"
+#include "randsim.h"
 #include "bvhardraw.h"
 
 //' Forecasting BVAR(p)

src/forecast-bvhar.cpp

@@ -1,5 +1,6 @@
 #include <RcppEigen.h>
 #include "bvharmisc.h"
+#include "randsim.h"
 #include "bvhardraw.h"
 
 //' Forecasting Bayesian VHAR