v4.1: Changes related to (posterior) 2nd order moments.

+ Efficiency changes:
++ Do not compute the variance decomposition if not needed.
++ Do not apply the Schur decomposition for each shock in case of a variance decomposition.
+ Bug corrections and cosmetic changes.
+ Corrections of some headers.


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

matlab/compute_moments_varendo.m

@@ -51,7 +51,7 @@ function oo_ = compute_moments_varendo(options_,M_,oo_,var_list_)
     end
     % VARIANCE DECOMPOSITION.
     for i=1:NumberOfEndogenousVariables
-        for j=i:NumberOfExogenousVariables
+        for j=1:NumberOfExogenousVariables
             oo_ = posterior_analysis('decomposition',var_list_(i,:),M_.exo_names(j,:),[],options_,M_,oo_);
         end
     end
\ No newline at end of file

matlab/dsge_posterior_theoretical_correlation.m

@@ -20,7 +20,7 @@ function [nvar,vartan,CorrFileNumber] = dsge_posterior_theoretical_correlation(S
 %                                                   th_autocovariances.m    
 %                                                   posterior_moments.m
 
-% Copyright (C) 2007-2008 Dynare Team
+% Copyright (C) 2007-2009 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -38,6 +38,7 @@ function [nvar,vartan,CorrFileNumber] = dsge_posterior_theoretical_correlation(S
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
 type = 'posterior';
+nodecomposition = 1;
     
 % Set varlist (vartan)
 [ivar,vartan] = set_stationary_variables_list;
@@ -80,7 +81,7 @@ for file = 1:NumberOfDrawsFiles
             set_parameters(pdraws{linee,1});
             [dr,info] = resol(oo_.steady_state,0);
         end
-        tmp = th_autocovariances(dr,ivar,M_,options_);
+        tmp = th_autocovariances(dr,ivar,M_,options_,nodecomposition);
         for i=1:nar
             Correlation_array(linea,:,:,i) = tmp{i+1};
         end

matlab/dsge_posterior_theoretical_covariance.m

@@ -19,7 +19,7 @@ function [nvar,vartan,CovarFileNumber] = dsge_posterior_theoretical_covariance(S
 %                                                   th_autocovariances.m    
 %                                                   posterior_moments.m
 
-% Copyright (C) 2007-2008 Dynare Team
+% Copyright (C) 2007-2009 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -37,7 +37,8 @@ function [nvar,vartan,CovarFileNumber] = dsge_posterior_theoretical_covariance(S
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
 type = 'posterior';
-    
+nodecomposition = 1;
+
 % Set varlist (vartan)
 [ivar,vartan] = set_stationary_variables_list;
 nvar = length(ivar);
@@ -79,7 +80,7 @@ for file = 1:NumberOfDrawsFiles
             set_parameters(pdraws{linee,1});
             [dr,info] = resol(oo_.steady_state,0);
         end
-        tmp = th_autocovariances(dr,ivar,M_,options_);
+        tmp = th_autocovariances(dr,ivar,M_,options_,nodecomposition);
         for i=1:nvar
             for j=i:nvar
                 Covariance_matrix(linea,symmetric_matrix_index(i,j,nvar)) = tmp{1}(i,j);
@@ -94,7 +95,7 @@ for file = 1:NumberOfDrawsFiles
                 Covariance_matrix = zeros(NumberOfLinesInTheLastCovarFile,nvar*(nvar+1)/2);
                 NumberOfCovarLines = NumberOfLinesInTheLastCovarFile;
                 CovarFileNumber = CovarFileNumber - 1;
-            elseif test<0;
+            elseif test<0
                 Covariance_matrix = zeros(MaXNumberOfCovarLines,nvar*(nvar+1)/2);
             else
                 clear('Covariance_matrix');

matlab/dsge_posterior_theoretical_variance_decomposition.m

@@ -18,7 +18,7 @@ function [nvar,vartan,NumberOfDecompFiles] = ...
 %                                                   th_autocovariances.m    
 %                                                   posterior_moments.m
 
-% Copyright (C) 2007-2008 Dynare Team
+% Copyright (C) 2007-2009 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -36,6 +36,7 @@ function [nvar,vartan,NumberOfDecompFiles] = ...
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
 type = 'posterior';% To be defined as a input argument later...
+nodecomposition = 0;
 
 % Set varlist (vartan)
 [ivar,vartan] = set_stationary_variables_list;
@@ -51,7 +52,6 @@ NumberOfDrawsFiles = length(DrawsFiles);
 
 nexo = M_.exo_nbr;
 
-
 NumberOfDrawsFiles = rows(DrawsFiles);
 NumberOfSavedElementsPerSimulation = nvar*(nexo+1);
 MaXNumberOfDecompLines = ceil(options_.MaxNumberOfBytes/NumberOfSavedElementsPerSimulation/8);
@@ -84,7 +84,7 @@ for file = 1:NumberOfDrawsFiles
             set_parameters(pdraws{linee,1});
             [dr,info] = resol(oo_.steady_state,0);
         end
-        tmp = th_autocovariances(dr,ivar,M_,options_);
+        tmp = th_autocovariances(dr,ivar,M_,options_,nodecomposition);
         for i=1:nvar
             for j=1:nexo
                 Decomposition_array(linea,(i-1)*nexo+j) = tmp{2}(i,j);

matlab/lyapunov_symm.m

-function [x,u]=lyapunov_symm(a,b,qz_criterium)
-% Solves the Lyapunov equation x-a*x*a' = b, for b (and then x) symmetrical
-% if a has some unit roots, the function computes only the solution of the stable subsystem
+function [x,u] = lyapunov_symm(a,b,qz_criterium,method)
+% Solves the Lyapunov equation x-a*x*a' = b, for b and x symmetric matrices.
+% If a has some unit roots, the function computes only the solution of the stable subsystem.
 %  
 % INPUTS:
-%   a:      coefficient matrix (n x n)  
-%   b:      coefficient square matrix (n x n)
-%
+%   a        [double]    n*n matrix.  
+%   b        [double]    n*n matrix
+%   method   [integer]   Scalar, if method=0 [default] then U, T, n and k are not persistent.  
+%                                   method=1 then U, T, n and k are declared as persistent variables and the schur decomposition is triggered.    
+%                                   method=2 then U, T, n and k are declared as persistent variables and the schur decomposition is not performed.
 % OUTPUTS
-%   x:      solution matrix (m x m)
-%   ns_var: vector of indices of non-stationary variables (p x 1)
-%           (m + p = n)
-%   u:      Schur vectors associated with unit roots  
+%   x:      [double]    m*m solution matrix of the lyapunov equation, where m is the dimension of the stable subsystem.
+%   u:      [double]    Schur vectors associated with unit roots  
 %
 % ALGORITHM
-%   uses reordered Schur decomposition
+%   Uses reordered Schur decomposition
 %
 % SPECIAL REQUIREMENTS
 %   needs Matlab >= 7.0.1 for ordeig function (otherwise uses my_ordeig)
 
-% Copyright (C) 2006-2008 Dynare Team
+% Copyright (C) 2006-2009 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -34,76 +34,90 @@ function [x,u]=lyapunov_symm(a,b,qz_criterium)
 %
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
-
-  ns_var = [];
-  u = [];
-  
-  info = 0;
-  if size(a,1) == 1
-    x=b/(1-a*a);
-    return
-  end
-  [u,t] = schur(a);
-  if exist('OCTAVE_VERSION') || matlab_ver_less_than('7.0.1')
-    e1 = abs(my_ordeig(t)) > 2-qz_criterium;
-  else
-    e1 = abs(ordeig(t)) > 2-qz_criterium;
-  end
-  k = sum(e1);
-  if exist('ordschur','builtin')
-    % selects stable roots
-    [u,t] = ordschur(u,t,e1); 
-    n = length(e1)-k;
-    b=u(:,k+1:end)'*b*u(:,k+1:end);
-    t = t(k+1:end,k+1:end);
-  elseif k > 0
-    % problem for Matlab version that don't have ordschur
-    error(['lyapunov_sym: you need a Matlab version > 6.5 to handle models' ...
-	   ' with unit roots'])
-  else
-    % no unit root
-    n = length(e1);
-    b=u'*b*u;
-  end
-  x=zeros(n,n);
-  i = n;
-  while i >= 2
-    if t(i,i-1) == 0
-      if i == n
-	c = zeros(n,1);
-      else
-	c = t(1:i,:)*(x(:,i+1:end)*t(i,i+1:end)')+...
-	    t(i,i)*t(1:i,i+1:end)*x(i+1:end,i);
-      end
-      q = eye(i)-t(1:i,1:i)*t(i,i);
-      x(1:i,i) = q\(b(1:i,i)+c);
-      x(i,1:i-1) = x(1:i-1,i)';
-      i = i - 1;
+    if nargin<4
+        method = 0;
+    end
+    if method
+        persistent U T k n
     else
-      if i == n
-	c = zeros(n,1);
-	c1 = zeros(n,1);
-      else
-	c = t(1:i,:)*(x(:,i+1:end)*t(i,i+1:end)')+...
-	    t(i,i)*t(1:i,i+1:end)*x(i+1:end,i)+...
-	    t(i,i-1)*t(1:i,i+1:end)*x(i+1:end,i-1);
-	c1 = t(1:i,:)*(x(:,i+1:end)*t(i-1,i+1:end)')+...
-	     t(i-1,i-1)*t(1:i,i+1:end)*x(i+1:end,i-1)+...
-	     t(i-1,i)*t(1:i,i+1:end)*x(i+1:end,i);
-      end
-      q = [eye(i)-t(1:i,1:i)*t(i,i) -t(1:i,1:i)*t(i,i-1);...
-	   -t(1:i,1:i)*t(i-1,i) eye(i)-t(1:i,1:i)*t(i-1,i-1)];
-      z =  q\[b(1:i,i)+c;b(1:i,i-1)+c1];
-      x(1:i,i) = z(1:i);
-      x(1:i,i-1) = z(i+1:end);
-      x(i,1:i-1)=x(1:i-1,i)';
-      x(i-1,1:i-2)=x(1:i-2,i-1)';
-      i = i - 2;
+        if exist('U','var')
+            clear('U','T','k','n')
+        end
+    end
+    
+    u = [];
+  
+    if size(a,1) == 1
+        x=b/(1-a*a);
+        return
+    end
+    
+    if method<2
+        [U,T] = schur(a);
+        if exist('OCTAVE_VERSION') || matlab_ver_less_than('7.0.1')
+            e1 = abs(my_ordeig(T)) > 2-qz_criterium;
+        else
+            e1 = abs(ordeig(T)) > 2-qz_criterium;
+        end
+        k = sum(e1);       % Number of unit roots. 
+        n = length(e1)-k;  % Number of stationary variables.
+        if exist('ordschur','builtin')
+            % Selects stable roots
+            [U,T] = ordschur(U,T,e1);
+            % B = U(:,k+1:end)'*b*U(:,k+1:end);
+            T = T(k+1:end,k+1:end);
+        elseif k > 0
+            % Problem for Matlab version that don't have ordschur
+            error(['lyapunov_sym: you need a Matlab version > 6.5 to handle models' ...
+                   ' with unit roots'])
+            %else
+            % No unit root
+            %n = length(e1);
+            %B=u'*b*u;
+        end
+    end
+    B = U(:,k+1:end)'*b*U(:,k+1:end);    
+    
+    x = zeros(n,n);
+    i = n;
+
+    while i >= 2
+        if T(i,i-1) == 0
+            if i == n
+                c = zeros(n,1);
+            else
+                c = T(1:i,:)*(x(:,i+1:end)*T(i,i+1:end)') + ...
+                    T(i,i)*T(1:i,i+1:end)*x(i+1:end,i);
+            end
+            q = eye(i)-T(1:i,1:i)*T(i,i);
+            x(1:i,i) = q\(B(1:i,i)+c);
+            x(i,1:i-1) = x(1:i-1,i)';
+            i = i - 1;
+        else
+            if i == n
+                c = zeros(n,1);
+                c1 = zeros(n,1);
+            else
+                c = T(1:i,:)*(x(:,i+1:end)*T(i,i+1:end)') + ...
+                    T(i,i)*T(1:i,i+1:end)*x(i+1:end,i) + ...
+                    T(i,i-1)*T(1:i,i+1:end)*x(i+1:end,i-1);
+                c1 = T(1:i,:)*(x(:,i+1:end)*T(i-1,i+1:end)') + ...
+                     T(i-1,i-1)*T(1:i,i+1:end)*x(i+1:end,i-1) + ...
+                     T(i-1,i)*T(1:i,i+1:end)*x(i+1:end,i);
+            end
+            q = [eye(i)-T(1:i,1:i)*T(i,i) - T(1:i,1:i)*T(i,i-1) ; ...
+                 -T(1:i,1:i)*T(i-1,i) eye(i) - T(1:i,1:i)*T(i-1,i-1) ];
+            z =  q\[ B(1:i,i)+c ; B(1:i,i-1) + c1 ];
+            x(1:i,i) = z(1:i);
+            x(1:i,i-1) = z(i+1:end);
+            x(i,1:i-1) = x(1:i-1,i)';
+            x(i-1,1:i-2) = x(1:i-2,i-1)';
+            i = i - 2;
+        end
+    end
+    if i == 1
+        c = T(1,:)*(x(:,2:end)*T(1,2:end)') + T(1,1)*T(1,2:end)*x(2:end,1);
+        x(1,1) = (B(1,1)+c)/(1-T(1,1)*T(1,1));
     end
-  end
-  if i == 1
-    c = t(1,:)*(x(:,2:end)*t(1,2:end)')+t(1,1)*t(1,2:end)*x(2:end,1);
-    x(1,1)=(b(1,1)+c)/(1-t(1,1)*t(1,1));
-  end
-  x = u(:,k+1:end)*x*u(:,k+1:end)';
-  u = u(:,1:k);
\ No newline at end of file
+    x = U(:,k+1:end)*x*U(:,k+1:end)';
+    u = U(:,1:k);
\ No newline at end of file

matlab/th_autocovariances.m

-function [Gamma_y,ivar]=th_autocovariances(dr,ivar,M_,options_)
+function [Gamma_y,ivar] = 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_)
-%  
+% Theoretical HPfiltering is available as an option
+%    
 % INPUTS
-%   dr:         structure of decisions rules for stochastic simulations
-%   ivar:       subset of variables
-%   M_
-%   options_
+%   dr:               [structure]    Reduced form solution of the DSGE model  (decisions rules)
+%   ivar:             [integer]      Vector of indices for a subset of variables.
+%   M_                [structure]    Global dynare's structure, description of the DSGE model.
+%   options_          [structure]    Global dynare's structure.
+%   nodecomposition   [integer]      Scalar, if different from zero the variance decomposition is not triggered.  
 %    
 % OUTPUTS
-%   Gamma_y:    theoritical auto-covariances
-%   ivar:       subset of variables
+%   Gamma_y           [cell]         Matlab cell of nar+3 (second order approximation) or nar+2 (first order approximation) arrays, 
+%                                    where nar is the order of the autocorrelation function.
+%                                      Gamma_y{1}       [double]  Covariance matrix.
+%                                      Gamma_y{i}       [double]  Autocorrelation function (for i=1,...,options_.nar).
+%                                      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.
 %
 % SPECIAL REQUIREMENTS
-%   Theoretical HP filtering is available as an option
+%   
 
-% Copyright (C) 2001-2008 Dynare Team
+% Copyright (C) 2001-2009 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -33,188 +41,181 @@ function [Gamma_y,ivar]=th_autocovariances(dr,ivar,M_,options_)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-  exo_names_orig_ord  = M_.exo_names_orig_ord;
-  if exist('OCTAVE_VERSION')
-    warning('off', 'Octave:divide-by-zero')
-  else
-    warning off MATLAB:dividebyzero
-  end
-  nar = options_.ar;
-  Gamma_y = cell(nar+1,1);
-  if isempty(ivar)
-    ivar = [1:M_.endo_nbr]';
-  end
-  nvar = size(ivar,1);
-  
-  ghx = dr.ghx;
-  ghu = dr.ghu;
-  npred = dr.npred;
-  nstatic = dr.nstatic;
-  kstate = dr.kstate;
-  order = dr.order_var;
-  iv(order) = [1:length(order)];
-  nx = size(ghx,2);
-  
-  ikx = [nstatic+1:nstatic+npred];
-  
-  k0 = kstate(find(kstate(:,2) <= M_.maximum_lag+1),:);
-  i0 = find(k0(:,2) == M_.maximum_lag+1);
-  i00 = i0;
-  n0 = length(i0);
-  AS = ghx(:,i0);
-  ghu1 = zeros(nx,M_.exo_nbr);
-  ghu1(i0,:) = ghu(ikx,:);
-  for i=M_.maximum_lag:-1:2
-    i1 = find(k0(:,2) == i);
-    n1 = size(i1,1);
-    j1 = zeros(n1,1);
-    for k1 = 1:n1
-      j1(k1) = find(k0(i00,1)==k0(i1(k1),1));
-    end
-    AS(:,j1) = AS(:,j1)+ghx(:,i1);
-    i0 = i1;
-  end
-  b = ghu1*M_.Sigma_e*ghu1';
-
-
-  ipred = nstatic+(1:npred)';
-  % state space representation for state variables only
-  [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);
-    iky = iv(ivar);
-    if ~isempty(u)
-      iky = iky(find(all(abs(ghx(iky,:)*u) < options_.Schur_vec_tol,2)));
-      ivar = dr.order_var(iky);
-    end
-    aa = ghx(iky,:);
-    bb = ghu(iky,:);
-    if options_.order == 2         % mean correction for 2nd order
-      Ex = (dr.ghs2(ikx)+dr.ghxx(ikx,:)*vx(:)+dr.ghuu(ikx,:)*M_.Sigma_e(:))/2;
-      Ex = (eye(n0)-AS(ikx,:))\Ex;
-      Gamma_y{nar+3} = AS(iky,:)*Ex+(dr.ghs2(iky)+dr.ghxx(iky,:)*vx(:)+...
-				     dr.ghuu(iky,:)*M_.Sigma_e(:))/2;
+    if nargin<5
+        nodecomposition = 0;
     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;
     
-    % autocorrelations
-    if nar > 0
-      vxy = (A*vx*aa'+ghu1*M_.Sigma_e*bb');
-
-      sy = sqrt(diag(Gamma_y{1}));
-      sy = sy *sy';
-      Gamma_y{2} = aa*vxy./sy;
-      
-      for i=2:nar
-	vxy = A*vxy;
-	Gamma_y{i+1} = aa*vxy./sy;
-      end
+    exo_names_orig_ord  = M_.exo_names_orig_ord;
+    if exist('OCTAVE_VERSION')
+        warning('off', 'Octave:divide-by-zero')
+    else
+        warning off MATLAB:dividebyzero
     end
-    
-    % variance decomposition
-    if M_.exo_nbr > 1
-      Gamma_y{nar+2} = zeros(length(ivar),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;
-      b1 = b1*cs;
-      b2(:,exo_names_orig_ord) = ghu(iky,:);
-      b2 = b2*cs;
-      vx  = lyapunov_symm(A,b1*b1',options_.qz_criterium);
-      vv = diag(aa*vx*aa'+b2*b2');
-      for i=1:M_.exo_nbr
-	vx1 = lyapunov_symm(A,b1(:,i)*b1(:,i)',options_.qz_criterium);
-	Gamma_y{nar+2}(:,i) = abs(diag(aa*vx1*aa'+b2(:,i)*b2(:,i)'))./vv;
-      end
+    nar = options_.ar;
+    Gamma_y = cell(nar+1,1);
+    if isempty(ivar)
+        ivar = [1:M_.endo_nbr]';
     end
-  else
-    if options_.order < 2
-      iky = iv(ivar);  
-      aa = ghx(iky,:);
-      bb = ghu(iky,:);
+    nvar = size(ivar,1);
+  
+    ghx = dr.ghx;
+    ghu = dr.ghu;
+    npred = dr.npred;
+    nstatic = dr.nstatic;
+    kstate = dr.kstate;
+    order = dr.order_var;
+    iv(order) = [1:length(order)];
+    nx = size(ghx,2);
+  
+    ikx = [nstatic+1:nstatic+npred];
+  
+    k0 = kstate(find(kstate(:,2) <= M_.maximum_lag+1),:);
+    i0 = find(k0(:,2) == M_.maximum_lag+1);
+    i00 = i0;
+    n0 = length(i0);
+    AS = ghx(:,i0);
+    ghu1 = zeros(nx,M_.exo_nbr);
+    ghu1(i0,:) = ghu(ikx,:);
+    for i=M_.maximum_lag:-1:2
+        i1 = find(k0(:,2) == i);
+        n1 = size(i1,1);
+        j1 = zeros(n1,1);
+        for k1 = 1:n1
+            j1(k1) = find(k0(i00,1)==k0(i1(k1),1));
+        end
+        AS(:,j1) = AS(:,j1)+ghx(:,i1);
+        i0 = i1;
     end
-    lambda = options_.hp_filter;
-    ngrid = options_.hp_ngrid;
-    freqs = 0 : ((2*pi)/ngrid) : (2*pi*(1 - .5/ngrid)); 
-    tpos  = exp( sqrt(-1)*freqs);
-    tneg  =  exp(-sqrt(-1)*freqs);
-    hp1 = 4*lambda*(1 - cos(freqs)).^2 ./ (1 + 4*lambda*(1 - cos(freqs)).^2);
+    b = ghu1*M_.Sigma_e*ghu1';
     
-    mathp_col = [];
-    IA = eye(size(A,1));
-    IE = eye(M_.exo_nbr);
-    for ig = 1:ngrid
-      f_omega  =(1/(2*pi))*( [inv(IA-A*tneg(ig))*ghu1;IE]...
-			     *M_.Sigma_e*[ghu1'*inv(IA-A'*tpos(ig)) ...
-		    IE]); % state variables
-      g_omega = [aa*tneg(ig) bb]*f_omega*[aa'*tpos(ig); bb']; % selected variables
-      f_hp = hp1(ig)^2*g_omega; % spectral density of selected filtered series
-      mathp_col = [mathp_col ; (f_hp(:))'];    % store as matrix row
-                                               % for ifft
-    end;
-
-    % covariance of filtered series
-    imathp_col = real(ifft(mathp_col))*(2*pi);
-
-    Gamma_y{1} = reshape(imathp_col(1,:),nvar,nvar);
     
-    % autocorrelations
-    if nar > 0
-      sy = sqrt(diag(Gamma_y{1}));
-      sy = sy *sy';
-      for i=1:nar
-	Gamma_y{i+1} = reshape(imathp_col(i+1,:),nvar,nvar)./sy;
-      end
+    ipred = nstatic+(1:npred)';
+    % state space representation for state variables only
+    [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);
+        iky = iv(ivar);
+        if ~isempty(u)
+            iky = iky(find(all(abs(ghx(iky,:)*u) < options_.Schur_vec_tol,2)));
+            ivar = dr.order_var(iky);
+        end
+        aa = ghx(iky,:);
+        bb = ghu(iky,:);
+        if options_.order == 2         % mean correction for 2nd order
+            Ex = (dr.ghs2(ikx)+dr.ghxx(ikx,:)*vx(:)+dr.ghuu(ikx,:)*M_.Sigma_e(:))/2;
+            Ex = (eye(n0)-AS(ikx,:))\Ex;
+            Gamma_y{nar+3} = AS(iky,:)*Ex+(dr.ghs2(iky)+dr.ghxx(iky,:)*vx(:)+...
+                                           dr.ghuu(iky,:)*M_.Sigma_e(:))/2;
+        end
     end
-    
-    %variance decomposition
-    if M_.exo_nbr > 1
-      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)';
-      SS = cs*cs';
-      b1(:,exo_names_orig_ord) = ghu1;
-      b2(:,exo_names_orig_ord) = ghu(iky,:);
-      mathp_col = [];
-      IA = eye(size(A,1));
-      IE = eye(M_.exo_nbr);
-      for ig = 1:ngrid
-	f_omega  =(1/(2*pi))*( [inv(IA-A*tneg(ig))*b1;IE]...
-			       *SS*[b1'*inv(IA-A'*tpos(ig)) ...
-		    IE]); % state variables
-	g_omega = [aa*tneg(ig) b2]*f_omega*[aa'*tpos(ig); b2']; % selected variables
-	f_hp = hp1(ig)^2*g_omega; % spectral density of selected filtered series
-	mathp_col = [mathp_col ; (f_hp(:))'];    % store as matrix row
-						 % for ifft
-      end;
-
-      imathp_col = real(ifft(mathp_col))*(2*pi);
-      vv = diag(reshape(imathp_col(1,:),nvar,nvar));
-      for i=1:M_.exo_nbr
-	mathp_col = [];
-	SSi = cs(:,i)*cs(:,i)';
-	for ig = 1:ngrid
-	  f_omega  =(1/(2*pi))*( [inv(IA-A*tneg(ig))*b1;IE]...
-				 *SSi*[b1'*inv(IA-A'*tpos(ig)) ...
-		    IE]); % state variables
-	  g_omega = [aa*tneg(ig) b2]*f_omega*[aa'*tpos(ig); b2']; % selected variables
-	  f_hp = hp1(ig)^2*g_omega; % spectral density of selected filtered series
-	  mathp_col = [mathp_col ; (f_hp(:))'];    % store as matrix row
-						   % for ifft
-	end;
-