diff --git a/R/RcppExports.R b/R/RcppExports.R
index c6d848f1..f669537c 100644
--- a/R/RcppExports.R
+++ b/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
diff --git a/R/bvar-sv.R b/R/bvar-sv.R
index 4ad1bda9..b7382ea3 100644
--- a/R/bvar-sv.R
+++ b/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)
diff --git a/R/bvhar-sv.R b/R/bvhar-sv.R
index f6db89d7..d2ac7a72 100644
--- a/R/bvhar-sv.R
+++ b/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)
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index ff3362db..26b526a1 100644
--- a/src/RcppExports.cpp
+++ b/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},
diff --git a/src/bvhardraw.h b/src/bvhardraw.h
index aecbd0e0..0947c5a0 100644
--- a/src/bvhardraw.h
+++ b/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);
 
diff --git a/src/bvharmisc.h b/src/bvharmisc.h
index c73b665d..9d10c381 100644
--- a/src/bvharmisc.h
+++ b/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);
diff --git a/src/estimate-hierarchical.cpp b/src/estimate-hierarchical.cpp
index 5740b8c0..b9512699 100644
--- a/src/estimate-hierarchical.cpp
+++ b/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>
 
diff --git a/src/estimate-sv.cpp b/src/estimate-sv.cpp
index 6efa0166..dce7a02b 100644
--- a/src/estimate-sv.cpp
+++ b/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)
diff --git a/src/forecast-bvar.cpp b/src/forecast-bvar.cpp
index 10de049c..c30f1864 100644
--- a/src/forecast-bvar.cpp
+++ b/src/forecast-bvar.cpp
@@ -1,5 +1,6 @@
 #include <RcppEigen.h>
 #include "bvharmisc.h"
+#include "randsim.h"
 #include "bvhardraw.h"
 
 //' Forecasting BVAR(p)
diff --git a/src/forecast-bvhar.cpp b/src/forecast-bvhar.cpp
index 26a07913..117d9aa8 100644
--- a/src/forecast-bvhar.cpp
+++ b/src/forecast-bvhar.cpp
@@ -1,5 +1,6 @@
 #include <RcppEigen.h>
 #include "bvharmisc.h"
+#include "randsim.h"
 #include "bvhardraw.h"
 
 //' Forecasting Bayesian VHAR
diff --git a/src/generate-distn.cpp b/src/generate-distn.cpp
index eabea3f2..9a8d4034 100644
--- a/src/generate-distn.cpp
+++ b/src/generate-distn.cpp
@@ -311,3 +311,173 @@ Eigen::MatrixXd sim_wishart(Eigen::MatrixXd mat_scale, double shape) {
   Eigen::MatrixXd chol_res = chol_scale * mat_bartlett;
   return chol_res * chol_res.transpose();
 }
+
+//' 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
+// [[Rcpp::export]]
+double dgig_quasi(double x, double lambda, double beta) {
+	return pow(x, lambda - 1) * exp(-beta * (x + 1 / x) / 2);
+}
+
+//' 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
+// [[Rcpp::export]]
+Eigen::VectorXd rgig_nonconcave(int num_sim, double lambda, double beta) {
+	Eigen::VectorXd res(num_sim);
+	double mode = beta / (sqrt((1 - lambda) * (1 - lambda) + beta * beta) + 1 - lambda); // argmax of g(x)
+	double x0, xstar, k1, k2, k3, A1, A2, A3;
+	x0 = beta / (1 - beta); // subdomain (0, x0)
+	xstar = std::max(x0, 2 / beta);
+	k1 = dgig_quasi(mode, lambda, beta);
+	A1 = k1 * x0;
+	if (x0 < 2 / beta) { // subdomain (x0, 2 / beta)
+		k2 = exp(-beta);
+		if (lambda == 0) {
+			A2 = k2 * log(2 / (beta * beta));
+		} else {
+			A2 = k2 * (pow(2 / beta, lambda) - pow(x0, lambda)) / lambda;
+		}
+	} else {
+		k2 = 0;
+		A2 = 0;
+	}
+	k3 = pow(xstar, lambda - 1);  // subdomain (xstar, inf)
+	A3 = 2 * k3 * exp(-xstar * beta / 2) / beta;
+	double A = A1 + A2 + A3;
+	bool rejected;
+	double draw_unif, draw_prop, cand, ar_const;
+	// double draw_x, draw_y, cand; // bounded rectangle
+	for (int i = 0; i < num_sim; i++) {
+		rejected = true;
+		while (rejected) {
+			draw_unif = unif_rand(0, 1);
+			draw_prop = unif_rand(0, A);
+			if (draw_prop <= A1) {
+				cand = x0 * draw_prop / A1;
+				ar_const = k1;
+			} else if (draw_prop <= A1 + A2) {
+				draw_prop -= A1;
+				cand = pow(pow(x0, lambda) + draw_prop * lambda / k2, 1 / lambda);
+				ar_const = k2 * pow(cand, lambda - 1);
+			} else {
+				draw_prop -= A1 + A2;
+				cand = -2 * log(exp(-xstar * beta / 2) - draw_prop * beta / (2 * k3));
+				ar_const = k3 * exp(-cand * beta / 2);
+			}
+			rejected = draw_unif * ar_const > dgig_quasi(cand, lambda, beta);
+		}
+		res[i] = cand;
+	}
+	return res;
+}
+
+//' 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
+// [[Rcpp::export]]
+Eigen::VectorXd rgig_without_mode(int num_sim, double lambda, double beta) {
+	Eigen::VectorXd res(num_sim);
+	double arg_y = beta / (sqrt((1 - lambda) * (1 - lambda) + beta * beta) + 1 - lambda); // argmax of g(x)
+	double arg_x = (1 + lambda + sqrt((1 + lambda) * (1 + lambda) + beta * beta)) / beta; // argmax of x g(x)
+	double bound_y = sqrt(dgig_quasi(arg_y, lambda, beta)); // max sqrt(g(x))
+	double bound_x = arg_x * sqrt(dgig_quasi(arg_x, lambda, beta)); // max x*sqrt(g(x))
+	bool rejected;
+	double draw_x, draw_y, cand; // bounded rectangle
+	for (int i = 0; i < num_sim; i++) {
+		rejected = true;
+		while (rejected) {
+			draw_x = unif_rand(0, bound_x);
+			draw_y = unif_rand(0, bound_y);
+			cand = draw_x / draw_y;
+			rejected = draw_y * draw_y > dgig_quasi(cand, lambda, beta); // Check if U <= g(y) / unif(y)
+		}
+		res[i] = cand;
+	}
+	return res;
+}
+
+//' 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
+// [[Rcpp::export]]
+Eigen::VectorXd rgig_with_mode(int num_sim, double lambda, double beta) {
+	Eigen::VectorXd res(num_sim);
+	double arg_y = (sqrt((1 - lambda) * (1 - lambda) + beta * beta) - 1 + lambda) / beta; // argmax of g(x)
+	double quad_coef = -2 * (lambda + 1) / beta - arg_y;
+	double lin_coef = 2 * arg_y * (lambda - 1) / beta - 1;
+	double p = lin_coef - quad_coef * quad_coef / 3;
+	double q = 2 * quad_coef * quad_coef * quad_coef / 27 - quad_coef * lin_coef * arg_y / 3 + arg_y;
+	double phi = acos(-q * sqrt(-27 / (p * p * p)) / 2);
+	double arg_x_neg = sqrt(-p * 4 / 3) * cos(phi / 3 + M_PI * 4 / 3) - quad_coef / 3;
+	double arg_x_pos = sqrt(-p * 4 / 3) * cos(phi / 3) - quad_coef / 3;
+	double bound_y = sqrt(dgig_quasi(arg_y, lambda, beta));
+	double bound_x_neg = (arg_x_neg - arg_y) * sqrt(dgig_quasi(arg_x_neg, lambda, beta));
+	double bound_x_pos = (arg_x_pos - arg_y) * sqrt(dgig_quasi(arg_x_pos, lambda, beta));
+	bool rejected;
+	double draw_x, draw_y, cand; // bounded rectangle
+	for (int i = 0; i < num_sim; i++) {
+		rejected = true;
+		while (rejected) {
+			draw_x = unif_rand(bound_x_neg, bound_x_pos);
+			draw_y = unif_rand(0, bound_y);
+			cand = draw_x / draw_y + arg_y;
+			rejected = draw_y * draw_y > dgig_quasi(cand, lambda, beta); // Check if U <= g(y) / unif(y)
+		}
+		res[i] = cand;
+	}
+	return res;
+}
+
+//' 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
+// [[Rcpp::export]]
+Eigen::VectorXd sim_gig(int num_sim, double lambda, double psi, double chi) {
+	if (psi <= 0 || chi <= 0) {
+		Rcpp::stop("Wrong 'psi' and 'chi' range.");
+	}
+	Eigen::VectorXd res(num_sim);
+	double abs_lam = abs(lambda); // If lambda < 0, use 1 / X as the result
+	double alpha = sqrt(psi / chi); // scaling parameter of quasi-density: scale the result
+	double beta = sqrt(psi * chi); // second parameter of quasi-density
+	double quasi_bound = sqrt(1 - lambda) * 2 / 3;
+	if (abs_lam < 1 && beta <= quasi_bound) {
+		res = rgig_nonconcave(num_sim, lambda, beta) * alpha; // non-T_(-1/2)-concave part
+	} else if (abs_lam <= 1 && beta >= std::min(.5, quasi_bound) && beta <= 1) {
+		res = rgig_without_mode(num_sim, lambda, beta) * alpha; // without mode shift
+	} else if (abs_lam > 1 && beta > 1) {
+		res = rgig_with_mode(num_sim, lambda, beta) * alpha; // with mode shift
+	} else {
+		Rcpp::stop("Wrong parameter ranges for quasi GIG density.");
+	}
+	if (lambda < 0) {
+		res = res.cwiseInverse();
+	}
+	return res;
+}
diff --git a/src/mcmc-draw.cpp b/src/mcmc-draw.cpp
index 83e74b7d..3bf4d59d 100644
--- a/src/mcmc-draw.cpp
+++ b/src/mcmc-draw.cpp
@@ -1,7 +1,7 @@
 #include "bvharomp.h"
 #include <RcppEigen.h>
 #include "bvharmisc.h"
-#include "bvharprob.h"
+#include "randsim.h"
 
 //' Building Spike-and-slab SD Diagonal Matrix
 //' 
diff --git a/src/randsim.h b/src/randsim.h
new file mode 100644
index 00000000..273656bf
--- /dev/null
+++ b/src/randsim.h
@@ -0,0 +1,22 @@
+#ifndef RANDSIM_H
+#define RANDSIM_H
+
+#include "bvharprob.h"
+
+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 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::VectorXd sim_gig(int num_sim, double lambda, double psi, double chi);
+
+#endif // RANDSIM_H
diff --git a/src/simulate-data.cpp b/src/simulate-data.cpp
index b3f3fb9b..8d309a1b 100644
--- a/src/simulate-data.cpp
+++ b/src/simulate-data.cpp
@@ -1,5 +1,6 @@
 #include <RcppEigen.h>
 #include "bvharmisc.h"
+#include "randsim.h"
 
 //' Generate Multivariate Time Series Process Following VAR(p)
 //'