diff --git a/R/RcppExports.R b/R/RcppExports.R
index 690e24d..5fa2cc0 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -430,15 +430,23 @@ post_coint_kls_sur <- function(y, beta, w, sigma_i, v_i, p_tau_i, g_i, x = NULL,
 
 #' Posterior Draws of Error Variances
 #' 
-#' Produces a draw of error variances from a gamma posterior density.
+#' Produces a draw of the constant diagonal error variance matrix of the
+#' state equation of a state space model using an inverse 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 a a \eqn{K \times T} matrix of time varying parameter draws.
+#' @param a_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.
+#' @param inverse logical. If \code{TRUE}, the function returns the precision matrix,
+#' i.e. the inverse of the variance matrix. Defaults to \code{FALSE}.
 #' 
-#' @details The function produces a posterior draw of the variaces vector \eqn{a} for the model
-#' Follow description in Chan eta al.
+#' @details For the state space model with state equation
+#' \deqn{a_t = a_{t-1} + v}
+#' and measurement equation
+#' \deqn{y_t = Z_{t} a_t + u_t}
+#' with \eqn{v_t \sim N(0, \Sigma_{v})} and \eqn{u_t \sim N(0, \Sigma_{u,t})}
+#' the function produces a draw of the constant diagonal error variances matrix of the
+#' state equation \eqn{\Simga_v}.
 #' 
 #' @references
 #' Chan, J., Koop, G., Poirier, D. J., & Tobias J. L. (2019). \emph{Bayesian econometric methods}
@@ -451,25 +459,25 @@ post_coint_kls_sur <- function(y, beta, w, sigma_i, v_i, p_tau_i, g_i, x = NULL,
 #' set.seed(1234) # Set RNG seed
 #' 
 #' # Generate artificial data according to a random walk
-#' phi <- matrix(rnorm(k), k, tt + 1)
+#' a <- matrix(rnorm(k), k, tt + 1)
 #' for (i in 2:(tt + 1)) {
-#'   phi[, i] <- phi[, i - 1] + rnorm(k, 0, sqrt(1 / 100))
+#'   a[, i] <- a[, 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
+#' a_init <- matrix(a[, 1]) # Define inital state
+#' a <- a[, -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)
+#' post_gamma_state_variance(a, a_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)
+post_gamma_state_variance <- function(a, a_init, shape_prior, rate_prior, inverse = FALSE) {
+    .Call(`_bvartools_post_gamma_state_variance`, a, a_init, shape_prior, rate_prior, inverse)
 }
 
 #' Posterior Draw from a Normal Distribution
diff --git a/man/post_gamma_state_variance.Rd b/man/post_gamma_state_variance.Rd
index 4bd385d..aa95682 100644
--- a/man/post_gamma_state_variance.Rd
+++ b/man/post_gamma_state_variance.Rd
@@ -4,26 +4,35 @@
 \alias{post_gamma_state_variance}
 \title{Posterior Draws of Error Variances}
 \usage{
-post_gamma_state_variance(phi, phi_init, shape_prior, rate_prior)
+post_gamma_state_variance(a, a_init, shape_prior, rate_prior, inverse = FALSE)
 }
 \arguments{
-\item{phi}{a \eqn{K \times T} matrix of time varying parameter draws.}
+\item{a}{a \eqn{K \times T} matrix of time varying parameter draws.}
 
-\item{phi_init}{a \eqn{K \times 1} vector of initial states.}
+\item{a_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.}
+
+\item{inverse}{logical. If \code{TRUE}, the function returns the precision matrix,
+i.e. the inverse of the variance matrix. Defaults to \code{FALSE}.}
 }
 \value{
 A matrix.
 }
 \description{
-Produces a draw of error variances from a gamma posterior density.
+Produces a draw of the constant diagonal error variance matrix of the
+state equation of a state space model using an inverse 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.
+For the state space model with state equation
+\deqn{a_t = a_{t-1} + v}
+and measurement equation
+\deqn{y_t = Z_{t} a_t + u_t}
+with \eqn{v_t \sim N(0, \Sigma_{v})} and \eqn{u_t \sim N(0, \Sigma_{u,t})}
+the function produces a draw of the constant diagonal error variances matrix of the
+state equation \eqn{\Simga_v}.
 }
 \examples{
 k <- 10 # Number of artificial coefficients
@@ -32,20 +41,20 @@ 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)
+a <- matrix(rnorm(k), k, tt + 1)
 for (i in 2:(tt + 1)) {
-  phi[, i] <- phi[, i - 1] + rnorm(k, 0, sqrt(1 / 100))
+  a[, i] <- a[, 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
+a_init <- matrix(a[, 1]) # Define inital state
+a <- a[, -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)
+post_gamma_state_variance(a, a_init, shape_prior, rate_prior)
 
 }
 \references{
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index e147392..65a1cdf 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -213,16 +213,17 @@ BEGIN_RCPP
 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) {
+arma::mat post_gamma_state_variance(arma::mat a, arma::vec a_init, arma::vec shape_prior, arma::vec rate_prior, bool inverse);
+RcppExport SEXP _bvartools_post_gamma_state_variance(SEXP aSEXP, SEXP a_initSEXP, SEXP shape_priorSEXP, SEXP rate_priorSEXP, SEXP inverseSEXP) {
 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::mat >::type a(aSEXP);
+    Rcpp::traits::input_parameter< arma::vec >::type a_init(a_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));
+    Rcpp::traits::input_parameter< bool >::type inverse(inverseSEXP);
+    rcpp_result_gen = Rcpp::wrap(post_gamma_state_variance(a, a_init, shape_prior, rate_prior, inverse));
     return rcpp_result_gen;
 END_RCPP
 }
@@ -450,7 +451,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_gamma_state_variance", (DL_FUNC) &_bvartools_post_gamma_state_variance, 5},
     {"_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
index 5162d8a..7e420a7 100644
--- a/src/post_gamma_state_variance.cpp
+++ b/src/post_gamma_state_variance.cpp
@@ -3,15 +3,23 @@
 
 //' Posterior Draws of Error Variances
 //' 
-//' Produces a draw of error variances from a gamma posterior density.
+//' Produces a draw of the constant diagonal error variance matrix of the
+//' state equation of a state space model using an inverse 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 a a \eqn{K \times T} matrix of time varying parameter draws.
+//' @param a_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.
+//' @param inverse logical. If \code{TRUE}, the function returns the precision matrix,
+//' i.e. the inverse of the variance matrix. Defaults to \code{FALSE}.
 //' 
-//' @details The function produces a posterior draw of the variaces vector \eqn{a} for the model
-//' Follow description in Chan eta al.
+//' @details For the state space model with state equation
+//' \deqn{a_t = a_{t-1} + v}
+//' and measurement equation
+//' \deqn{y_t = Z_{t} a_t + u_t}
+//' with \eqn{v_t \sim N(0, \Sigma_{v})} and \eqn{u_t \sim N(0, \Sigma_{u,t})}
+//' the function produces a draw of the constant diagonal error variances matrix of the
+//' state equation \eqn{\Simga_v}.
 //' 
 //' @references
 //' Chan, J., Koop, G., Poirier, D. J., & Tobias J. L. (2019). \emph{Bayesian econometric methods}
@@ -24,40 +32,43 @@
 //' set.seed(1234) # Set RNG seed
 //' 
 //' # Generate artificial data according to a random walk
-//' phi <- matrix(rnorm(k), k, tt + 1)
+//' a <- matrix(rnorm(k), k, tt + 1)
 //' for (i in 2:(tt + 1)) {
-//'   phi[, i] <- phi[, i - 1] + rnorm(k, 0, sqrt(1 / 100))
+//'   a[, i] <- a[, 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
+//' a_init <- matrix(a[, 1]) # Define inital state
+//' a <- a[, -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)
+//' post_gamma_state_variance(a, a_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) {
+arma::mat post_gamma_state_variance(arma::mat a, arma::vec a_init, arma::vec shape_prior, arma::vec rate_prior, bool inverse = false) {
   
-  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);
+  int k = a.n_rows;
+  int tt = a.n_cols;
+  arma::mat a_lag = a;
+  a_lag.col(0) = a_init;
+  a_lag.cols(1, tt - 1) = a.cols(0, tt - 2);
+  arma::mat a_v = arma::trans(a - a_lag);
+  arma::vec psi_sigma_v_post_scale = 1 / (rate_prior + arma::vectorise(arma::sum(arma::pow(a_v, 2))) * 0.5);
+  arma::mat psi_sigma = 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)));
+    psi_sigma(i, i) = arma::randg<double>(arma::distr_param(shape_post(i), psi_sigma_v_post_scale(i)));
+    if (inverse) {
+      psi_sigma(i, i) = 1 / psi_sigma(i, i);
+    }
   }
   
-  return psi_sigma_i;
+  return psi_sigma;
 }
 
 /*** R
@@ -67,19 +78,19 @@ 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)
+a <- matrix(rnorm(k), k, tt + 1)
 for (i in 2:(tt + 1)) {
-  phi[, i] <- phi[, i - 1] + rnorm(k, 0, sqrt(1 / 100))
+  a[, i] <- a[, 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
+a_init <- matrix(a[, 1]) # Define inital state
+a <- a[, -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)
+post_gamma_state_variance(a, a_init, shape_prior, rate_prior, inverse = FALSE)
 
 */
\ No newline at end of file