diff --git a/R/RcppExports.R b/R/RcppExports.R
index 0f1e6f4..690e24d 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -428,6 +428,50 @@ post_coint_kls_sur <- function(y, beta, w, sigma_i, v_i, p_tau_i, g_i, x = NULL,
     .Call(`_bvartools_post_coint_kls_sur`, y, beta, w, sigma_i, v_i, p_tau_i, g_i, x, gamma_mu_prior, gamma_v_i_prior, svd)
 }
 
+#' Posterior Draws of Error Variances
+#' 
+#' Produces a draw of error variances from a gamma posterior density.
+#' 
+#' @param phi a \eqn{K \times T} matrix of time varying parameter draws.
+#' @param phi_init a \eqn{K \times 1} vector of initial states.
+#' @param shape_prior a \eqn{K \times 1} vector of prior shape parameters.
+#' @param rate_prior a \eqn{K \times 1} vector of prior rate parameters.
+#' 
+#' @details The function produces a posterior draw of the variaces vector \eqn{a} for the model
+#' Follow description in Chan eta al.
+#' 
+#' @references
+#' Chan, J., Koop, G., Poirier, D. J., & Tobias J. L. (2019). \emph{Bayesian econometric methods}
+#' (2nd ed.). Cambridge: Cambridge University Press.
+#' 
+#' @examples
+#' k <- 10 # Number of artificial coefficients
+#' tt <- 1000 # Number of observations
+#' 
+#' set.seed(1234) # Set RNG seed
+#' 
+#' # Generate artificial data according to a random walk
+#' phi <- matrix(rnorm(k), k, tt + 1)
+#' for (i in 2:(tt + 1)) {
+#'   phi[, i] <- phi[, i - 1] + rnorm(k, 0, sqrt(1 / 100))
+#' }
+#' 
+#' phi_init <- matrix(phi[, 1]) # Define inital state
+#' phi <- phi[, -1] # Drop initial state from main sample
+#' 
+#' # Define priors
+#' shape_prior <- matrix(1, k)
+#' rate_prior <- matrix(.0001, k)
+#' 
+#' # Obtain posterior draw
+#' post_gamma_state_variance(phi, phi_init, shape_prior, rate_prior)
+#' 
+#' @return A matrix.
+#' 
+post_gamma_state_variance <- function(phi, phi_init, shape_prior, rate_prior) {
+    .Call(`_bvartools_post_gamma_state_variance`, phi, phi_init, shape_prior, rate_prior)
+}
+
 #' Posterior Draw from a Normal Distribution
 #' 
 #' Produces a draw of coefficients from a normal posterior density.
diff --git a/man/post_gamma_state_variance.Rd b/man/post_gamma_state_variance.Rd
new file mode 100644
index 0000000..4bd385d
--- /dev/null
+++ b/man/post_gamma_state_variance.Rd
@@ -0,0 +1,54 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RcppExports.R
+\name{post_gamma_state_variance}
+\alias{post_gamma_state_variance}
+\title{Posterior Draws of Error Variances}
+\usage{
+post_gamma_state_variance(phi, phi_init, shape_prior, rate_prior)
+}
+\arguments{
+\item{phi}{a \eqn{K \times T} matrix of time varying parameter draws.}
+
+\item{phi_init}{a \eqn{K \times 1} vector of initial states.}
+
+\item{shape_prior}{a \eqn{K \times 1} vector of prior shape parameters.}
+
+\item{rate_prior}{a \eqn{K \times 1} vector of prior rate parameters.}
+}
+\value{
+A matrix.
+}
+\description{
+Produces a draw of error variances from a gamma posterior density.
+}
+\details{
+The function produces a posterior draw of the variaces vector \eqn{a} for the model
+Follow description in Chan eta al.
+}
+\examples{
+k <- 10 # Number of artificial coefficients
+tt <- 1000 # Number of observations
+
+set.seed(1234) # Set RNG seed
+
+# Generate artificial data according to a random walk
+phi <- matrix(rnorm(k), k, tt + 1)
+for (i in 2:(tt + 1)) {
+  phi[, i] <- phi[, i - 1] + rnorm(k, 0, sqrt(1 / 100))
+}
+
+phi_init <- matrix(phi[, 1]) # Define inital state
+phi <- phi[, -1] # Drop initial state from main sample
+
+# Define priors
+shape_prior <- matrix(1, k)
+rate_prior <- matrix(.0001, k)
+
+# Obtain posterior draw
+post_gamma_state_variance(phi, phi_init, shape_prior, rate_prior)
+
+}
+\references{
+Chan, J., Koop, G., Poirier, D. J., & Tobias J. L. (2019). \emph{Bayesian econometric methods}
+(2nd ed.). Cambridge: Cambridge University Press.
+}
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 5222404..e147392 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -212,6 +212,20 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// post_gamma_state_variance
+arma::mat post_gamma_state_variance(arma::mat phi, arma::vec phi_init, arma::vec shape_prior, arma::vec rate_prior);
+RcppExport SEXP _bvartools_post_gamma_state_variance(SEXP phiSEXP, SEXP phi_initSEXP, SEXP shape_priorSEXP, SEXP rate_priorSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< arma::mat >::type phi(phiSEXP);
+    Rcpp::traits::input_parameter< arma::vec >::type phi_init(phi_initSEXP);
+    Rcpp::traits::input_parameter< arma::vec >::type shape_prior(shape_priorSEXP);
+    Rcpp::traits::input_parameter< arma::vec >::type rate_prior(rate_priorSEXP);
+    rcpp_result_gen = Rcpp::wrap(post_gamma_state_variance(phi, phi_init, shape_prior, rate_prior));
+    return rcpp_result_gen;
+END_RCPP
+}
 // post_normal
 arma::vec post_normal(arma::mat y, arma::mat x, arma::mat sigma_i, arma::vec a_prior, arma::mat v_i_prior);
 RcppExport SEXP _bvartools_post_normal(SEXP ySEXP, SEXP xSEXP, SEXP sigma_iSEXP, SEXP a_priorSEXP, SEXP v_i_priorSEXP) {
@@ -436,6 +450,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_bvartools_loglik_normal", (DL_FUNC) &_bvartools_loglik_normal, 2},
     {"_bvartools_post_coint_kls", (DL_FUNC) &_bvartools_post_coint_kls, 10},
     {"_bvartools_post_coint_kls_sur", (DL_FUNC) &_bvartools_post_coint_kls_sur, 11},
+    {"_bvartools_post_gamma_state_variance", (DL_FUNC) &_bvartools_post_gamma_state_variance, 4},
     {"_bvartools_post_normal", (DL_FUNC) &_bvartools_post_normal, 5},
     {"_bvartools_post_normal_sur", (DL_FUNC) &_bvartools_post_normal_sur, 6},
     {"_bvartools_prep_covar_data", (DL_FUNC) &_bvartools_prep_covar_data, 4},
diff --git a/src/post_gamma_state_variance.cpp b/src/post_gamma_state_variance.cpp
new file mode 100644
index 0000000..5162d8a
--- /dev/null
+++ b/src/post_gamma_state_variance.cpp
@@ -0,0 +1,85 @@
+#include <RcppArmadillo.h>
+// [[Rcpp::depends(RcppArmadillo)]]
+
+//' Posterior Draws of Error Variances
+//' 
+//' Produces a draw of error variances from a gamma posterior density.
+//' 
+//' @param phi a \eqn{K \times T} matrix of time varying parameter draws.
+//' @param phi_init a \eqn{K \times 1} vector of initial states.
+//' @param shape_prior a \eqn{K \times 1} vector of prior shape parameters.
+//' @param rate_prior a \eqn{K \times 1} vector of prior rate parameters.
+//' 
+//' @details The function produces a posterior draw of the variaces vector \eqn{a} for the model
+//' Follow description in Chan eta al.
+//' 
+//' @references
+//' Chan, J., Koop, G., Poirier, D. J., & Tobias J. L. (2019). \emph{Bayesian econometric methods}
+//' (2nd ed.). Cambridge: Cambridge University Press.
+//' 
+//' @examples
+//' k <- 10 # Number of artificial coefficients
+//' tt <- 1000 # Number of observations
+//' 
+//' set.seed(1234) # Set RNG seed
+//' 
+//' # Generate artificial data according to a random walk
+//' phi <- matrix(rnorm(k), k, tt + 1)
+//' for (i in 2:(tt + 1)) {
+//'   phi[, i] <- phi[, i - 1] + rnorm(k, 0, sqrt(1 / 100))
+//' }
+//' 
+//' phi_init <- matrix(phi[, 1]) # Define inital state
+//' phi <- phi[, -1] # Drop initial state from main sample
+//' 
+//' # Define priors
+//' shape_prior <- matrix(1, k)
+//' rate_prior <- matrix(.0001, k)
+//' 
+//' # Obtain posterior draw
+//' post_gamma_state_variance(phi, phi_init, shape_prior, rate_prior)
+//' 
+//' @return A matrix.
+//' 
+// [[Rcpp::export]]
+arma::mat post_gamma_state_variance(arma::mat phi, arma::vec phi_init, arma::vec shape_prior, arma::vec rate_prior) {
+  
+  int k = phi.n_rows;
+  int tt = phi.n_cols;
+  arma::mat phi_lag = phi;
+  phi_lag.col(0) = phi_init;
+  phi_lag.cols(1, tt - 1) = phi.cols(0, tt - 2);
+  arma::mat phi_v = arma::trans(phi - phi_lag);
+  arma::vec psi_sigma_v_post_scale = 1 / (rate_prior + arma::vectorise(arma::sum(arma::pow(phi_v, 2))) * 0.5);
+  arma::mat psi_sigma_i = arma::zeros<arma::mat>(k, k);
+  arma::vec shape_post = shape_prior + tt * 0.5;
+  for (int i = 0; i < k; i++) {
+    psi_sigma_i(i, i) = arma::randg<double>(arma::distr_param(shape_post(i), psi_sigma_v_post_scale(i)));
+  }
+  
+  return psi_sigma_i;
+}
+
+/*** R
+k <- 10 # Number of artificial coefficients
+tt <- 1000 # Number of observations
+
+set.seed(1234) # Set RNG seed
+
+# Generate artificial data according to a random walk
+phi <- matrix(rnorm(k), k, tt + 1)
+for (i in 2:(tt + 1)) {
+  phi[, i] <- phi[, i - 1] + rnorm(k, 0, sqrt(1 / 100))
+}
+
+phi_init <- matrix(phi[, 1]) # Define inital state
+phi <- phi[, -1] # Drop initial state from main sample
+
+# Define priors
+shape_prior <- matrix(1, k)
+rate_prior <- matrix(.0001, k)
+
+# Obtain posterior draw
+post_gamma_state_variance(phi, phi_init, shape_prior, rate_prior)
+
+*/
\ No newline at end of file