From f76a657ceac93b488e0b242f9d30b20c0a9d5856 Mon Sep 17 00:00:00 2001
From: Michel Juillard <michel.juillard@ens.fr>
Date: Sat, 3 Apr 2010 11:27:49 +0200
Subject: [PATCH] computes now variance decomposition relative to the sum of
 the effects of individual shocks rather than aggregate variance. When the
 aggregate variance differs from the shock of the sum of the effects of
 individual shocks by more than 0.01% a warning is displayed.

This behaviour is documented in the reference manual.
(cherry picked from commit c2f7f0a555882acd50cbb93b4c7c79603fac58ba)
---
 doc/manual.xml              |  3 +++
 matlab/th_autocovariances.m | 13 ++++++++++++-
 2 files changed, 15 insertions(+), 1 deletion(-)

diff --git a/doc/manual.xml b/doc/manual.xml
index 02dc94fa33..d49a134e98 100644
--- a/doc/manual.xml
+++ b/doc/manual.xml
@@ -2014,6 +2014,9 @@ steady(homotopy_mode = 1, homotopy_steps = 50);
 <para>
   Variance decomposition, correlation, autocorrelation are only displayed for variables with positive variance. Impulse response functions are only plotted for variables with response larger than 10<superscript>-10</superscript>.
 </para>
+<para>
+  Variance decomposition is computed relative to the sum of the contribution of each shock. Normally, this is of course equal to aggregate variance, but if a model generates very large variances, it may happen that, due to numerical error, the two differ by a significant amount. Dynare issues a warning if the maximum relative difference between the sum of the contribution of each shock and aggregate variance is larger than 0.01 %.
+</para>
 <para>
   Currently, the IRFs are only plotted for 12 variables. Select the ones you want to see, if your model contains more than 12 endogenous variables.
 </para>
diff --git a/matlab/th_autocovariances.m b/matlab/th_autocovariances.m
index 8262e3948c..4da9ca22a8 100644
--- a/matlab/th_autocovariances.m
+++ b/matlab/th_autocovariances.m
@@ -147,9 +147,20 @@ if options_.hp_filter == 0
         b2 = b2*cs;
         vx  = lyapunov_symm(A,b1*b1',options_.qz_criterium,options_.lyapunov_complex_threshold,1);
         vv = diag(aa*vx*aa'+b2*b2');
+        vv2 = 0;
         for i=1:M_.exo_nbr
             vx1 = lyapunov_symm(A,b1(:,i)*b1(:,i)',options_.qz_criterium,options_.lyapunov_complex_threshold,2);
-            Gamma_y{nar+2}(stationary_vars,i) = abs(diag(aa*vx1*aa'+b2(:,i)*b2(:,i)'))./vv;
+            vx2 = abs(diag(aa*vx1*aa'+b2(:,i)*b2(:,i)'));
+            Gamma_y{nar+2}(stationary_vars,i) = vx2;
+            vv2 = vv2 +vx2;
+        end
+        if max(abs(vv2-vv)./vv) > 1e-4
+            warning(['Aggregate variance and sum of variances by shocks ' ...
+                     'differ by more than 0.01 %'])
+        end
+        for i=1:M_.exo_nbr
+            Gamma_y{nar+2}(stationary_vars,i) = Gamma_y{nar+ ...
+                                2}(stationary_vars,i)./vv2;
         end
     end
 else% ==> Theoretical HP filter.
-- 
GitLab