1function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dSPECTRUM_NO_MEAN, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params, M, oo, options, options_ident, indpmodel, indpstderr, indpcorr, indvobs) 2% function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params, M, oo, options, options_ident, indpmodel, indpstderr, indpcorr, indvobs) 3% previously getJJ.m in Dynare 4.5 4% Sets up the Jacobians needed for identification analysis 5% ========================================================================= 6% INPUTS 7% estim_params: [structure] storing the estimation information 8% M: [structure] storing the model information 9% oo: [structure] storing the reduced-form solution results 10% options: [structure] storing the options 11% options_ident: [structure] storing the options for identification 12% indpmodel: [modparam_nbr by 1] index of estimated parameters in M_.params; 13% corresponds to model parameters (no stderr and no corr) 14% in estimated_params block; if estimated_params block is 15% not available, then all model parameters are selected 16% indpstderr: [stderrparam_nbr by 1] index of estimated standard errors, 17% i.e. for all exogenous variables where "stderr" is given 18% in the estimated_params block; if estimated_params block 19% is not available, then all stderr parameters are selected 20% indpcorr: [corrparam_nbr by 2] matrix of estimated correlations, 21% i.e. for all exogenous variables where "corr" is given 22% in the estimated_params block; if estimated_params block 23% is not available, then no corr parameters are selected 24% indvobs: [obs_nbr by 1] index of observed (VAROBS) variables 25% ------------------------------------------------------------------------- 26% OUTPUTS 27% 28% MEAN [endo_nbr by 1], in DR order. Expectation of all model variables 29% * order==1: corresponds to steady state 30% * order==2|3: corresponds to mean computed from pruned state space system (as in Andreasen, Fernandez-Villaverde, Rubio-Ramirez, 2018) 31% dMEAN [endo_nbr by totparam_nbr], in DR Order, Jacobian (wrt all params) of MEAN 32% 33% REDUCEDFORM [rowredform_nbr by 1] in DR order. Steady state and reduced-form model solution matrices for all model variables 34% * order==1: [Yss' vec(ghx)' vech(ghu*Sigma_e*ghu')']', 35% where rowredform_nbr = endo_nbr*(1+nspred+(endo_nbr+1)/2) 36% * order==2: [Yss' vec(ghx)' vech(ghu*Sigma_e*ghu')' vec(ghxx)' vec(ghxu)' vec(ghuu)' vec(ghs2)']', 37% where rowredform_nbr = endo_nbr*(1+nspred+(endo_nbr+1)/2+nspred^2+nspred*exo_nr+exo_nbr^2+1) 38% * order==3: [Yss' vec(ghx)' vech(ghu*Sigma_e*ghu')' vec(ghxx)' vec(ghxu)' vec(ghuu)' vec(ghs2)' vec(ghxxx)' vec(ghxxu)' vec(ghxuu)' vec(ghuuu)' vec(ghxss)' vec(ghuss)']', 39% where rowredform_nbr = endo_nbr*(1+nspred+(endo_nbr+1)/2+nspred^2+nspred*exo_nr+exo_nbr^2+1+nspred^3+nspred^2*exo_nbr+nspred*exo_nbr^2+exo_nbr^3+nspred+exo_nbr) 40% dREDUCEDFORM: [rowredform_nbr by totparam_nbr] in DR order, Jacobian (wrt all params) of REDUCEDFORM 41% * order==1: corresponds to Iskrev (2010)'s J_2 matrix 42% * order==2: corresponds to Mutschler (2015)'s J matrix 43% 44% DYNAMIC [rowdyn_nbr by 1] in declaration order. Steady state and dynamic model derivatives for all model variables 45% * order==1: [ys' vec(g1)']', rowdyn_nbr=endo_nbr+length(g1) 46% * order==2: [ys' vec(g1)' vec(g2)']', rowdyn_nbr=endo_nbr+length(g1)+length(g2) 47% * order==3: [ys' vec(g1)' vec(g2)' vec(g3)']', rowdyn_nbr=endo_nbr+length(g1)+length(g2)+length(g3) 48% dDYNAMIC [rowdyn_nbr by modparam_nbr] in declaration order. Jacobian (wrt model parameters) of DYNAMIC 49% * order==1: corresponds to Ratto and Iskrev (2011)'s J_\Gamma matrix (or LRE) 50% 51% MOMENTS: [obs_nbr+obs_nbr*(obs_nbr+1)/2+nlags*obs_nbr^2 by 1] in DR order. First two theoretical moments for VAROBS variables, i.e. 52% [E[varobs]' vech(E[varobs*varobs'])' vec(E[varobs*varobs(-1)'])' ... vec(E[varobs*varobs(-nlag)'])'] 53% dMOMENTS: [obs_nbr+obs_nbr*(obs_nbr+1)/2+nlags*obs_nbr^2 by totparam_nbr] in DR order. Jacobian (wrt all params) of MOMENTS 54% * order==1: corresponds to Iskrev (2010)'s J matrix 55% * order==2: corresponds to Mutschler (2015)'s \bar{M}_2 matrix, i.e. theoretical moments from the pruned state space system 56% 57% dSPECTRUM: [totparam_nbr by totparam_nbr] in DR order. Gram matrix of Jacobian (wrt all params) of mean and of spectral density for VAROBS variables, where 58% spectral density at frequency w: f(w) = (2*pi)^(-1)*H(exp(-i*w))*E[Inov*Inov']*ctranspose(H(exp(-i*w)) with H being the Transfer function 59% dSPECTRUM = dMEAN*dMEAN + int_{-\pi}^\pi transpose(df(w)/dp')*(df(w)/dp') dw 60% * order==1: corresponds to Qu and Tkachenko (2012)'s G matrix, where Inov and H are computed from linear state space system 61% * order==2: corresponds to Mutschler (2015)'s G_2 matrix, where Inov and H are computed from second-order pruned state space system 62% * order==3: Inov and H are computed from third-order pruned state space system 63% 64% dSPECTRUM_NO_MEAN:[totparam_nbr by totparam_nbr] in DR order. Gram matrix of Jacobian (wrt all params) of spectral density for VAROBS variables, where 65% spectral density at frequency w: f(w) = (2*pi)^(-1)*H(exp(-i*w))*E[Inov*Inov']*ctranspose(H(exp(-i*w)) with H being the Transfer function 66% dSPECTRUM = int_{-\pi}^\pi transpose(df(w)/dp')*(df(w)/dp') dw 67% * order==1: corresponds to Qu and Tkachenko (2012)'s G matrix, where Inov and H are computed from linear state space system 68% * order==2: corresponds to Mutschler (2015)'s G_2 matrix, where Inov and H are computed from second-order pruned state space system 69% * order==3: Inov and H are computed from third-order pruned state space system 70% 71% dMINIMAL: [obs_nbr+minx_nbr^2+minx_nbr*exo_nbr+obs_nbr*minx_nbr+obs_nbr*exo_nbr+exo_nbr*(exo_nbr+1)/2 by totparam_nbr+minx_nbr^2+exo_nbr^2] 72% Jacobian (wrt all params, and similarity_transformation_matrices (T and U)) of observational equivalent minimal ABCD system, 73% corresponds to Komunjer and Ng (2011)'s Deltabar matrix, where 74% MINIMAL = [vec(E[varobs]' vec(minA)' vec(minB)' vec(minC)' vec(minD)' vech(Sigma_e)']' 75% minA, minB, minC and minD is the minimal state space system computed in get_minimal_state_representation 76% * order==1: E[varobs] is equal to steady state 77% * order==2|3: E[varobs] is computed from the pruned state space system (second|third-order accurate), as noted in section 5 of Komunjer and Ng (2011) 78% 79% derivatives_info [structure] for use in dsge_likelihood to compute Hessian analytically. Only used at order==1. 80% Contains dA, dB, and d(B*Sigma_e*B'), where A and B are Kalman filter transition matrice. 81% 82% ------------------------------------------------------------------------- 83% This function is called by 84% * identification_analysis.m 85% ------------------------------------------------------------------------- 86% This function calls 87% * commutation 88% * get_minimal_state_representation 89% * duplication 90% * dyn_vech 91% * fjaco 92% * get_perturbation_params_derivs (previously getH) 93% * get_all_parameters 94% * identification_numerical_objective (previously thet2tau) 95% * pruned_state_space_system 96% * vec 97% ========================================================================= 98% Copyright (C) 2010-2020 Dynare Team 99% 100% This file is part of Dynare. 101% 102% Dynare is free software: you can redistribute it and/or modify 103% it under the terms of the GNU General Public License as published by 104% the Free Software Foundation, either version 3 of the License, or 105% (at your option) any later version. 106% 107% Dynare is distributed in the hope that it will be useful, 108% but WITHOUT ANY WARRANTY; without even the implied warranty of 109% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 110% GNU General Public License for more details. 111% 112% You should have received a copy of the GNU General Public License 113% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 114% ========================================================================= 115 116%get fields from options_ident 117no_identification_moments = options_ident.no_identification_moments; 118no_identification_minimal = options_ident.no_identification_minimal; 119no_identification_spectrum = options_ident.no_identification_spectrum; 120order = options_ident.order; 121nlags = options_ident.ar; 122useautocorr = options_ident.useautocorr; 123grid_nbr = options_ident.grid_nbr; 124kronflag = options_ident.analytic_derivation_mode; 125 126% get fields from M 127endo_nbr = M.endo_nbr; 128exo_nbr = M.exo_nbr; 129fname = M.fname; 130lead_lag_incidence = M.lead_lag_incidence; 131nspred = M.nspred; 132nstatic = M.nstatic; 133params = M.params; 134Sigma_e = M.Sigma_e; 135stderr_e = sqrt(diag(Sigma_e)); 136 137% set all selected values: stderr and corr come first, then model parameters 138xparam1 = get_all_parameters(estim_params, M); %try using estimated_params block 139if isempty(xparam1) 140 %if there is no estimated_params block, consider all stderr and all model parameters, but no corr parameters 141 xparam1 = [stderr_e', params']; 142end 143 144%get numbers/lengths of vectors 145modparam_nbr = length(indpmodel); 146stderrparam_nbr = length(indpstderr); 147corrparam_nbr = size(indpcorr,1); 148totparam_nbr = stderrparam_nbr + corrparam_nbr + modparam_nbr; 149obs_nbr = length(indvobs); 150d2flag = 0; % do not compute second parameter derivatives 151 152% Get Jacobians (wrt selected params) of steady state, dynamic model derivatives and perturbation solution matrices for all endogenous variables 153oo.dr.derivs = get_perturbation_params_derivs(M, options, estim_params, oo, indpmodel, indpstderr, indpcorr, d2flag); 154 155[I,~] = find(lead_lag_incidence'); %I is used to select nonzero columns of the Jacobian of endogenous variables in dynamic model files 156yy0 = oo.dr.ys(I); %steady state of dynamic (endogenous and auxiliary variables) in lead_lag_incidence order 157Yss = oo.dr.ys(oo.dr.order_var); % steady state in DR order 158if order == 1 159 [~, g1 ] = feval([fname,'.dynamic'], yy0, oo.exo_steady_state', params, oo.dr.ys, 1); 160 %g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order 161 DYNAMIC = [Yss; 162 vec(g1(oo.dr.order_var,:))]; %add steady state and put rows of g1 in DR order 163 dDYNAMIC = [oo.dr.derivs.dYss; 164 reshape(oo.dr.derivs.dg1(oo.dr.order_var,:,:),size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2),size(oo.dr.derivs.dg1,3)) ]; %reshape dg1 in DR order and add steady state 165 REDUCEDFORM = [Yss; 166 vec(oo.dr.ghx); 167 dyn_vech(oo.dr.ghu*Sigma_e*transpose(oo.dr.ghu))]; %in DR order 168 dREDUCEDFORM = zeros(endo_nbr*nspred+endo_nbr*(endo_nbr+1)/2, totparam_nbr); 169 for j=1:totparam_nbr 170 dREDUCEDFORM(:,j) = [vec(oo.dr.derivs.dghx(:,:,j)); 171 dyn_vech(oo.dr.derivs.dOm(:,:,j))]; 172 end 173 dREDUCEDFORM = [ [zeros(endo_nbr, stderrparam_nbr+corrparam_nbr) oo.dr.derivs.dYss]; dREDUCEDFORM ]; % add steady state 174 175elseif order == 2 176 [~, g1, g2 ] = feval([fname,'.dynamic'], yy0, oo.exo_steady_state', params, oo.dr.ys, 1); 177 %g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order 178 %g2 is [endo_nbr by yy0ex0_nbr^2] second derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order 179 DYNAMIC = [Yss; 180 vec(g1(oo.dr.order_var,:)); 181 vec(g2(oo.dr.order_var,:))]; %add steady state and put rows of g1 and g2 in DR order 182 dDYNAMIC = [oo.dr.derivs.dYss; 183 reshape(oo.dr.derivs.dg1(oo.dr.order_var,:,:),size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2),size(oo.dr.derivs.dg1,3)); %reshape dg1 in DR order 184 reshape(oo.dr.derivs.dg2(oo.dr.order_var,:),size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2)^2,size(oo.dr.derivs.dg1,3))]; %reshape dg2 in DR order 185 REDUCEDFORM = [Yss; 186 vec(oo.dr.ghx); 187 dyn_vech(oo.dr.ghu*Sigma_e*transpose(oo.dr.ghu)); 188 vec(oo.dr.ghxx); 189 vec(oo.dr.ghxu); 190 vec(oo.dr.ghuu); 191 vec(oo.dr.ghs2)]; %in DR order 192 dREDUCEDFORM = zeros(endo_nbr*nspred+endo_nbr*(endo_nbr+1)/2+endo_nbr*nspred^2+endo_nbr*nspred*exo_nbr+endo_nbr*exo_nbr^2+endo_nbr, totparam_nbr); 193 for j=1:totparam_nbr 194 dREDUCEDFORM(:,j) = [vec(oo.dr.derivs.dghx(:,:,j)); 195 dyn_vech(oo.dr.derivs.dOm(:,:,j)); 196 vec(oo.dr.derivs.dghxx(:,:,j)); 197 vec(oo.dr.derivs.dghxu(:,:,j)); 198 vec(oo.dr.derivs.dghuu(:,:,j)); 199 vec(oo.dr.derivs.dghs2(:,j))]; 200 end 201 dREDUCEDFORM = [ [zeros(endo_nbr, stderrparam_nbr+corrparam_nbr) oo.dr.derivs.dYss]; dREDUCEDFORM ]; % add steady state 202elseif order == 3 203 [~, g1, g2, g3 ] = feval([fname,'.dynamic'], yy0, oo.exo_steady_state', params, oo.dr.ys, 1); 204 %g1 is [endo_nbr by yy0ex0_nbr first derivative (wrt all dynamic variables) of dynamic model equations, i.e. df/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order 205 %g2 is [endo_nbr by yy0ex0_nbr^2] second derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order 206 DYNAMIC = [Yss; 207 vec(g1(oo.dr.order_var,:)); 208 vec(g2(oo.dr.order_var,:)); 209 vec(g3(oo.dr.order_var,:))]; %add steady state and put rows of g1 and g2 in DR order 210 dDYNAMIC = [oo.dr.derivs.dYss; 211 reshape(oo.dr.derivs.dg1(oo.dr.order_var,:,:),size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2),size(oo.dr.derivs.dg1,3)); %reshape dg1 in DR order 212 reshape(oo.dr.derivs.dg2(oo.dr.order_var,:),size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2)^2,size(oo.dr.derivs.dg1,3)); 213 reshape(oo.dr.derivs.dg2(oo.dr.order_var,:),size(oo.dr.derivs.dg1,1)*size(oo.dr.derivs.dg1,2)^2,size(oo.dr.derivs.dg1,3))]; %reshape dg3 in DR order 214 REDUCEDFORM = [Yss; 215 vec(oo.dr.ghx); 216 dyn_vech(oo.dr.ghu*Sigma_e*transpose(oo.dr.ghu)); 217 vec(oo.dr.ghxx); vec(oo.dr.ghxu); vec(oo.dr.ghuu); vec(oo.dr.ghs2); 218 vec(oo.dr.ghxxx); vec(oo.dr.ghxxu); vec(oo.dr.ghxuu); vec(oo.dr.ghuuu); vec(oo.dr.ghxss); vec(oo.dr.ghuss)]; %in DR order 219 dREDUCEDFORM = zeros(size(REDUCEDFORM,1)-endo_nbr, totparam_nbr); 220 for j=1:totparam_nbr 221 dREDUCEDFORM(:,j) = [vec(oo.dr.derivs.dghx(:,:,j)); 222 dyn_vech(oo.dr.derivs.dOm(:,:,j)); 223 vec(oo.dr.derivs.dghxx(:,:,j)); vec(oo.dr.derivs.dghxu(:,:,j)); vec(oo.dr.derivs.dghuu(:,:,j)); vec(oo.dr.derivs.dghs2(:,j)) 224 vec(oo.dr.derivs.dghxxx(:,:,j)); vec(oo.dr.derivs.dghxxu(:,:,j)); vec(oo.dr.derivs.dghxuu(:,:,j)); vec(oo.dr.derivs.dghuuu(:,:,j)); vec(oo.dr.derivs.dghxss(:,:,j)); vec(oo.dr.derivs.dghuss(:,:,j))]; 225 end 226 dREDUCEDFORM = [ [zeros(endo_nbr, stderrparam_nbr+corrparam_nbr) oo.dr.derivs.dYss]; dREDUCEDFORM ]; % add steady state 227end 228 229% Get (pruned) state space representation: 230pruned = pruned_state_space_system(M, options, oo.dr, indvobs, nlags, useautocorr, 1); 231MEAN = pruned.E_y; 232dMEAN = pruned.dE_y; 233%storage for Jacobians used in dsge_likelihood.m for analytical Gradient and Hession of likelihood (only at order=1) 234derivatives_info = struct(); 235if order == 1 236 dT = zeros(endo_nbr,endo_nbr,totparam_nbr); 237 dT(:,(nstatic+1):(nstatic+nspred),:) = oo.dr.derivs.dghx; 238 derivatives_info.DT = dT; 239 derivatives_info.DOm = oo.dr.derivs.dOm; 240 derivatives_info.DYss = oo.dr.derivs.dYss; 241end 242 243%% Compute dMOMENTS 244if ~no_identification_moments 245 E_yy = pruned.Var_y; dE_yy = pruned.dVar_y; 246 if useautocorr 247 E_yyi = pruned.Corr_yi; dE_yyi = pruned.dCorr_yi; 248 else 249 E_yyi = pruned.Var_yi; dE_yyi = pruned.dVar_yi; 250 end 251 MOMENTS = [MEAN; dyn_vech(E_yy)]; 252 for i=1:nlags 253 MOMENTS = [MOMENTS; vec(E_yyi(:,:,i))]; 254 end 255 256 if kronflag == -1 257 %numerical derivative of autocovariogram 258 dMOMENTS = fjaco(str2func('identification_numerical_objective'), xparam1, 1, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvobs, useautocorr, nlags, grid_nbr); %[outputflag=1] 259 dMOMENTS = [dMEAN; dMOMENTS]; %add Jacobian of steady state of VAROBS variables 260 else 261 dMOMENTS = zeros(obs_nbr + obs_nbr*(obs_nbr+1)/2 + nlags*obs_nbr^2 , totparam_nbr); 262 dMOMENTS(1:obs_nbr,:) = dMEAN; %add Jacobian of first moments of VAROBS variables 263 for jp = 1:totparam_nbr 264 dMOMENTS(obs_nbr+1 : obs_nbr+obs_nbr*(obs_nbr+1)/2 , jp) = dyn_vech(dE_yy(:,:,jp)); 265 for i = 1:nlags 266 dMOMENTS(obs_nbr + obs_nbr*(obs_nbr+1)/2 + (i-1)*obs_nbr^2 + 1 : obs_nbr + obs_nbr*(obs_nbr+1)/2 + i*obs_nbr^2, jp) = vec(dE_yyi(:,:,i,jp)); 267 end 268 end 269 end 270else 271 MOMENTS = []; 272 dMOMENTS = []; 273end 274 275%% Compute dSPECTRUM 276%Note that our state space system for computing the spectrum is the following: 277% zhat = A*zhat(-1) + B*xi, where zhat = z - E(z) 278% yhat = C*zhat(-1) + D*xi, where yhat = y - E(y) 279if ~no_identification_spectrum 280 %Some info on the spectral density: Dynare's state space system is given for states (z_{t} = A*z_{t-1} + B*u_{t}) and observables (y_{t} = C*z_{t-1} + D*u_{t}) 281 %The spectral density for y_{t} can be computed using different methods, which are numerically equivalent 282 %See the following code example for the different ways; 283% freqs = 2*pi*(-(grid_nbr/2):1:(grid_nbr/2))'/grid_nbr; %divides the interval [-pi;pi] into ngrid+1 points 284% tpos = exp( sqrt(-1)*freqs); %positive Fourier frequencies 285% tneg = exp(-sqrt(-1)*freqs); %negative Fourier frequencies 286% IA = eye(size(A,1)); 287% IE = eye(exo_nbr); 288% mathp_col1 = NaN(length(freqs),obs_nbr^2); mathp_col2 = mathp_col1; mathp_col3 = mathp_col1; mathp_col4 = mathp_col1; 289% for ig = 1:length(freqs) 290% %method 1: as in UnivariateSpectralDensity.m 291% f_omega =(1/(2*pi))*( [(IA-A*tneg(ig))\B;IE]*Sigma_e*[B'/(IA-A'*tpos(ig)) IE]); % state variables 292% g_omega1 = [C*tneg(ig) D]*f_omega*[C'*tpos(ig); D']; % selected variables 293% %method 2: as in UnivariateSpectralDensity.m but simplified algebraically 294% g_omega2 = (1/(2*pi))*( C*((tpos(ig)*IA-A)\(B*Sigma_e*B'))*((tneg(ig)*IA-A')\(C')) + D*Sigma_e*B'*((tneg(ig)*IA-A')\(C')) + C* ((tpos(ig)*IA-A)\(B*Sigma_e*D')) + D*Sigma_e*D' ); 295% %method 3: use transfer function note that ' is the complex conjugate transpose operator i.e. transpose(ffneg')==ffpos 296% Transferfct = D+C*((tpos(ig)*IA-A)\B); 297% g_omega3 = (1/(2*pi))*(Transferfct*Sigma_e*Transferfct'); 298% %method 4: kronecker products 299% g_omega4 = (1/(2*pi))*( kron( D+C*((tneg(ig)^(-1)*IA-A)\B) , D+C*((tneg(ig)*IA-A)\B) )*Sigma_e(:)); 300% % store as matrix row 301% mathp_col1(ig,:) = (g_omega1(:))'; mathp_col2(ig,:) = (g_omega2(:))'; mathp_col3(ig,:) = (g_omega3(:))'; mathp_col4(ig,:) = g_omega4; 302% end 303% disp([norm(mathp_col1 - mathp_col2); norm(mathp_col1 - mathp_col3); norm(mathp_col1 - mathp_col4); norm(mathp_col2 - mathp_col3); norm(mathp_col2 - mathp_col4); norm(mathp_col3 - mathp_col4);]) 304 %In the following we focus on method 3 305 %Symmetry: 306 % Note that for the compuation of the G matrix we focus only on positive Fourier frequencies due to symmetry of the real part of the spectral density and, hence, the G matrix (which is real by construction). 307 % E.g. if grid_nbr=4, then we subdivide the intervall [-pi;pi] into [-3.1416;-1.5708;0;1.5708;3.1416], but focus only on [0;1.5708;3.1416] for the computation of the G matrix, 308 % keeping in mind that the frequencies [1.5708;3.1416] need to be added twice, whereas the 0 frequency is only added once. 309 freqs = (0 : pi/(grid_nbr/2):pi); % we focus only on positive frequencies 310 tpos = exp( sqrt(-1)*freqs); %positive Fourier frequencies 311 tneg = exp(-sqrt(-1)*freqs); %negative Fourier frequencies 312 IA = eye(size(pruned.A,1)); 313 if kronflag == -1 314 %numerical derivative of spectral density 315 dOmega_tmp = fjaco(str2func('identification_numerical_objective'), xparam1, 2, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvobs, useautocorr, nlags, grid_nbr); %[outputflag=2] 316 kk = 0; 317 for ig = 1:length(freqs) 318 kk = kk+1; 319 dOmega = dOmega_tmp(1 + (kk-1)*obs_nbr^2 : kk*obs_nbr^2,:); 320 if ig == 1 % add zero frequency once 321 dSPECTRUM_NO_MEAN = dOmega'*dOmega; 322 else % due to symmetry to negative frequencies we can add positive frequencies twice 323 dSPECTRUM_NO_MEAN = dSPECTRUM_NO_MEAN + 2*(dOmega'*dOmega); 324 end 325 end 326 elseif kronflag == 1 327 %use Kronecker products 328 dA = reshape(pruned.dA,size(pruned.dA,1)*size(pruned.dA,2),size(pruned.dA,3)); 329 dB = reshape(pruned.dB,size(pruned.dB,1)*size(pruned.dB,2),size(pruned.dB,3)); 330 dC = reshape(pruned.dC,size(pruned.dC,1)*size(pruned.dC,2),size(pruned.dC,3)); 331 dD = reshape(pruned.dD,size(pruned.dD,1)*size(pruned.dD,2),size(pruned.dD,3)); 332 dVarinov = reshape(pruned.dVarinov,size(pruned.dVarinov,1)*size(pruned.dVarinov,2),size(pruned.dVarinov,3)); 333 K_obs_exo = commutation(obs_nbr,size(pruned.Varinov,1)); 334 for ig=1:length(freqs) 335 z = tneg(ig); 336 zIminusA = (z*IA - pruned.A); 337 zIminusAinv = zIminusA\IA; 338 Transferfct = pruned.D + pruned.C*zIminusAinv*pruned.B; % Transfer function 339 dzIminusA = -dA; 340 dzIminusAinv = kron(-(transpose(zIminusA)\IA),zIminusAinv)*dzIminusA; %this takes long 341 dTransferfct = dD + DerivABCD(pruned.C,dC,zIminusAinv,dzIminusAinv,pruned.B,dB); %this takes long 342 dTransferfct_conjt = K_obs_exo*conj(dTransferfct); 343 dOmega = (1/(2*pi))*DerivABCD(Transferfct,dTransferfct,pruned.Varinov,dVarinov,Transferfct',dTransferfct_conjt); %also long 344 if ig == 1 % add zero frequency once 345 dSPECTRUM_NO_MEAN = dOmega'*dOmega; 346 else % due to symmetry to negative frequencies we can add positive frequencies twice 347 dSPECTRUM_NO_MEAN = dSPECTRUM_NO_MEAN + 2*(dOmega'*dOmega); 348 end 349 end 350 elseif (kronflag==0) || (kronflag==-2) 351 for ig = 1:length(freqs) 352 IzminusA = tpos(ig)*IA - pruned.A; 353 invIzminusA = IzminusA\eye(size(pruned.A,1)); 354 Transferfct = pruned.D + pruned.C*invIzminusA*pruned.B; 355 dOmega = zeros(obs_nbr^2,totparam_nbr); 356 for j = 1:totparam_nbr 357 if j <= stderrparam_nbr+corrparam_nbr %stderr and corr parameters: only dSig is nonzero 358 dOmega_tmp = Transferfct*pruned.dVarinov(:,:,j)*Transferfct'; 359 else %model parameters 360 dinvIzminusA = -invIzminusA*(-pruned.dA(:,:,j))*invIzminusA; 361 dTransferfct = pruned.dD(:,:,j) + pruned.dC(:,:,j)*invIzminusA*pruned.B + pruned.C*dinvIzminusA*pruned.B + pruned.C*invIzminusA*pruned.dB(:,:,j); 362 dOmega_tmp = dTransferfct*pruned.Varinov*Transferfct' + Transferfct*pruned.dVarinov(:,:,j)*Transferfct' + Transferfct*pruned.Varinov*dTransferfct'; 363 end 364 dOmega(:,j) = (1/(2*pi))*dOmega_tmp(:); 365 end 366 if ig == 1 % add zero frequency once 367 dSPECTRUM_NO_MEAN = dOmega'*dOmega; 368 else % due to symmetry to negative frequencies we can add positive frequencies twice 369 dSPECTRUM_NO_MEAN = dSPECTRUM_NO_MEAN + 2*(dOmega'*dOmega); 370 end 371 end 372 end 373 % Normalize Matrix and add steady state Jacobian, note that G is real and symmetric by construction 374 dSPECTRUM_NO_MEAN = real(2*pi*dSPECTRUM_NO_MEAN./(2*length(freqs)-1)); 375 dSPECTRUM = dSPECTRUM_NO_MEAN + dMEAN'*dMEAN; 376else 377 dSPECTRUM_NO_MEAN = []; 378 dSPECTRUM = []; 379end 380 381%% Compute dMINIMAL 382if ~no_identification_minimal 383 if obs_nbr < exo_nbr 384 % Check whether criteria can be used 385 warning_KomunjerNg = 'WARNING: Komunjer and Ng (2011) failed:\n'; 386 warning_KomunjerNg = [warning_KomunjerNg ' There are more shocks and measurement errors than observables, this is not implemented (yet).\n']; 387 warning_KomunjerNg = [warning_KomunjerNg ' Skip identification analysis based on minimal state space system.\n']; 388 fprintf(warning_KomunjerNg); 389 dMINIMAL = []; 390 else 391 % Derive and check minimal state vector of first-order 392 SYS.A = oo.dr.ghx(pruned.indx,:); 393 SYS.dA = oo.dr.derivs.dghx(pruned.indx,:,:); 394 SYS.B = oo.dr.ghu(pruned.indx,:); 395 SYS.dB = oo.dr.derivs.dghu(pruned.indx,:,:); 396 SYS.C = oo.dr.ghx(pruned.indy,:); 397 SYS.dC = oo.dr.derivs.dghx(pruned.indy,:,:); 398 SYS.D = oo.dr.ghu(pruned.indy,:); 399 SYS.dD = oo.dr.derivs.dghu(pruned.indy,:,:); 400 [CheckCO,minnx,SYS] = get_minimal_state_representation(SYS,1); 401 402 if CheckCO == 0 403 warning_KomunjerNg = 'WARNING: Komunjer and Ng (2011) failed:\n'; 404 warning_KomunjerNg = [warning_KomunjerNg ' Conditions for minimality are not fullfilled:\n']; 405 warning_KomunjerNg = [warning_KomunjerNg ' Skip identification analysis based on minimal state space system.\n']; 406 fprintf(warning_KomunjerNg); %use sprintf to have line breaks 407 dMINIMAL = []; 408 else 409 minA = SYS.A; dminA = SYS.dA; 410 minB = SYS.B; dminB = SYS.dB; 411 minC = SYS.C; dminC = SYS.dC; 412 minD = SYS.D; dminD = SYS.dD; 413 %reshape into Magnus-Neudecker Jacobians, i.e. dvec(X)/dp 414 dminA = reshape(dminA,size(dminA,1)*size(dminA,2),size(dminA,3)); 415 dminB = reshape(dminB,size(dminB,1)*size(dminB,2),size(dminB,3)); 416 dminC = reshape(dminC,size(dminC,1)*size(dminC,2),size(dminC,3)); 417 dminD = reshape(dminD,size(dminD,1)*size(dminD,2),size(dminD,3)); 418 dvechSig = reshape(oo.dr.derivs.dSigma_e,exo_nbr*exo_nbr,totparam_nbr); 419 indvechSig= find(tril(ones(exo_nbr,exo_nbr))); 420 dvechSig = dvechSig(indvechSig,:); 421 Inx = eye(minnx); 422 Inu = eye(exo_nbr); 423 [~,Enu] = duplication(exo_nbr); 424 KomunjerNg_DL = [dminA; dminB; dminC; dminD; dvechSig]; 425 KomunjerNg_DT = [kron(transpose(minA),Inx) - kron(Inx,minA); 426 kron(transpose(minB),Inx); 427 -1*kron(Inx,minC); 428 zeros(obs_nbr*exo_nbr,minnx^2); 429 zeros(exo_nbr*(exo_nbr+1)/2,minnx^2)]; 430 KomunjerNg_DU = [zeros(minnx^2,exo_nbr^2); 431 kron(Inu,minB); 432 zeros(obs_nbr*minnx,exo_nbr^2); 433 kron(Inu,minD); 434 -2*Enu*kron(Sigma_e,Inu)]; 435 dMINIMAL = full([KomunjerNg_DL KomunjerNg_DT KomunjerNg_DU]); 436 %add Jacobian of steady state (here we also allow for higher-order perturbation, i.e. only the mean provides additional restrictions 437 dMINIMAL = [dMEAN zeros(obs_nbr,minnx^2+exo_nbr^2); dMINIMAL]; 438 end 439 end 440else 441 dMINIMAL = []; 442end 443 444 445function [dX] = DerivABCD(X1,dX1,X2,dX2,X3,dX3,X4,dX4) 446% function [dX] = DerivABCD(X1,dX1,X2,dX2,X3,dX3,X4,dX4) 447% ------------------------------------------------------------------------- 448% Derivative of X(p)=X1(p)*X2(p)*X3(p)*X4(p) w.r.t to p 449% See Magnus and Neudecker (1999), p. 175 450% ------------------------------------------------------------------------- 451% Inputs: Matrices X1,X2,X3,X4, and the corresponding derivatives w.r.t p. 452% Output: Derivative of product of X1*X2*X3*X4 w.r.t. p 453% ========================================================================= 454nparam = size(dX1,2); 455% If one or more matrices are left out, they are set to zero 456if nargin == 4 457 X3=speye(size(X2,2)); dX3=spalloc(numel(X3),nparam,0); 458 X4=speye(size(X3,2)); dX4=spalloc(numel(X4),nparam,0); 459elseif nargin == 6 460 X4=speye(size(X3,2)); dX4=spalloc(numel(X4),nparam,0); 461end 462 463dX = kron(transpose(X4)*transpose(X3)*transpose(X2),speye(size(X1,1)))*dX1... 464 + kron(transpose(X4)*transpose(X3),X1)*dX2... 465 + kron(transpose(X4),X1*X2)*dX3... 466 + kron(speye(size(X4,2)),X1*X2*X3)*dX4; 467end %DerivABCD end 468 469end%main function end