From 1224cec429a532be64c843051f2ce3f0db2d6897 Mon Sep 17 00:00:00 2001
From: adjemian <adjemian@ac1d8469-bf42-47a9-8791-bf33cf982152>
Date: Fri, 21 Sep 2007 12:17:56 +0000
Subject: [PATCH] Explosive VARs (the companion matrix has eigenvalues outside
 the unit circle) are discarded from the posterior distribution (by default).

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1388 ac1d8469-bf42-47a9-8791-bf33cf982152
---
 matlab/bvar_forecast.m | 29 ++++++++++++++++++++++++++---
 1 file changed, 26 insertions(+), 3 deletions(-)

diff --git a/matlab/bvar_forecast.m b/matlab/bvar_forecast.m
index 6c52e5cb8..65a535c50 100644
--- a/matlab/bvar_forecast.m
+++ b/matlab/bvar_forecast.m
@@ -2,12 +2,14 @@ function bvar_forecast(nlags)
 
     global options_
     
+    nvar = size(options_.varobs,1);
+    
     options_ = set_default_option(options_, 'bvar_replic', 2000);
     if options_.forecast == 0
         error('bvar_forecast: you must specify "forecast" option')
     end
-
     [ny, nx, posterior, prior, forecast_data] = bvar_toolbox(nlags);
+
         
     sims_no_shock = NaN(options_.forecast, ny, options_.bvar_replic);
     sims_with_shocks = NaN(options_.forecast, ny, options_.bvar_replic);
@@ -20,7 +22,14 @@ function bvar_forecast(nlags)
     
     k = ny*nlags+nx;
     
-    for d = 1:options_.bvar_replic
+    % Declaration of the companion matrix:
+    Companion_matrix = diag(ones(nvar*(nlags-1),1),-nvar);
+    
+    p = 0;
+    d = 0;
+    while d<=options_.bvar_replic  
+        d = d+1;
+        
         Sigma = rand_inverse_wishart(ny, posterior.df, S_inv_upper_chol);
 
         % Option 'lower' of chol() not available in old versions of
@@ -29,6 +38,15 @@ function bvar_forecast(nlags)
 
         Phi = rand_matrix_normal(k, ny, posterior.PhiHat, XXi_lower_chol, Sigma_lower_chol);
         
+        % All the eigenvalues of the companion matrix have to be on or inside the unit circle
+        Companion_matrix(1:nvar,:) = Phi(1:nvar*nlags,:)'; 
+        test = (abs(eig(Companion_matrix)));
+        if any(test>1.0000000000001)
+            p = p+1;
+            d = d-1;
+            continue
+        end
+        
         % Without shocks
         lags_data = forecast_data.initval;
         for t = 1:options_.forecast
@@ -48,7 +66,12 @@ function bvar_forecast(nlags)
             sims_with_shocks(t, :, d) = y;
         end
     end
-
+    
+    disp('')
+    disp(['Some of the VAR models sampled from the posterior distribution'])
+    disp(['were found to be explosive (' int2str(p) ').'])
+    disp('')
+    
     % Plot graphs
     sims_no_shock_mean = mean(sims_no_shock, 3);