included a centred SV process closes #22

bsvars · Nov 25, 2022 · 66ea974 · 66ea974
1 parent 917c66a
commit 66ea974
Show file tree

Hide file tree

Showing 20 changed files with 412 additions and 72 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -20,4 +20,4 @@ Suggests:
     tinytest
 LinkingTo: Rcpp, RcppProgress, RcppArmadillo, RcppTN
 BugReports: https://github.com/donotdespair/bsvars/issues
-RoxygenNote: 7.2.1
+RoxygenNote: 7.2.2
diff --git a/NEWS.md b/NEWS.md
@@ -15,7 +15,7 @@
 6.  Include citation info for the package [#12](https://github.com/donotdespair/bsvars/issues/12)
 7.  Corrected sampler for AR parameter of the SV equations [#19](https://github.com/donotdespair/bsvars/issues/19)
 8.  Added samplers from joint predictive densities [#15](https://github.com/donotdespair/bsvars/issues/15)
-
+9.  A new centred Stochastic Volatility heteroskedastic process is implemented [#22](https://github.com/donotdespair/bsvars/issues/22)
 
 # bsvars 1.0.0
 

diff --git a/R/bsvars-package.R b/R/bsvars-package.R
@@ -56,7 +56,8 @@
 #' \itemize{
 #'   \item homoskedastic model with unit variances
 #'   \item heteroskedastic model with stationary Markov switching in the variances
-#'   \item heteroskedastic model with Stochastic Volatility process for variances
+#'   \item heteroskedastic model with non-centred Stochastic Volatility process for variances
+#'   \item heteroskedastic model with centred Stochastic Volatility process for variances
 #'   \item non-normal model with a finite mixture of normal components and component-specific variances
 #'   \item heteroskedastic model with sparse Markov switching in the variances where the number of heteroskedastic components is estimated
 #'   \item non-normal model with a sparse mixture of normal components and component-specific variances where the number of heteroskedastic components is estimated

diff --git a/R/estimate.BSVARSV.R b/R/estimate.BSVARSV.R
@@ -2,13 +2,17 @@
 #' @title Bayesian estimation of a Structural Vector Autoregression with 
 #' Stochastic Volatility heteroskedasticity via Gibbs sampler
 #'
-#' @description Estimates the SVAR with Stochastic Volatility (SV) heteroskedasticity proposed by Lütkepohl, Shang, Uzeda, and Woźniak (2022).
+#' @description Estimates the SVAR with Stochastic Volatility (SV) heteroskedasticity 
+#' proposed by Lütkepohl, Shang, Uzeda, and Woźniak (2022).
 #' Implements the Gibbs sampler proposed by Waggoner & Zha (2003)
-#' for the structural matrix \eqn{B} and the equation-by-equation sampler by Chan, Koop, & Yu (2021)
-#' for the autoregressive slope parameters \eqn{A}. Additionally, the parameter matrices \eqn{A} and \eqn{B}
-#' follow a Minnesota prior and generalised-normal prior distributions respectively with the matrix-specific
-#' overall shrinkage parameters estimated thanks to a 3-level hierarchical prior distribution. The SV model 
-#' is estimated in a non-centred parameterisation using a range of techniques including: 
+#' for the structural matrix \eqn{B} and the equation-by-equation sampler 
+#' by Chan, Koop, & Yu (2021)
+#' for the autoregressive slope parameters \eqn{A}. Additionally, 
+#' the parameter matrices \eqn{A} and \eqn{B}
+#' follow a Minnesota prior and generalised-normal prior distributions 
+#' respectively with the matrix-specific
+#' overall shrinkage parameters estimated thanks to a 3-level hierarchical prior distribution. 
+#' The SV model is estimated using a range of techniques including: 
 #' simulation smoother, auxiliary mixture, ancillarity-sufficiency interweaving strategy, 
 #' and generalised inverse Gaussian distribution summarised by Kastner & Frühwirth-Schnatter (2014). 
 #' See section \bold{Details} for the model equations.
@@ -23,15 +27,30 @@
 #' \deqn{BE = U}
 #' where \eqn{U} is an \code{NxT} matrix of structural form error terms, and
 #' \eqn{B} is an \code{NxN} matrix of contemporaneous relationships.
-#' 
 #' Finally, the structural shocks, \eqn{U}, are temporally and contemporaneously independent and jointly normally distributed with zero mean.
-#' The conditional variance of the \code{n}th shock at time \code{t} is given by:
+#' 
+#' Two alternative specifications of the conditional variance of the \code{n}th shock at time \code{t} 
+#' can be estimated: non-centred Stochastic Volatility by Lütkepohl, Shang, Uzeda, and Woźniak (2022) 
+#' or centred Stochastic Volatility by Chan, Koop, & Yu (2021).
+#' 
+#' The non-centred Stochastic Volatility by Lütkepohl, Shang, Uzeda, and Woźniak (2022) 
+#' is selected by setting argument \code{centred_sv} of function \code{specify_bsvar_sv$new()} to value \code{FALSE}.
+#' It has the conditional variances given by:
 #' \deqn{Var_{t-1}[u_{n.t}] = exp(w_n h_{n.t})}
 #' where \eqn{w_n} is the estimated conditional standard deviation of the log-conditional variance
 #' and the log-volatility process \eqn{h_{n.t}} follows an autoregressive process:
 #' \deqn{h_{n.t} = g_n h_{n.t-1} + v_{n.t}}
 #' where \eqn{h_{n.0}=0}, \eqn{g_n} is an autoregressive parameter and \eqn{v_{n.t}} is a standard normal error term.
 #' 
+#' The centred Stochastic Volatility by Chan, Koop, & Yu (2021)
+#' is selected by setting argument \code{centred_sv} of function \code{specify_bsvar_sv$new()} to value \code{TRUE}.
+#' Its conditional variances are given by:
+#' \deqn{Var_{t-1}[u_{n.t}] = exp(h_{n.t})}
+#' where the log-conditional variances \eqn{h_{n.t}} follow an autoregressive process:
+#' \deqn{h_{n.t} = g_n h_{n.t-1} + v_{n.t}}
+#' where \eqn{h_{n.0}=0}, \eqn{g_n} is an autoregressive parameter and \eqn{v_{n.t}} is a zero-mean normal error term
+#' with variance \eqn{s_{v.n}^2}.
+#' 
 #' @param specification an object of class BSVARSV generated using the \code{specify_bsvar_sv$new()} function.
 #' @param S a positive integer, the number of posterior draws to be generated
 #' @param thin a positive integer, specifying the frequency of MCMC output thinning
@@ -112,9 +131,10 @@ estimate.BSVARSV <- function(specification, S, thin = 10, show_progress = TRUE)
   starting_values     = specification$starting_values$get_starting_values()
   VB                  = specification$identification$get_identification()
   data_matrices       = specification$data_matrices$get_data_matrices()
+  centred_sv          = specification$centred_sv
 
   # estimation
-  qqq                 = .Call(`_bsvars_bsvar_sv_cpp`, S, data_matrices$Y, data_matrices$X, prior, VB, starting_values, thin, TRUE, show_progress)
+  qqq                 = .Call(`_bsvars_bsvar_sv_cpp`, S, data_matrices$Y, data_matrices$X, prior, VB, starting_values, thin, centred_sv, show_progress)
 
   specification$starting_values$set_starting_values(qqq$last_draw)
   output              = specify_posterior_bsvar_sv$new(specification, qqq$posterior)
@@ -171,9 +191,10 @@ estimate.PosteriorBSVARSV <- function(specification, S, thin = 10, show_progress
   starting_values     = specification$last_draw$starting_values$get_starting_values()
   VB                  = specification$last_draw$identification$get_identification()
   data_matrices       = specification$last_draw$data_matrices$get_data_matrices()
+  centred_sv          = specification$last_draw$centred_sv
 
   # estimation
-  qqq                 = .Call(`_bsvars_bsvar_sv_cpp`, S, data_matrices$Y, data_matrices$X, prior, VB, starting_values, thin, TRUE, show_progress)
+  qqq                 = .Call(`_bsvars_bsvar_sv_cpp`, S, data_matrices$Y, data_matrices$X, prior, VB, starting_values, thin, centred_sv, show_progress)
 
   specification$last_draw$starting_values$set_starting_values(qqq$last_draw)
   output              = specify_posterior_bsvar_sv$new(specification$last_draw, qqq$posterior)

diff --git a/R/forecast.R b/R/forecast.R
@@ -227,8 +227,9 @@ forecast.PosteriorBSVARSV = function(posterior, horizon) {
   T                 = ncol(posterior$last_draw$data_matrices$X)
   X_T               = posterior$last_draw$data_matrices$X[,T]
   posterior_h_T     = posterior$posterior$h[,T,]
+  centred_sv        = posterior$last_draw$centred_sv
 
-  fore            = .Call(`_bsvars_forecast_bsvar_sv`, posterior_B, posterior_A, posterior_h_T, posterior_rho, posterior_omega, X_T, horizon)
+  fore            = .Call(`_bsvars_forecast_bsvar_sv`, posterior_B, posterior_A, posterior_h_T, posterior_rho, posterior_omega, X_T, horizon, centred_sv)
   class(fore)     = "Forecasts"
 
   return(fore)

diff --git a/R/specify_bsvar_sv.R b/R/specify_bsvar_sv.R
@@ -125,6 +125,9 @@ specify_starting_values_bsvar_sv = R6::R6Class(
     #' @field omega an \code{N}-vector with values of SV process conditional standard deviations.
     omega         = numeric(),
 
+    #' @field sigma2v an \code{N}-vector with values of SV process conditional variances.
+    sigma2v       = numeric(),
+
     #' @field S an \code{NxT} integer matrix with the auxiliary mixture component indicators.
     S             = matrix(),
 
@@ -150,6 +153,7 @@ specify_starting_values_bsvar_sv = R6::R6Class(
       self$h              = matrix(rnorm(N * T, sd = .01), N, T)
       self$rho            = rep(.5, N)
       self$omega          = rep(.1, N)
+      self$sigma2v        = rep(.1^2, N)
       self$S              = matrix(1, N, T)
       self$sigma2_omega   = rep(1, N)
       self$s_             = rep(0.05, N)
@@ -171,6 +175,7 @@ specify_starting_values_bsvar_sv = R6::R6Class(
         h                 = self$h,
         rho               = self$rho,
         omega             = self$omega,
+        sigma2v           = self$sigma2v,
         S                 = self$S,
         sigma2_omega      = self$sigma2_omega,
         s_                = self$s_
@@ -198,6 +203,7 @@ specify_starting_values_bsvar_sv = R6::R6Class(
       self$h              = last_draw$h
       self$rho            = last_draw$rho
       self$omega          = last_draw$omega
+      self$sigma2v        = last_draw$sigma2v
       self$S              = last_draw$S
       self$sigma2_omega   = last_draw$sigma2_omega
       self$s_             = last_draw$s_
@@ -242,17 +248,22 @@ specify_bsvar_sv = R6::R6Class(
     #' @field starting_values an object StartingValuesBSVARSV with the starting values.
     starting_values        = list(),
 
+    #' @field centred_sv a logical value - if true a centred parameterisation of the Stochastic Volatility process is estimated. Otherwise, its non-centred parameterisation is estimated. See Lütkepohl, Shang, Uzeda, Woźniak (2022) for more info.
+    centred_sv             = logical(),
+
     #' @description
     #' Create a new specification of the BSVAR model with Stochastic Volatility heteroskedasticity, BSVARSV.
     #' @param data a \code{(T+p)xN} matrix with time series data.
     #' @param p a positive integer providing model's autoregressive lag order.
     #' @param B a logical \code{NxN} matrix containing value \code{TRUE} for the elements of the structural matrix \eqn{B} to be estimated and value \code{FALSE} for exclusion restrictions to be set to zero.
+    #' @param centred_sv a logical value. If \code{FALSE} a non-centred Stochastic Volatility processes for conditional variances are estimated. Otherwise, a centred process is estimated.
     #' @param stationary an \code{N} logical vector - its element set to \code{FALSE} sets the prior mean for the autoregressive parameters of the \code{N}th equation to the white noise process, otherwise to random walk.
     #' @return A new complete specification for the bsvar model with Stochastic Volatility heteroskedasticity, BSVARSV.
     initialize = function(
     data,
     p = 1L,
     B,
+    centred_sv = FALSE,
     stationary = rep(FALSE, ncol(data))
     ) {
       stopifnot("Argument p has to be a positive integer." = ((p %% 1) == 0 & p > 0))
@@ -273,6 +284,7 @@ specify_bsvar_sv = R6::R6Class(
       self$identification  = specify_identification_bsvars$new(N, B)
       self$prior           = specify_prior_bsvar_sv$new(N, p, stationary)
       self$starting_values = specify_starting_values_bsvar_sv$new(N, self$p, T)
+      self$centred_sv      = centred_sv
     }, # END initialize
 
     #' @description

diff --git a/inst/include/bsvars_RcppExports.h b/inst/include/bsvars_RcppExports.h
@@ -193,7 +193,7 @@ namespace bsvars {
         return Rcpp::as<Rcpp::List >(rcpp_result_gen);
     }
 
-    inline Rcpp::List bsvar_sv_cpp(const int& S, const arma::mat& Y, const arma::mat& X, const Rcpp::List& prior, const arma::field<arma::mat>& VB, const Rcpp::List& starting_values, const int thin = 100, const bool sample_s_ = true, const bool show_progress = true) {
+    inline Rcpp::List bsvar_sv_cpp(const int& S, const arma::mat& Y, const arma::mat& X, const Rcpp::List& prior, const arma::field<arma::mat>& VB, const Rcpp::List& starting_values, const int thin = 100, const bool centred_sv = false, const bool show_progress = true) {
         typedef SEXP(*Ptr_bsvar_sv_cpp)(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP);
         static Ptr_bsvar_sv_cpp p_bsvar_sv_cpp = NULL;
         if (p_bsvar_sv_cpp == NULL) {
@@ -203,7 +203,7 @@ namespace bsvars {
         RObject rcpp_result_gen;
         {
             RNGScope RCPP_rngScope_gen;
-            rcpp_result_gen = p_bsvar_sv_cpp(Shield<SEXP>(Rcpp::wrap(S)), Shield<SEXP>(Rcpp::wrap(Y)), Shield<SEXP>(Rcpp::wrap(X)), Shield<SEXP>(Rcpp::wrap(prior)), Shield<SEXP>(Rcpp::wrap(VB)), Shield<SEXP>(Rcpp::wrap(starting_values)), Shield<SEXP>(Rcpp::wrap(thin)), Shield<SEXP>(Rcpp::wrap(sample_s_)), Shield<SEXP>(Rcpp::wrap(show_progress)));
+            rcpp_result_gen = p_bsvar_sv_cpp(Shield<SEXP>(Rcpp::wrap(S)), Shield<SEXP>(Rcpp::wrap(Y)), Shield<SEXP>(Rcpp::wrap(X)), Shield<SEXP>(Rcpp::wrap(prior)), Shield<SEXP>(Rcpp::wrap(VB)), Shield<SEXP>(Rcpp::wrap(starting_values)), Shield<SEXP>(Rcpp::wrap(thin)), Shield<SEXP>(Rcpp::wrap(centred_sv)), Shield<SEXP>(Rcpp::wrap(show_progress)));
         }
         if (rcpp_result_gen.inherits("interrupted-error"))
             throw Rcpp::internal::InterruptedException();
@@ -256,17 +256,17 @@ namespace bsvars {
         return Rcpp::as<Rcpp::List >(rcpp_result_gen);
     }
 
-    inline Rcpp::List forecast_bsvar_sv(arma::cube& posterior_B, arma::cube& posterior_A, arma::vec& posterior_h_T, arma::mat& posterior_rho, arma::mat& posterior_omega, arma::vec& X_T, const int horizon) {
-        typedef SEXP(*Ptr_forecast_bsvar_sv)(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP);
+    inline Rcpp::List forecast_bsvar_sv(arma::cube& posterior_B, arma::cube& posterior_A, arma::vec& posterior_h_T, arma::mat& posterior_rho, arma::mat& posterior_omega, arma::vec& X_T, const int horizon, const bool centred_sv = FALSE) {
+        typedef SEXP(*Ptr_forecast_bsvar_sv)(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP);
         static Ptr_forecast_bsvar_sv p_forecast_bsvar_sv = NULL;
         if (p_forecast_bsvar_sv == NULL) {
-            validateSignature("Rcpp::List(*forecast_bsvar_sv)(arma::cube&,arma::cube&,arma::vec&,arma::mat&,arma::mat&,arma::vec&,const int)");
+            validateSignature("Rcpp::List(*forecast_bsvar_sv)(arma::cube&,arma::cube&,arma::vec&,arma::mat&,arma::mat&,arma::vec&,const int,const bool)");
             p_forecast_bsvar_sv = (Ptr_forecast_bsvar_sv)R_GetCCallable("bsvars", "_bsvars_forecast_bsvar_sv");
         }
         RObject rcpp_result_gen;
         {
             RNGScope RCPP_rngScope_gen;
-            rcpp_result_gen = p_forecast_bsvar_sv(Shield<SEXP>(Rcpp::wrap(posterior_B)), Shield<SEXP>(Rcpp::wrap(posterior_A)), Shield<SEXP>(Rcpp::wrap(posterior_h_T)), Shield<SEXP>(Rcpp::wrap(posterior_rho)), Shield<SEXP>(Rcpp::wrap(posterior_omega)), Shield<SEXP>(Rcpp::wrap(X_T)), Shield<SEXP>(Rcpp::wrap(horizon)));
+            rcpp_result_gen = p_forecast_bsvar_sv(Shield<SEXP>(Rcpp::wrap(posterior_B)), Shield<SEXP>(Rcpp::wrap(posterior_A)), Shield<SEXP>(Rcpp::wrap(posterior_h_T)), Shield<SEXP>(Rcpp::wrap(posterior_rho)), Shield<SEXP>(Rcpp::wrap(posterior_omega)), Shield<SEXP>(Rcpp::wrap(X_T)), Shield<SEXP>(Rcpp::wrap(horizon)), Shield<SEXP>(Rcpp::wrap(centred_sv)));
         }
         if (rcpp_result_gen.inherits("interrupted-error"))
             throw Rcpp::internal::InterruptedException();
@@ -793,17 +793,38 @@ namespace bsvars {
         return Rcpp::as<arma::vec >(rcpp_result_gen);
     }
 
-    inline Rcpp::List svar_nc1(arma::rowvec& aux_h_n, double& aux_rho_n, double& aux_omega_n, double& aux_sigma2_omega_n, double& aux_s_n, arma::urowvec& aux_S_n, const arma::rowvec& u, const Rcpp::List& prior, bool sample_s_ = true) {
-        typedef SEXP(*Ptr_svar_nc1)(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP);
+    inline Rcpp::List svar_nc1(arma::rowvec& aux_h_n, double& aux_rho_n, double& aux_omega_n, double& aux_sigma2v_n, double& aux_sigma2_omega_n, double& aux_s_n, arma::urowvec& aux_S_n, const arma::rowvec& u, const Rcpp::List& prior, bool sample_s_ = true) {
+        typedef SEXP(*Ptr_svar_nc1)(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP);
         static Ptr_svar_nc1 p_svar_nc1 = NULL;
         if (p_svar_nc1 == NULL) {
-            validateSignature("Rcpp::List(*svar_nc1)(arma::rowvec&,double&,double&,double&,double&,arma::urowvec&,const arma::rowvec&,const Rcpp::List&,bool)");
+            validateSignature("Rcpp::List(*svar_nc1)(arma::rowvec&,double&,double&,double&,double&,double&,arma::urowvec&,const arma::rowvec&,const Rcpp::List&,bool)");
             p_svar_nc1 = (Ptr_svar_nc1)R_GetCCallable("bsvars", "_bsvars_svar_nc1");
         }
         RObject rcpp_result_gen;
         {
             RNGScope RCPP_rngScope_gen;
-            rcpp_result_gen = p_svar_nc1(Shield<SEXP>(Rcpp::wrap(aux_h_n)), Shield<SEXP>(Rcpp::wrap(aux_rho_n)), Shield<SEXP>(Rcpp::wrap(aux_omega_n)), Shield<SEXP>(Rcpp::wrap(aux_sigma2_omega_n)), Shield<SEXP>(Rcpp::wrap(aux_s_n)), Shield<SEXP>(Rcpp::wrap(aux_S_n)), Shield<SEXP>(Rcpp::wrap(u)), Shield<SEXP>(Rcpp::wrap(prior)), Shield<SEXP>(Rcpp::wrap(sample_s_)));
+            rcpp_result_gen = p_svar_nc1(Shield<SEXP>(Rcpp::wrap(aux_h_n)), Shield<SEXP>(Rcpp::wrap(aux_rho_n)), Shield<SEXP>(Rcpp::wrap(aux_omega_n)), Shield<SEXP>(Rcpp::wrap(aux_sigma2v_n)), Shield<SEXP>(Rcpp::wrap(aux_sigma2_omega_n)), Shield<SEXP>(Rcpp::wrap(aux_s_n)), Shield<SEXP>(Rcpp::wrap(aux_S_n)), Shield<SEXP>(Rcpp::wrap(u)), Shield<SEXP>(Rcpp::wrap(prior)), Shield<SEXP>(Rcpp::wrap(sample_s_)));
+        }
+        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);
+    }
+
+    inline Rcpp::List svar_ce1(arma::rowvec& aux_h_n, double& aux_rho_n, double& aux_omega_n, double& aux_sigma2v_n, double& aux_sigma2_omega_n, double& aux_s_n, arma::urowvec& aux_S_n, const arma::rowvec& u, const Rcpp::List& prior, bool sample_s_ = true) {
+        typedef SEXP(*Ptr_svar_ce1)(SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP,SEXP);
+        static Ptr_svar_ce1 p_svar_ce1 = NULL;
+        if (p_svar_ce1 == NULL) {
+            validateSignature("Rcpp::List(*svar_ce1)(arma::rowvec&,double&,double&,double&,double&,double&,arma::urowvec&,const arma::rowvec&,const Rcpp::List&,bool)");
+            p_svar_ce1 = (Ptr_svar_ce1)R_GetCCallable("bsvars", "_bsvars_svar_ce1");
+        }
+        RObject rcpp_result_gen;
+        {
+            RNGScope RCPP_rngScope_gen;
+            rcpp_result_gen = p_svar_ce1(Shield<SEXP>(Rcpp::wrap(aux_h_n)), Shield<SEXP>(Rcpp::wrap(aux_rho_n)), Shield<SEXP>(Rcpp::wrap(aux_omega_n)), Shield<SEXP>(Rcpp::wrap(aux_sigma2v_n)), Shield<SEXP>(Rcpp::wrap(aux_sigma2_omega_n)), Shield<SEXP>(Rcpp::wrap(aux_s_n)), Shield<SEXP>(Rcpp::wrap(aux_S_n)), Shield<SEXP>(Rcpp::wrap(u)), Shield<SEXP>(Rcpp::wrap(prior)), Shield<SEXP>(Rcpp::wrap(sample_s_)));
         }
         if (rcpp_result_gen.inherits("interrupted-error"))
             throw Rcpp::internal::InterruptedException();

diff --git a/inst/tinytest/test_estimate_bsvar_sv.R b/inst/tinytest/test_estimate_bsvar_sv.R
@@ -35,3 +35,22 @@ expect_identical(
   run_no3$last_draw$starting_values$B[1,1],
   info = "estimate_bsvar_sv: the last_draw(s) of a normal and pipe run to be identical."
 )
+
+
+set.seed(1)
+suppressMessages(
+  specification_no1 <- specify_bsvar_sv$new(us_fiscal_lsuw, centred_sv = TRUE)
+)
+run_no1             <- estimate(specification_no1, 3, 1, show_progress = FALSE)
+
+set.seed(1)
+suppressMessages(
+  specification_no2 <- specify_bsvar_sv$new(us_fiscal_lsuw, centred_sv = TRUE)
+)
+run_no2             <- estimate(specification_no2, 3, 1, show_progress = FALSE)
+
+expect_identical(
+  run_no1$last_draw$starting_values$B[1,1],
+  run_no2$last_draw$starting_values$B[1,1],
+  info = "estimate_bsvar_sv centred: the last_draw(s) of two runs to be identical."
+)