zero logvars for homoscedastic factors

luisgruber · Feb 2, 2024 · 86895cd · 86895cd
1 parent 19573f6
commit 86895cd
Show file tree

Hide file tree

Showing 3 changed files with 21 additions and 3 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -1,5 +1,7 @@
 # bayesianVARs (development version)
 
+* bugfix concerning VAR with factor structure on errors with homoscedastic factors.
+
 # bayesianVARs 0.1.2
 
 * Added minimum version to factorstochvol in the Imports field of the DESCRIPTION file in order to avoid unnecessary building errors. Thanks to Sergey Fedorov for pointing this out.

diff --git a/R/bvar_wrapper.R b/R/bvar_wrapper.R
@@ -551,7 +551,7 @@ bvar <- function(data,
                       phi = c(rep(.8, M), rep(.8, prior_sigma$factor_factors)) + pmin(rnorm(M + prior_sigma$factor_factors, sd=.06), .095),
                       sigma = rep(.1, M + prior_sigma$factor_factors) + rgamma(M + prior_sigma$factor_factors, 1, 10))
     startlogvar <- matrix(startpara["mu",][1] + rnorm(Tobs*(M + prior_sigma$factor_factors)), Tobs, M + prior_sigma$factor_factors)
-    startlogvar[,M+which(isFALSE(prior_sigma$sv_heteroscedastic[-c(1:M)]))] <- 0 # !!! important, needed for factorstochvol: if factor is assumed to be homoscedastic, the corresponding column in logvar has to be 0!!!
+    startlogvar[,M+which(prior_sigma$sv_heteroscedastic[-c(1:M)]==FALSE)] <- 0 # !!! important, needed for factorstochvol: if factor is assumed to be homoscedastic, the corresponding column in logvar has to be 0!!!
     startlogvar0 <- startpara["mu",][1] + rnorm(M + prior_sigma$factor_factors)
     starttau2 <- if(!prior_sigma$factor_ngprior){ # if prior is 'normal'
       prior_sigma$factor_starttau2

diff --git a/tests/testthat/test-bvar.R b/tests/testthat/test-bvar.R
@@ -104,7 +104,7 @@ test_that("flat prior factor", {
                                factor_priorfacloadtype = "normal",
                                factor_priorfacload = 1e6,
                                factor_heteroskedastic = c(FALSE, FALSE),
-                               factor_priorhomoskedastic = matrix(c(0.01,.01), ncol(data), 2))
+                               factor_priorhomoskedastic = matrix(c(1e-6,1e-6), ncol(data), 2))
  set.seed(123)
   mod <- bvar(data, lags = 2, prior_intercept = 1e6, prior_phi = phi,
               prior_sigma = sigma, draws = 10000)
@@ -122,7 +122,7 @@ test_that("flat prior cholesky", {
                                cholesky_U_prior = "normal",
                                cholesky_normal_sds = 1e6,
                                cholesky_heteroscedastic = FALSE,
-                               cholesky_priorhomoscedastic = matrix(c(0.01,0.01), ncol(data), 2))
+                               cholesky_priorhomoscedastic = matrix(c(1e-6,1e-6), ncol(data), 2))
   set.seed(123)
   mod <- bvar(data, lags = 2, prior_intercept = 1e6, prior_phi = phi,
               prior_sigma = sigma, draws = 10000)
@@ -131,3 +131,19 @@ test_that("flat prior cholesky", {
 
   expect_lt(max(abs(ols-phi_post_mean)),0.01)
 })
+
+test_that("homoscedastic factor VAR input", {
+  data <- usmacro_growth[,c("GDPC1","CPIAUCSL","FEDFUNDS")]
+  phi <- specify_prior_phi(data = data, lags = 2L, prior = "normal",
+                           normal_sds = 1e6)
+  sigma <- specify_prior_sigma(data = data, type = "factor",
+                               factor_factors = 2,
+                               factor_priorfacloadtype = "normal",
+                               factor_priorfacload = 1e6,
+                               factor_heteroskedastic = c(FALSE, FALSE),
+                               factor_priorhomoskedastic = matrix(c(1e-6,1e-6), ncol(data), 2))
+
+  mod <- bvar(data, lags = 2, prior_intercept = 1e6, prior_phi = phi,
+              prior_sigma = sigma, draws = 10, burnin = 10)
+  expect_true(all(mod$logvar[,-seq_len(ncol(mod$Y)),]==0))
+})