changed handling of nonstationary variables:

- oo_.mean, oo_.var, oo_.gamma_y contains all selected variables
- moments of non-stationary variables are set to NaN
options_.Schur_vec_tolerance lowered to 10^-11


git-svn-id: https://www.dynare.org/svn/dynare/trunk@2810 ac1d8469-bf42-47a9-8791-bf33cf982152

matlab/disp_th_moments.m

@@ -35,8 +35,10 @@ function disp_th_moments(dr,var_list)
     end
   end
   
-  [oo_.gamma_y,ivar] = th_autocovariances(dr,ivar,M_,options_);
+  [oo_.gamma_y,stationary_vars] = th_autocovariances(dr,ivar,M_,options_);
   m = dr.ys(ivar);
+  non_stationary_vars = setdiff(1:length(ivar),stationary_vars);
+  m(ivar(non_stationary_vars)) = NaN;
 
   
   i1 = find(abs(diag(oo_.gamma_y{1})) > 1e-12);
@@ -50,14 +52,15 @@ function disp_th_moments(dr,var_list)
   oo_.mean = m;
   oo_.var = oo_.gamma_y{1};
   
-  lh = size(deblank(M_.endo_names(ivar,:)),2)+2;
   if options_.nomoments == 0
     title='THEORETICAL MOMENTS';
     if options_.hp_filter
       title = [title ' (HP filter, lambda = ' int2str(options_.hp_filter) ')'];
     end
     headers=strvcat('VARIABLE','MEAN','STD. DEV.','VARIANCE');
-    dyntable(title,headers,deblank(M_.endo_names(ivar,:)),z,lh,11,4);
+    labels = deblank(M_.endo_names(ivar,:));
+    lh = size(labels,2)+2;
+    dyntable(title,headers,labels,z,lh,11,4);
     if M_.exo_nbr > 1
       disp(' ')
       title='VARIANCE DECOMPOSITION (in percent)';
@@ -68,8 +71,9 @@ function disp_th_moments(dr,var_list)
       headers = M_.exo_names;
       headers(M_.exo_names_orig_ord,:) = headers;
       headers = strvcat(' ',headers);
-      dyntable(title,headers,deblank(M_.endo_names(ivar(i1),:)),100*oo_.gamma_y{options_.ar+2}(i1,:), ...
-	    lh,8,2);
+      lh = size(deblank(M_.endo_names(ivar(stationary_vars),:)),2)+2;
+      dyntable(title,headers,deblank(M_.endo_names(ivar(stationary_vars), ...
+                                                   :)),100*oo_.gamma_y{options_.ar+2}(stationary_vars,:),lh,8,2);
     end
   end
   
@@ -79,10 +83,11 @@ function disp_th_moments(dr,var_list)
     if options_.hp_filter
       title = [title ' (HP filter, lambda = ' int2str(options_.hp_filter) ')'];
     end
-    labels = deblank(M_.endo_names(ivar,:));
-    headers = strvcat('Variables',labels(i1,:));
+    labels = deblank(M_.endo_names(ivar(i1),:));
+    headers = strvcat('Variables',labels);
     corr = oo_.gamma_y{1}(i1,i1)./(sd(i1)*sd(i1)');
-    dyntable(title,headers,labels(i1,:),corr,lh,8,4);
+    lh = size(labels,2)+2;
+    dyntable(title,headers,labels,corr,lh,8,4);
   end
   
   if options_.ar > 0
@@ -98,7 +103,8 @@ function disp_th_moments(dr,var_list)
       oo_.autocorr{i} = oo_.gamma_y{i+1};
       z(:,i) = diag(oo_.gamma_y{i+1}(i1,i1));
     end
-    dyntable(title,headers,labels,z,0,8,4);
+    lh = size(labels,2)+2;
+    dyntable(title,headers,labels,z,lh,8,4);
   end
   
 % 10/09/02 MJ 

matlab/global_initialization.m

@@ -41,9 +41,9 @@ function global_initialization()
   options_.gstep = 1e-2;
   options_.debug = 0;
   options_.initval_file = 0;
-  options_.Schur_vec_tol = 1e-8; % used to find nonstationary variables
-                                 % in Schur decomposition of the
-                                 % transition matrix
+  options_.Schur_vec_tol = 1e-11; % used to find nonstationary variables
+                                  % in Schur decomposition of the
+                                  % transition matrix
   options_.qz_criterium = 1.000001;
   options_.lyapunov_complex_threshold = 1e-15;
   options_.solve_tolf = eps^(1/3);

matlab/th_autocovariances.m

@@ -1,4 +1,4 @@
-function [Gamma_y,ivar] = th_autocovariances(dr,ivar,M_,options_,nodecomposition)
+function [Gamma_y,stationary_vars] = th_autocovariances(dr,ivar,M_,options_,nodecomposition)
 % Computes the theoretical auto-covariances, Gamma_y, for an AR(p) process 
 % with coefficients dr.ghx and dr.ghu and shock variances Sigma_e_
 % for a subset of variables ivar (indices in lgy_)
@@ -19,7 +19,8 @@ function [Gamma_y,ivar] = th_autocovariances(dr,ivar,M_,options_,nodecomposition
 %                                      Gamma_y{nar+2}   [double]  Variance decomposition.  
 %                                      Gamma_y{nar+3}   [double]  Expectation of the endogenous variables associated with a second 
 %                                                                 order approximation.    
-%   ivar              [integer]      Vector of indices for a subset of variables.
+%   stationary_vars   [integer]      Vector of indices of stationary
+%                                           variables in declaration order
 %
 % SPECIAL REQUIREMENTS
 %   
@@ -45,6 +46,7 @@ function [Gamma_y,ivar] = th_autocovariances(dr,ivar,M_,options_,nodecomposition
         nodecomposition = 0;
     end
     
+    endo_nbr = M_.endo_nbr;
     exo_names_orig_ord  = M_.exo_names_orig_ord;
     if exist('OCTAVE_VERSION')
         warning('off', 'Octave:divide-by-zero')
@@ -54,7 +56,7 @@ function [Gamma_y,ivar] = th_autocovariances(dr,ivar,M_,options_,nodecomposition
     nar = options_.ar;
     Gamma_y = cell(nar+1,1);
     if isempty(ivar)
-        ivar = [1:M_.endo_nbr]';
+        ivar = [1:endo_nbr]';
     end
     nvar = size(ivar,1);
   
@@ -63,8 +65,8 @@ function [Gamma_y,ivar] = th_autocovariances(dr,ivar,M_,options_,nodecomposition
     npred = dr.npred;
     nstatic = dr.nstatic;
     kstate = dr.kstate;
-    order = dr.order_var;
-    iv(order) = [1:length(order)];
+    order_var = dr.order_var;
+    inv_order_var = dr.inv_order_var;
     nx = size(ghx,2);
   
     ikx = [nstatic+1:nstatic+npred];
@@ -94,10 +96,12 @@ function [Gamma_y,ivar] = th_autocovariances(dr,ivar,M_,options_,nodecomposition
     [A,B] = kalman_transition_matrix(dr,ipred,1:nx,dr.transition_auxiliary_variables,M_.exo_nbr);
     if options_.order == 2 | options_.hp_filter == 0
         [vx, u] =  lyapunov_symm(A,B*M_.Sigma_e*B',options_.qz_criterium,options_.lyapunov_complex_threshold);
-        iky = iv(ivar);
+        iky = inv_order_var(ivar);
+        stationary_vars = (1:length(ivar))';
         if ~isempty(u)
-            iky = iky(find(all(abs(ghx(iky,:)*u) < options_.Schur_vec_tol,2)));
-            ivar = dr.order_var(iky);
+            x = abs(ghx*u);
+            iky = iky(find(all(x(iky,:) < options_.Schur_vec_tol,2)));
+            stationary_vars = find(all(x(inv_order_var(ivar(stationary_vars)),:) < options_.Schur_vec_tol,2));
         end
         aa = ghx(iky,:);
         bb = ghu(iky,:);
@@ -109,23 +113,28 @@ function [Gamma_y,ivar] = th_autocovariances(dr,ivar,M_,options_,nodecomposition
         end
     end
     if options_.hp_filter == 0
-        Gamma_y{1} = aa*vx*aa'+ bb*M_.Sigma_e*bb';
-        k = find(abs(Gamma_y{1}) < 1e-12);
-        Gamma_y{1}(k) = 0;    
+        v = NaN*ones(nvar,nvar);
+        v(stationary_vars,stationary_vars) = aa*vx*aa'+ bb*M_.Sigma_e*bb';
+        k = find(abs(v) < 1e-12);
+        v(k) = 0;
+        Gamma_y{1} = v;
         % autocorrelations
         if nar > 0
             vxy = (A*vx*aa'+ghu1*M_.Sigma_e*bb');    
             sy = sqrt(diag(Gamma_y{1}));
+            sy = sy(stationary_vars);
             sy = sy *sy';
-            Gamma_y{2} = aa*vxy./sy;
+            Gamma_y{2} = NaN*ones(nvar,nvar);
+            Gamma_y{2}(stationary_vars,stationary_vars) = aa*vxy./sy;
             for i=2:nar
                 vxy = A*vxy;
-                Gamma_y{i+1} = aa*vxy./sy;
+                Gamma_y{i+1} = NaN*ones(nvar,nvar);
+                Gamma_y{i+1}(stationary_vars,stationary_vars) = aa*vxy./sy;
             end
         end
         % variance decomposition
         if ~nodecomposition && M_.exo_nbr > 1
-            Gamma_y{nar+2} = zeros(length(ivar),M_.exo_nbr);
+            Gamma_y{nar+2} = zeros(nvar,M_.exo_nbr);
             SS(exo_names_orig_ord,exo_names_orig_ord)=M_.Sigma_e+1e-14*eye(M_.exo_nbr);
             cs = chol(SS)';
             b1(:,exo_names_orig_ord) = ghu1;
@@ -136,12 +145,12 @@ function [Gamma_y,ivar] = th_autocovariances(dr,ivar,M_,options_,nodecomposition
             vv = diag(aa*vx*aa'+b2*b2');
             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}(:,i) = abs(diag(aa*vx1*aa'+b2(:,i)*b2(:,i)'))./vv;
+                Gamma_y{nar+2}(stationary_vars,i) = abs(diag(aa*vx1*aa'+b2(:,i)*b2(:,i)'))./vv;
             end
         end
     else% ==> Theoretical HP filter.
         if options_.order < 2
-            iky = iv(ivar);  
+            iky = inv_order_var(ivar);  
             aa = ghx(iky,:);
             bb = ghu(iky,:);
         end