th_autocovariances

 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]    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           [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+1}     [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.    
   stationary_vars   [integer]      Vector of indices of stationary variables (as a subset of 1:length(ivar))

 SPECIAL REQUIREMENTS

0001 function [Gamma_y,stationary_vars] = th_autocovariances(dr,ivar,M_,options_,nodecomposition)
0002 % Computes the theoretical auto-covariances, Gamma_y, for an AR(p) process
0003 % with coefficients dr.ghx and dr.ghu and shock variances Sigma_e_
0004 % for a subset of variables ivar (indices in lgy_)
0005 % Theoretical HPfiltering is available as an option
0006 %
0007 % INPUTS
0008 %   dr:               [structure]    Reduced form solution of the DSGE model  (decisions rules)
0009 %   ivar:             [integer]      Vector of indices for a subset of variables.
0010 %   M_                [structure]    Global dynare's structure, description of the DSGE model.
0011 %   options_          [structure]    Global dynare's structure.
0012 %   nodecomposition   [integer]      Scalar, if different from zero the variance decomposition is not triggered.
0013 %
0014 % OUTPUTS
0015 %   Gamma_y           [cell]         Matlab cell of nar+3 (second order approximation) or nar+2 (first order approximation) arrays,
0016 %                                    where nar is the order of the autocorrelation function.
0017 %                                      Gamma_y{1}       [double]  Covariance matrix.
0018 %                                      Gamma_y{i+1}     [double]  Autocorrelation function (for i=1,...,options_.nar).
0019 %                                      Gamma_y{nar+2}   [double]  Variance decomposition.
0020 %                                      Gamma_y{nar+3}   [double]  Expectation of the endogenous variables associated with a second
0021 %                                                                 order approximation.
0022 %   stationary_vars   [integer]      Vector of indices of stationary variables (as a subset of 1:length(ivar))
0023 %
0024 % SPECIAL REQUIREMENTS
0025 %
0026 
0027 % Copyright (C) 2001-2011 Dynare Team
0028 %
0029 % This file is part of Dynare.
0030 %
0031 % Dynare is free software: you can redistribute it and/or modify
0032 % it under the terms of the GNU General Public License as published by
0033 % the Free Software Foundation, either version 3 of the License, or
0034 % (at your option) any later version.
0035 %
0036 % Dynare is distributed in the hope that it will be useful,
0037 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0038 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0039 % GNU General Public License for more details.
0040 %
0041 % You should have received a copy of the GNU General Public License
0042 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0043 
0044 if nargin<5
0045     nodecomposition = 0;
0046 end
0047 
0048 if options_.order >= 3
0049     error('Theoretical moments not implemented above 2nd order')
0050 end
0051 
0052 endo_nbr = M_.endo_nbr;
0053 exo_names_orig_ord  = M_.exo_names_orig_ord;
0054 if exist('OCTAVE_VERSION')
0055     warning('off', 'Octave:divide-by-zero')
0056 else
0057     warning off MATLAB:dividebyzero
0058 end
0059 nar = options_.ar;
0060 Gamma_y = cell(nar+1,1);
0061 if isempty(ivar)
0062     ivar = [1:endo_nbr]';
0063 end
0064 nvar = size(ivar,1);
0065 
0066 ghx = dr.ghx;
0067 ghu = dr.ghu;
0068 npred = dr.npred;
0069 nstatic = dr.nstatic;
0070 
0071 nx = size(ghx,2);
0072 if options_.block == 0
0073     %order_var = dr.order_var;
0074     inv_order_var = dr.inv_order_var;
0075     kstate = dr.kstate;
0076     ikx = [nstatic+1:nstatic+npred];
0077     k0 = kstate(find(kstate(:,2) <= M_.maximum_lag+1),:);
0078     i0 = find(k0(:,2) == M_.maximum_lag+1);
0079     i00 = i0;
0080     n0 = length(i0);
0081     AS = ghx(:,i0);
0082     ghu1 = zeros(nx,M_.exo_nbr);
0083     ghu1(i0,:) = ghu(ikx,:);
0084     for i=M_.maximum_lag:-1:2
0085         i1 = find(k0(:,2) == i);
0086         n1 = size(i1,1);
0087         j1 = zeros(n1,1);
0088         for k1 = 1:n1
0089             j1(k1) = find(k0(i00,1)==k0(i1(k1),1));
0090         end
0091         AS(:,j1) = AS(:,j1)+ghx(:,i1);
0092         i0 = i1;
0093     end
0094 else
0095     ghu1 = zeros(nx,M_.exo_nbr);
0096     trend = 1:M_.endo_nbr;
0097     inv_order_var = trend(M_.block_structure.variable_reordered);
0098     ghu1(1:length(dr.state_var),:) = ghu(dr.state_var,:);
0099     npred = npred + dr.nboth;
0100 end;
0101 b = ghu1*M_.Sigma_e*ghu1';
0102 
0103 
0104 if options_.block == 0
0105     ipred = nstatic+(1:npred)';
0106 else
0107     ipred = dr.state_var;
0108 end;
0109 % state space representation for state variables only
0110 [A,B] = kalman_transition_matrix(dr,ipred,1:nx,M_.exo_nbr);
0111 % Compute stationary variables (before HP filtering),
0112 % and compute 2nd order mean correction on stationary variables (in case of
0113 % HP filtering, this mean correction is computed *before* filtering)
0114 if options_.order == 2 || options_.hp_filter == 0
0115     [vx, u] =  lyapunov_symm(A,B*M_.Sigma_e*B',options_.qz_criterium,options_.lyapunov_complex_threshold);
0116     if options_.block == 0
0117         iky = inv_order_var(ivar);
0118     else
0119         iky = ivar;
0120     end;
0121     stationary_vars = (1:length(ivar))';
0122     if ~isempty(u)
0123         x = abs(ghx*u);
0124         iky = iky(find(all(x(iky,:) < options_.Schur_vec_tol,2)));
0125         stationary_vars = find(all(x(inv_order_var(ivar(stationary_vars)),:) < options_.Schur_vec_tol,2));
0126     end
0127     aa = ghx(iky,:);
0128     bb = ghu(iky,:);
0129     if options_.order == 2         % mean correction for 2nd order
0130         Ex = (dr.ghs2(ikx)+dr.ghxx(ikx,:)*vx(:)+dr.ghuu(ikx,:)*M_.Sigma_e(:))/2;
0131         Ex = (eye(n0)-AS(ikx,:))\Ex;
0132         Gamma_y{nar+3} = NaN*ones(nvar, 1);
0133         Gamma_y{nar+3}(stationary_vars) = AS(iky,:)*Ex+(dr.ghs2(iky)+dr.ghxx(iky,:)*vx(:)+...
0134                                                         dr.ghuu(iky,:)*M_.Sigma_e(:))/2;
0135     end
0136 end
0137 if options_.hp_filter == 0
0138     v = NaN*ones(nvar,nvar);
0139     v(stationary_vars,stationary_vars) = aa*vx*aa'+ bb*M_.Sigma_e*bb';
0140     k = find(abs(v) < 1e-12);
0141     v(k) = 0;
0142     Gamma_y{1} = v;
0143     % autocorrelations
0144     if nar > 0
0145         vxy = (A*vx*aa'+ghu1*M_.Sigma_e*bb');    
0146         sy = sqrt(diag(Gamma_y{1}));
0147         sy = sy(stationary_vars);
0148         sy = sy *sy';
0149         Gamma_y{2} = NaN*ones(nvar,nvar);
0150         Gamma_y{2}(stationary_vars,stationary_vars) = aa*vxy./sy;
0151         for i=2:nar
0152             vxy = A*vxy;
0153             Gamma_y{i+1} = NaN*ones(nvar,nvar);
0154             Gamma_y{i+1}(stationary_vars,stationary_vars) = aa*vxy./sy;
0155         end
0156     end
0157     % variance decomposition
0158     if ~nodecomposition && M_.exo_nbr > 1 && size(stationary_vars, 1) > 0
0159         Gamma_y{nar+2} = zeros(nvar,M_.exo_nbr);
0160         SS(exo_names_orig_ord,exo_names_orig_ord)=M_.Sigma_e+1e-14*eye(M_.exo_nbr);
0161         cs = chol(SS)';
0162         b1(:,exo_names_orig_ord) = ghu1;
0163         b1 = b1*cs;
0164         b2(:,exo_names_orig_ord) = ghu(iky,:);
0165         b2 = b2*cs;
0166         vx  = lyapunov_symm(A,b1*b1',options_.qz_criterium,options_.lyapunov_complex_threshold,1);
0167         vv = diag(aa*vx*aa'+b2*b2');
0168         vv2 = 0;
0169         for i=1:M_.exo_nbr
0170             vx1 = lyapunov_symm(A,b1(:,i)*b1(:,i)',options_.qz_criterium,options_.lyapunov_complex_threshold,2);
0171             vx2 = abs(diag(aa*vx1*aa'+b2(:,i)*b2(:,i)'));
0172             Gamma_y{nar+2}(stationary_vars,i) = vx2;
0173             vv2 = vv2 +vx2;
0174         end
0175         if max(abs(vv2-vv)./vv) > 1e-4
0176             warning(['Aggregate variance and sum of variances by shocks ' ...
0177                      'differ by more than 0.01 %'])
0178         end
0179         for i=1:M_.exo_nbr
0180             Gamma_y{nar+2}(stationary_vars,i) = Gamma_y{nar+ ...
0181                                 2}(stationary_vars,i)./vv2;
0182         end
0183     end
0184 else% ==> Theoretical HP filter.
0185     % By construction, all variables are stationary when HP filtered
0186     iky = inv_order_var(ivar);  
0187     stationary_vars = (1:length(ivar))';
0188     aa = ghx(iky,:);
0189     bb = ghu(iky,:);
0190 
0191     lambda = options_.hp_filter;
0192     ngrid = options_.hp_ngrid;
0193     freqs = 0 : ((2*pi)/ngrid) : (2*pi*(1 - .5/ngrid)); 
0194     tpos  = exp( sqrt(-1)*freqs);
0195     tneg  =  exp(-sqrt(-1)*freqs);
0196     hp1 = 4*lambda*(1 - cos(freqs)).^2 ./ (1 + 4*lambda*(1 - cos(freqs)).^2);
0197     mathp_col = [];
0198     IA = eye(size(A,1));
0199     IE = eye(M_.exo_nbr);
0200     for ig = 1:ngrid
0201         f_omega  =(1/(2*pi))*( [inv(IA-A*tneg(ig))*ghu1;IE]...
0202                                *M_.Sigma_e*[ghu1'*inv(IA-A'*tpos(ig)) ...
0203                             IE]); % state variables
0204         g_omega = [aa*tneg(ig) bb]*f_omega*[aa'*tpos(ig); bb']; % selected variables
0205         f_hp = hp1(ig)^2*g_omega; % spectral density of selected filtered series
0206         mathp_col = [mathp_col ; (f_hp(:))'];    % store as matrix row
0207                                                  % for ifft
0208     end;
0209     % Covariance of filtered series
0210     imathp_col = real(ifft(mathp_col))*(2*pi);
0211     Gamma_y{1} = reshape(imathp_col(1,:),nvar,nvar);
0212     % Autocorrelations
0213     if nar > 0
0214         sy = sqrt(diag(Gamma_y{1}));
0215         sy = sy *sy';
0216         for i=1:nar
0217             Gamma_y{i+1} = reshape(imathp_col(i+1,:),nvar,nvar)./sy;
0218         end
0219     end
0220     % Variance decomposition
0221     if ~nodecomposition && M_.exo_nbr > 1
0222         Gamma_y{nar+2} = zeros(nvar,M_.exo_nbr);
0223         SS(exo_names_orig_ord,exo_names_orig_ord) = M_.Sigma_e+1e-14*eye(M_.exo_nbr);
0224         cs = chol(SS)';
0225         SS = cs*cs';
0226         b1(:,exo_names_orig_ord) = ghu1;
0227         b2(:,exo_names_orig_ord) = ghu(iky,:);
0228         mathp_col = [];
0229         IA = eye(size(A,1));
0230         IE = eye(M_.exo_nbr);
0231         for ig = 1:ngrid
0232             f_omega  =(1/(2*pi))*( [inv(IA-A*tneg(ig))*b1;IE]...
0233                                    *SS*[b1'*inv(IA-A'*tpos(ig)) ...
0234                                 IE]); % state variables
0235             g_omega = [aa*tneg(ig) b2]*f_omega*[aa'*tpos(ig); b2']; % selected variables
0236             f_hp = hp1(ig)^2*g_omega; % spectral density of selected filtered series
0237             mathp_col = [mathp_col ; (f_hp(:))'];    % store as matrix row
0238                                                      % for ifft
0239         end;  
0240         imathp_col = real(ifft(mathp_col))*(2*pi);
0241         vv = diag(reshape(imathp_col(1,:),nvar,nvar));
0242         for i=1:M_.exo_nbr
0243             mathp_col = [];
0244             SSi = cs(:,i)*cs(:,i)';
0245             for ig = 1:ngrid
0246                 f_omega  =(1/(2*pi))*( [inv(IA-A*tneg(ig))*b1;IE]...
0247                                        *SSi*[b1'*inv(IA-A'*tpos(ig)) ...
0248                                     IE]); % state variables
0249                 g_omega = [aa*tneg(ig) b2]*f_omega*[aa'*tpos(ig); b2']; % selected variables
0250                 f_hp = hp1(ig)^2*g_omega; % spectral density of selected filtered series
0251                 mathp_col = [mathp_col ; (f_hp(:))'];    % store as matrix row
0252                                                          % for ifft
0253             end;
0254             imathp_col = real(ifft(mathp_col))*(2*pi);
0255             Gamma_y{nar+2}(:,i) = abs(diag(reshape(imathp_col(1,:),nvar,nvar)))./vv;
0256         end
0257     end
0258 end
0259 if exist('OCTAVE_VERSION')
0260     warning('on', 'Octave:divide-by-zero')
0261 else
0262     warning on MATLAB:dividebyzero
0263 end