UnivariateSpectralDensity

 This function computes the theoretical spectral density of each
 endogenous variable declared in var_list. Results are stored in 
 oo_ and may be plotted.
 
 Adapted from th_autocovariances.m.

0001 function [omega,f] = UnivariateSpectralDensity(dr,var_list)
0002 % This function computes the theoretical spectral density of each
0003 % endogenous variable declared in var_list. Results are stored in
0004 % oo_ and may be plotted.
0005 %
0006 % Adapted from th_autocovariances.m.
0007 
0008 % Copyright (C) 2006-2011 Dynare Team
0009 %
0010 % This file is part of Dynare.
0011 %
0012 % Dynare is free software: you can redistribute it and/or modify
0013 % it under the terms of the GNU General Public License as published by
0014 % the Free Software Foundation, either version 3 of the License, or
0015 % (at your option) any later version.
0016 %
0017 % Dynare is distributed in the hope that it will be useful,
0018 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0019 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0020 % GNU General Public License for more details.
0021 %
0022 % You should have received a copy of the GNU General Public License
0023 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0024 
0025 global options_ oo_ M_
0026 
0027 
0028 omega = []; f = [];
0029 
0030 if options_.order > 1
0031     disp('UnivariateSpectralDensity :: I Cannot compute the theoretical spectral density') 
0032     disp('with a second order approximation of the DSGE model!')
0033     disp('Please set order = 1. I abort')
0034     return
0035 end
0036 pltinfo  = options_.SpectralDensity.plot; 
0037 cutoff   = options_.SpectralDensity.cutoff; 
0038 sdl      = options_.SpectralDensity.sdl; 
0039 omega    = (0:sdl:pi)';
0040 GridSize = length(omega);
0041 exo_names_orig_ord  = M_.exo_names_orig_ord;
0042 if exist('OCTAVE_VERSION')
0043     warning('off', 'Octave:divide-by-zero')
0044 else
0045     warning off MATLAB:dividebyzero
0046 end
0047 if nargin<2
0048     var_list = [];
0049 end
0050 if size(var_list,1) == 0
0051     var_list = M_.endo_names(1:M_.orig_endo_nbr, :);
0052 end
0053 nvar = size(var_list,1);
0054 ivar=zeros(nvar,1);
0055 for i=1:nvar
0056     i_tmp = strmatch(var_list(i,:),M_.endo_names,'exact');
0057     if isempty(i_tmp)
0058         error (['One of the variables specified does not exist']) ;
0059     else
0060         ivar(i) = i_tmp;
0061     end
0062 end
0063 f = zeros(nvar,GridSize);
0064 ghx = dr.ghx;
0065 ghu = dr.ghu;
0066 npred = dr.npred;
0067 nstatic = dr.nstatic;
0068 kstate = dr.kstate;
0069 order = dr.order_var;
0070 iv(order) = [1:length(order)];
0071 nx = size(ghx,2);
0072 ikx = [nstatic+1:nstatic+npred];
0073 A = zeros(nx,nx);
0074 k0 = kstate(find(kstate(:,2) <= M_.maximum_lag+1),:);
0075 i0 = find(k0(:,2) == M_.maximum_lag+1);
0076 i00 = i0;
0077 n0 = length(i0);
0078 A(i0,:) = ghx(ikx,:);
0079 AS = ghx(:,i0);
0080 ghu1 = zeros(nx,M_.exo_nbr);
0081 ghu1(i0,:) = ghu(ikx,:);
0082 for i=M_.maximum_lag:-1:2
0083     i1 = find(k0(:,2) == i);
0084     n1 = size(i1,1);
0085     j1 = zeros(n1,1);
0086     j2 = j1;
0087     for k1 = 1:n1
0088         j1(k1) = find(k0(i00,1)==k0(i1(k1),1));
0089         j2(k1) = find(k0(i0,1)==k0(i1(k1),1));
0090     end
0091     AS(:,j1) = AS(:,j1)+ghx(:,i1);
0092     i0 = i1;
0093 end
0094 Gamma = zeros(nvar,cutoff+1);
0095 [A,B] = kalman_transition_matrix(dr,ikx',1:nx,M_.exo_nbr);
0096 [vx, u] =  lyapunov_symm(A,B*M_.Sigma_e*B',options_.qz_criterium,options_.lyapunov_complex_threshold);
0097 iky = iv(ivar);
0098 if ~isempty(u)
0099     iky = iky(find(any(abs(ghx(iky,:)*u) < options_.Schur_vec_tol,2)));
0100     ivar = dr.order_var(iky);
0101 end
0102 aa = ghx(iky,:);
0103 bb = ghu(iky,:);
0104 
0105 if options_.hp_filter == 0
0106     tmp = aa*vx*aa'+ bb*M_.Sigma_e*bb';
0107     k = find(abs(tmp) < 1e-12);
0108     tmp(k) = 0;
0109     Gamma(:,1) = diag(tmp);
0110     vxy = (A*vx*aa'+ghu1*M_.Sigma_e*bb');
0111     tmp = aa*vxy;
0112     k = find(abs(tmp) < 1e-12);
0113     tmp(k) = 0;
0114     Gamma(:,2) = diag(tmp);
0115     for i=2:cutoff
0116         vxy = A*vxy;
0117         tmp = aa*vxy;
0118         k = find(abs(tmp) < 1e-12);
0119         tmp(k) = 0;
0120         Gamma(:,i+1) = diag(tmp);
0121     end
0122 else
0123     iky = iv(ivar);  
0124     aa = ghx(iky,:);
0125     bb = ghu(iky,:);
0126     lambda = options_.hp_filter;
0127     ngrid = options_.hp_ngrid;
0128     freqs = 0 : ((2*pi)/ngrid) : (2*pi*(1 - .5/ngrid)); 
0129     tpos  = exp( sqrt(-1)*freqs);
0130     tneg  =  exp(-sqrt(-1)*freqs);
0131     hp1 = 4*lambda*(1 - cos(freqs)).^2 ./ (1 + 4*lambda*(1 - cos(freqs)).^2);
0132     mathp_col = [];
0133     IA = eye(size(A,1));
0134     IE = eye(M_.exo_nbr);
0135     for ig = 1:ngrid
0136         f_omega  =(1/(2*pi))*( [inv(IA-A*tneg(ig))*ghu1;IE]...
0137                                *M_.Sigma_e*[ghu1'*inv(IA-A'*tpos(ig)) IE]); % state variables
0138         g_omega = [aa*tneg(ig) bb]*f_omega*[aa'*tpos(ig); bb']; % selected variables
0139         f_hp = hp1(ig)^2*g_omega; % spectral density of selected filtered series
0140         mathp_col = [mathp_col ; (f_hp(:))'];    % store as matrix row
0141                                                  % for ifft
0142     end;
0143     imathp_col = real(ifft(mathp_col))*(2*pi);
0144     tmp = reshape(imathp_col(1,:),nvar,nvar);
0145     k = find(abs(tmp)<1e-12);
0146     tmp(k) = 0;
0147     Gamma(:,1) = diag(tmp);
0148     sy = sqrt(Gamma(:,1));
0149     sy = sy *sy';
0150     for i=1:cutoff-1
0151         tmp = reshape(imathp_col(i+1,:),nvar,nvar)./sy;
0152         k = find(abs(tmp) < 1e-12);
0153         tmp(k) = 0;
0154         Gamma(:,i+1) = diag(tmp);
0155     end
0156 end
0157 
0158 H = 1:cutoff;
0159 for i=1:nvar
0160     f(i,:) = Gamma(i,1)/(2*pi) + Gamma(i,H+1)*cos(H'*omega')/pi;
0161 end  
0162 
0163 if exist('OCTAVE_VERSION')
0164     warning('on', 'Octave:divide-by-zero')
0165 else
0166     warning on MATLAB:dividebyzero
0167 end
0168 
0169 if pltinfo
0170     for i= 1:nvar
0171         figure('Name',['Spectral Density of ' deblank(M_.endo_names(ivar(i),:)) '.'])
0172         plot(omega,f(i,:),'-k','linewidth',2)
0173         xlabel('0 \leq \omega \leq \pi')
0174         ylabel('f(\omega)')
0175         box on
0176         axis tight
0177     end
0178 end