1function DERIVS = get_perturbation_params_derivs(M, options, estim_params, oo, indpmodel, indpstderr, indpcorr, d2flag) 2% DERIVS = get_perturbation_params_derivs(M, options, estim_params, oo, indpmodel, indpstderr, indpcorr, d2flag) 3% previously getH.m in dynare 4.5 4% ------------------------------------------------------------------------- 5% Computes derivatives (with respect to parameters) of 6% (1) steady-state (ys) and covariance matrix of shocks (Sigma_e) 7% (2) dynamic model jacobians (g1, g2, g3) 8% (3) perturbation solution matrices: 9% * order==1: ghx,ghu 10% * order==2: ghx,ghu,ghxx,ghxu,ghuu,ghs2 11% * order==3: ghx,ghu,ghxx,ghxu,ghuu,ghs2,ghxxx,ghxxu,ghxuu,ghuuu,ghxss,ghuss 12% Note that the order in the parameter Jacobians is the following: 13% (1) stderr parameters (indpstderr) 14% (2) corr parameters (indpcorr) 15% (3) model parameters (indpmodel) 16% 17% ========================================================================= 18% INPUTS 19% M: [structure] storing the model information 20% options: [structure] storing the options 21% estim_params: [structure] storing the estimation information 22% oo: [structure] storing the solution results 23% indpmodel: [modparam_nbr by 1] index of selected (estimated) parameters in M.params; 24% corresponds to model parameters (no stderr and no corr) in estimated_params block 25% indpstderr: [stderrparam_nbr by 1] index of selected (estimated) standard errors, 26% i.e. for all exogenous variables where 'stderr' is given in the estimated_params block 27% indpcorr: [corrparam_nbr by 2] matrix of selected (estimated) correlations, 28% i.e. for all exogenous variables where 'corr' is given in the estimated_params block 29% d2flag: [boolean] flag to compute second-order parameter derivatives of steady state and first-order Kalman transition matrices 30% ------------------------------------------------------------------------- 31% OUTPUTS 32% DERIVS: Structure with the following fields: 33% dYss: [endo_nbr by modparam_nbr] in DR order 34% Jacobian (wrt model parameters only) of steady state, i.e. ys(order_var,:) 35% dSigma_e: [exo_nbr by exo_nbr by totparam_nbr] in declaration order 36% Jacobian (wrt to all paramters) of covariance matrix of shocks, i.e. Sigma_e 37% dg1: [endo_nbr by yy0ex0_nbr by modparam_nbr] in DR order 38% Parameter Jacobian of first derivative (wrt dynamic model variables) of dynamic model (wrt to model parameters only) 39% dg2: [endo_nbr by yy0ex0_nbr^2*modparam_nbr] in DR order 40% Parameter Jacobian of second derivative (wrt dynamic model variables) of dynamic model (wrt to model parameters only) 41% Note that instead of tensors we use matrix notation with blocks: dg2 = [dg2_dp1 dg2_dp2 ...], 42% where dg2_dpj is [endo_nbr by yy0ex0_nbr^2] and represents the derivative of g2 wrt parameter pj 43% dg3: [endo_nbr by yy0ex0_nbr^3*modparam_nbr] in DR order 44% Parameter Jacobian of third derivative (wrt dynamic model variables) of dynamic model (wrt to model parameters only) 45% Note that instead of tensors we use matrix notation with blocks: dg3 = [dg3_dp1 dg3_dp2 ...], 46% where dg3_dpj is [endo_nbr by yy0ex0_nbr^3] and represents the derivative of g3 wrt parameter pj 47% dghx: [endo_nbr by nspred by totparam_nbr] in DR order 48% Jacobian (wrt to all parameters) of first-order perturbation solution matrix ghx 49% dghu: [endo_nbr by exo_nbr by totparam_nbr] in DR order 50% Jacobian (wrt to all parameters) of first-order perturbation solution matrix ghu 51% dOm: [endo_nbr by endo_nbr by totparam_nbr] in DR order 52% Jacobian (wrt to all paramters) of Om = ghu*Sigma_e*transpose(ghu) 53% dghxx [endo_nbr by nspred*nspred by totparam_nbr] in DR order 54% Jacobian (wrt to all parameters) of second-order perturbation solution matrix ghxx 55% dghxu [endo_nbr by nspred*exo_nbr by totparam_nbr] in DR order 56% Jacobian (wrt to all parameters) of second-order perturbation solution matrix ghxu 57% dghuu [endo_nbr by exo_nbr*exo_nbr by totparam_nbr] in DR order 58% Jacobian (wrt to all parameters) of second-order perturbation solution matrix ghuu 59% dghs2 [endo_nbr by totparam_nbr] in DR order 60% Jacobian (wrt to all parameters) of second-order perturbation solution matrix ghs2 61% dghxxx [endo_nbr by nspred*nspred*nspred by totparam_nbr] in DR order 62% Jacobian (wrt to all parameters) of third-order perturbation solution matrix ghxxx 63% dghxxu [endo_nbr by nspred*nspred*exo_nbr by totparam_nbr] in DR order 64% Jacobian (wrt to all parameters) of third-order perturbation solution matrix ghxxu 65% dghxuu [endo_nbr by nspred*exo_nbr*exo_nbr by totparam_nbr] in DR order 66% Jacobian (wrt to all parameters) of third-order perturbation solution matrix ghxuu 67% dghuuu [endo_nbr by exo_nbr*exo_nbr*exo_nbr by totparam_nbr] in DR order 68% Jacobian (wrt to all parameters) of third-order perturbation solution matrix ghuuu 69% dghxss [endo_nbr by nspred by totparam_nbr] in DR order 70% Jacobian (wrt to all parameters) of third-order perturbation solution matrix ghxss 71% dghuss [endo_nbr by exo_nbr by totparam_nbr] in DR order 72% Jacobian (wrt to all parameters) of third-order perturbation solution matrix ghuss 73% if d2flag==true, we additional output: 74% d2KalmanA: [endo_nbr*endo_nbr by totparam_nbr*(totparam_nbr+1)/2] in DR order 75% Unique entries of Hessian (wrt all parameters) of Kalman transition matrix A 76% d2Om: [endo_nbr*(endo_nbr+1)/2 by totparam_nbr*(totparam_nbr+1)/2] in DR order 77% Unique entries of Hessian (wrt all parameters) of Om=ghu*Sigma_e*transpose(ghu) 78% d2Yss: [endo_nbr by modparam_nbr by modparam_nbr] in DR order 79% Unique entries of Hessian (wrt model parameters only) of steady state ys(order_var,:) 80% 81% ------------------------------------------------------------------------- 82% This function is called by 83% * dsge_likelihood.m 84% * get_identification_jacobians.m 85% ------------------------------------------------------------------------- 86% This function calls 87% * [fname,'.dynamic'] 88% * [fname,'.dynamic_params_derivs'] 89% * [fname,'.static'] 90% * [fname,'.static_params_derivs'] 91% * commutation 92% * dyn_vech 93% * dyn_unvech 94% * fjaco 95% * get_2nd_deriv (embedded) 96% * get_2nd_deriv_mat(embedded) 97% * get_all_parameters 98% * get_all_resid_2nd_derivs (embedded) 99% * get_hess_deriv (embedded) 100% * hessian_sparse 101% * sylvester3 102% * sylvester3a 103% * get_perturbation_params_derivs_numerical_objective 104% ========================================================================= 105% Copyright (C) 2019 Dynare Team 106% 107% This file is part of Dynare. 108% 109% Dynare is free software: you can redistribute it and/or modify 110% it under the terms of the GNU General Public License as published by 111% the Free Software Foundation, either version 3 of the License, or 112% (at your option) any later version. 113% 114% Dynare is distributed in the hope that it will be useful, 115% but WITHOUT ANY WARRANTY; without even the implied warranty of 116% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 117% GNU General Public License for more details. 118% 119% You should have received a copy of the GNU General Public License 120% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 121% ========================================================================= 122% Get fields from M 123Correlation_matrix = M.Correlation_matrix; 124dname = M.dname; 125dynamic_tmp_nbr = M.dynamic_tmp_nbr; 126endo_nbr = M.endo_nbr; 127exo_nbr = M.exo_nbr; 128exo_det_nbr = M.exo_det_nbr; 129fname = M.fname; 130lead_lag_incidence = M.lead_lag_incidence; 131nfwrd = M.nfwrd; 132npred = M.npred; 133nspred = M.nspred; 134nstatic = M.nstatic; 135params = M.params; 136param_nbr = M.param_nbr; 137Sigma_e = M.Sigma_e; 138 139% Get fields from options 140analytic_derivation_mode = options.analytic_derivation_mode; 141 % analytic_derivation_mode: select method to compute Jacobians, default is 0 142 % * 0: efficient sylvester equation method to compute analytical derivatives as in Ratto & Iskrev (2012) 143 % * 1: kronecker products method to compute analytical derivatives as in Iskrev (2010), only for order=1 144 % * -1: numerical two-sided finite difference method to compute numerical derivatives of all output arguments using function get_perturbation_params_derivs_numerical_objective.m 145 % * -2: numerical two-sided finite difference method to compute numerically dYss, dg1, dg2, dg3, d2Yss and d2g1, the other output arguments are computed analytically as in kronflag=0 146gstep = options.gstep; 147order = options.order; 148if isempty(options.qz_criterium) 149 % set default value for qz_criterium: if there are no unit roots one can use 1.0 150 % If they are possible, you may have have multiple unit roots and the accuracy 151 % decreases when computing the eigenvalues in lyapunov_symm. Hence, we normally use 1+1e-6 152 options = select_qz_criterium_value(options); 153end 154qz_criterium = options.qz_criterium; 155threads_BC = options.threads.kronecker.sparse_hessian_times_B_kronecker_C; 156 157% Get fields from oo 158exo_steady_state = oo.exo_steady_state; 159ghx = oo.dr.ghx; 160ghu = oo.dr.ghu; 161if order > 1 162 ghxx = oo.dr.ghxx; 163 ghxu = oo.dr.ghxu; 164 ghuu = oo.dr.ghuu; 165 ghs2 = oo.dr.ghs2; 166end 167if order > 2 168 ghxxx = oo.dr.ghxxx; 169 ghxxu = oo.dr.ghxxu; 170 ghxuu = oo.dr.ghxuu; 171 ghuuu = oo.dr.ghuuu; 172 ghxss = oo.dr.ghxss; 173 ghuss = oo.dr.ghuss; 174end 175order_var = oo.dr.order_var; 176ys = oo.dr.ys; 177 178% Some checks 179if exo_det_nbr > 0 180 error('''get_perturbation_params_derivs'': not compatible with deterministic exogenous variables, please declare as endogenous.') 181end 182if order > 1 && analytic_derivation_mode == 1 183 %analytic derivatives using Kronecker products is implemented only at first-order, at higher order we reset to analytic derivatives with sylvester equations 184 %options.analytic_derivation_mode = 0; fprintf('As order > 1, reset ''analytic_derivation_mode'' to 0\n'); 185 analytic_derivation_mode = 0; fprintf('As order > 1, reset ''analytic_derivation_mode'' to 0\n'); 186end 187 188numerical_objective_fname = str2func('get_perturbation_params_derivs_numerical_objective'); 189idx_states = nstatic+(1:nspred); %index for state variables, in DR order 190modparam_nbr = length(indpmodel); %number of selected model parameters 191stderrparam_nbr = length(indpstderr); %number of selected stderr parameters 192corrparam_nbr = size(indpcorr,1); %number of selected corr parameters 193totparam_nbr = modparam_nbr + stderrparam_nbr + corrparam_nbr; %total number of selected parameters 194[I,~] = find(lead_lag_incidence'); %I is used to select nonzero columns of the Jacobian of endogenous variables in dynamic model files 195yy0_nbr = length(ys(I)); %number of dynamic variables 196yy0ex0_nbr = yy0_nbr+exo_nbr; %number of dynamic variables + exogenous variables 197kyy0 = nonzeros(lead_lag_incidence(:,order_var)'); %index for nonzero entries in dynamic files at t-1,t,t+1 in DR order 198kyy0ex0 = [kyy0; length(kyy0)+(1:exo_nbr)']; %dynamic files include derivatives wrt exogenous variables, note that exo_det is always 0 199if order > 1 200 k2yy0ex0 = transpose(reshape(1:yy0ex0_nbr^2,yy0ex0_nbr,yy0ex0_nbr)); %index for the second dynamic derivatives, i.e. to evaluate the derivative of f wrt to yy0ex0(i) and yy0ex0(j), in DR order 201end 202if order > 2 203 k3yy0ex0 = permute(reshape(transpose(reshape(1:yy0ex0_nbr^3,yy0ex0_nbr,yy0ex0_nbr^2)),yy0ex0_nbr,yy0ex0_nbr,yy0ex0_nbr),[2 1 3]); %index for the third dynamic derivative, i.e. df/(dyyex0_i*dyyex0_j*dyyex0_k) 204end 205 206% Check for purely backward or forward looking models 207if size(lead_lag_incidence,1)<3 208 if nfwrd == 0 %purely backward models 209 klag = lead_lag_incidence(1,order_var); %indices of lagged (i.e. t-1) variables in dynamic files, columns are in DR order 210 kcurr = lead_lag_incidence(2,order_var); %indices of current (i.e. t) variables in dynamic files, columns are in DR order 211 klead = zeros(1,size(lead_lag_incidence,2)); %indices of lead (i.e. t+1) variables in dynamic files, columns are in DR order 212 elseif npred == 0 %purely forward models 213 klag = zeros(1,size(lead_lag_incidence,2)); %indices of lagged (i.e. t-1) variables in dynamic files, columns are in DR order 214 kcurr = lead_lag_incidence(1,order_var); %indices of current (i.e. t) variables in dynamic files, columns are in DR order 215 klead = lead_lag_incidence(2,order_var); %indices of lead (i.e. t+1) variables in dynamic files, columns are in DR order 216 end 217else %normal models 218 klag = lead_lag_incidence(1,order_var); %indices of lagged (i.e. t-1) variables in dynamic files, columns are in DR order 219 kcurr = lead_lag_incidence(2,order_var); %indices of current (i.e. t) variables in dynamic files, columns are in DR order 220 klead = lead_lag_incidence(3,order_var); %indices of lead (i.e. t+1) variables in dynamic files, columns are in DR order 221end 222 223if analytic_derivation_mode < 0 224 %Create auxiliary estim_params blocks if not available for numerical derivatives, estim_params_model contains only model parameters 225 estim_params_model.np = length(indpmodel); 226 estim_params_model.param_vals(:,1) = indpmodel; 227 estim_params_model.nvx = 0; estim_params_model.ncx = 0; estim_params_model.nvn = 0; estim_params_model.ncn = 0; 228 modparam1 = get_all_parameters(estim_params_model, M); %get all selected model parameters 229 if ~isempty(indpstderr) && isempty(estim_params.var_exo) %if there are stderr parameters but no estimated_params_block 230 %provide temporary necessary information for stderr parameters 231 estim_params.nvx = length(indpstderr); 232 estim_params.var_exo(:,1) = indpstderr; 233 end 234 if ~isempty(indpcorr) && isempty(estim_params.corrx) %if there are corr parameters but no estimated_params_block 235 %provide temporary necessary information for stderr parameters 236 estim_params.ncx = size(indpcorr,1); 237 estim_params.corrx(:,1:2) = indpcorr; 238 end 239 if ~isfield(estim_params,'nvn') %stderr of measurement errors not yet 240 estim_params.nvn = 0; 241 estim_params.var_endo = []; 242 end 243 if ~isfield(estim_params,'ncn') %corr of measurement errors not yet 244 estim_params.ncn = 0; 245 estim_params.corrn = []; 246 end 247 if ~isempty(indpmodel) && isempty(estim_params.param_vals) %if there are model parameters but no estimated_params_block 248 %provide temporary necessary information for model parameters 249 estim_params.np = length(indpmodel); 250 estim_params.param_vals(:,1) = indpmodel; 251 end 252 xparam1 = get_all_parameters(estim_params, M); %get all selected stderr, corr, and model parameters in estimated_params block, stderr and corr come first, then model parameters 253end 254if d2flag 255 modparam_nbr2 = modparam_nbr*(modparam_nbr+1)/2; %number of unique entries of selected model parameters only in second-order derivative matrix 256 totparam_nbr2 = totparam_nbr*(totparam_nbr+1)/2; %number of unique entries of all selected parameters in second-order derivative matrix 257 %get indices of elements in second derivatives of parameters 258 indp2tottot = reshape(1:totparam_nbr^2,totparam_nbr,totparam_nbr); 259 indp2stderrstderr = indp2tottot(1:stderrparam_nbr , 1:stderrparam_nbr); 260 indp2stderrcorr = indp2tottot(1:stderrparam_nbr , stderrparam_nbr+1:stderrparam_nbr+corrparam_nbr); 261 indp2modmod = indp2tottot(stderrparam_nbr+corrparam_nbr+1:stderrparam_nbr+corrparam_nbr+modparam_nbr , stderrparam_nbr+corrparam_nbr+1:stderrparam_nbr+corrparam_nbr+modparam_nbr); 262 if totparam_nbr ~=1 263 indp2tottot2 = dyn_vech(indp2tottot); %index of unique second-order derivatives 264 else 265 indp2tottot2 = indp2tottot; 266 end 267 if modparam_nbr ~= 1 268 indp2modmod2 = dyn_vech(indp2modmod); %get rid of cross derivatives 269 else 270 indp2modmod2 = indp2modmod; 271 end 272 %Kalman transition matrices, as in kalman_transition_matrix.m 273 KalmanA = zeros(endo_nbr,endo_nbr); 274 KalmanA(:,idx_states) = ghx; 275 KalmanB = ghu; 276end 277 278% Store some objects 279DERIVS.indpmodel = indpmodel; 280DERIVS.indpstderr = indpstderr; 281DERIVS.indpcorr = indpcorr; 282 283if analytic_derivation_mode == -1 284%% numerical two-sided finite difference method using function get_perturbation_params_derivs_numerical_objective.m (previously thet2tau.m in Dynare 4.5) for 285% Jacobian (wrt selected stderr, corr and model parameters) of 286% - dynamic model derivatives: dg1, dg2, dg3 287% - steady state (in DR order): dYss 288% - covariance matrix of shocks: dSigma_e 289% - perturbation solution matrices: dghx, dghu, dghxx, dghxu, dghuu, dghs2, dghxxx, dghxxu, dghxuu, dghuuu, dghxss, dghuss 290 291 %Parameter Jacobian of covariance matrix and solution matrices (wrt selected stderr, corr and model paramters) 292 dSig_gh = fjaco(numerical_objective_fname, xparam1, 'perturbation_solution', estim_params, M, oo, options); 293 ind_Sigma_e = (1:exo_nbr^2); 294 ind_ghx = ind_Sigma_e(end) + (1:endo_nbr*nspred); 295 ind_ghu = ind_ghx(end) + (1:endo_nbr*exo_nbr); 296 DERIVS.dSigma_e = reshape(dSig_gh(ind_Sigma_e,:),[exo_nbr exo_nbr totparam_nbr]); %in tensor notation, wrt selected parameters 297 DERIVS.dghx = reshape(dSig_gh(ind_ghx,:),[endo_nbr nspred totparam_nbr]); %in tensor notation, wrt selected parameters 298 DERIVS.dghu = reshape(dSig_gh(ind_ghu,:),[endo_nbr exo_nbr totparam_nbr]); %in tensor notation, wrt selected parameters 299 if order > 1 300 ind_ghxx = ind_ghu(end) + (1:endo_nbr*nspred^2); 301 ind_ghxu = ind_ghxx(end) + (1:endo_nbr*nspred*exo_nbr); 302 ind_ghuu = ind_ghxu(end) + (1:endo_nbr*exo_nbr*exo_nbr); 303 ind_ghs2 = ind_ghuu(end) + (1:endo_nbr); 304 DERIVS.dghxx = reshape(dSig_gh(ind_ghxx,:), [endo_nbr nspred^2 totparam_nbr]); %in tensor notation, wrt selected parameters 305 DERIVS.dghxu = reshape(dSig_gh(ind_ghxu,:), [endo_nbr nspred*exo_nbr totparam_nbr]); %in tensor notation, wrt selected parameters 306 DERIVS.dghuu = reshape(dSig_gh(ind_ghuu,:), [endo_nbr exo_nbr*exo_nbr totparam_nbr]); %in tensor notation, wrt selected parameters 307 DERIVS.dghs2 = reshape(dSig_gh(ind_ghs2,:), [endo_nbr totparam_nbr]); %in tensor notation, wrt selected parameters 308 end 309 if order > 2 310 ind_ghxxx = ind_ghs2(end) + (1:endo_nbr*nspred^3); 311 ind_ghxxu = ind_ghxxx(end) + (1:endo_nbr*nspred^2*exo_nbr); 312 ind_ghxuu = ind_ghxxu(end) + (1:endo_nbr*nspred*exo_nbr^2); 313 ind_ghuuu = ind_ghxuu(end) + (1:endo_nbr*exo_nbr^3); 314 ind_ghxss = ind_ghuuu(end) + (1:endo_nbr*nspred); 315 ind_ghuss = ind_ghxss(end) + (1:endo_nbr*exo_nbr); 316 DERIVS.dghxxx = reshape(dSig_gh(ind_ghxxx,:), [endo_nbr nspred^3 totparam_nbr]); %in tensor notation, wrt selected parameters 317 DERIVS.dghxxu = reshape(dSig_gh(ind_ghxxu,:), [endo_nbr nspred^2*exo_nbr totparam_nbr]); %in tensor notation, wrt selected parameters 318 DERIVS.dghxuu = reshape(dSig_gh(ind_ghxuu,:), [endo_nbr nspred*exo_nbr^2 totparam_nbr]); %in tensor notation, wrt selected parameters 319 DERIVS.dghuuu = reshape(dSig_gh(ind_ghuuu,:), [endo_nbr exo_nbr^3 totparam_nbr]); %in tensor notation, wrt selected parameters 320 DERIVS.dghxss = reshape(dSig_gh(ind_ghxss,:), [endo_nbr nspred totparam_nbr]); %in tensor notation, wrt selected parameters 321 DERIVS.dghuss = reshape(dSig_gh(ind_ghuss,:), [endo_nbr exo_nbr totparam_nbr]); %in tensor notation, wrt selected parameters 322 end 323 % Parameter Jacobian of Om=ghu*Sigma_e*ghu' and Correlation_matrix (wrt selected stderr, corr and model paramters) 324 DERIVS.dOm = zeros(endo_nbr,endo_nbr,totparam_nbr); %initialize in tensor notation 325 DERIVS.dCorrelation_matrix = zeros(exo_nbr,exo_nbr,totparam_nbr); %initialize in tensor notation 326 if ~isempty(indpstderr) %derivatives of ghu wrt stderr parameters are zero by construction 327 for jp=1:stderrparam_nbr 328 DERIVS.dOm(:,:,jp) = ghu*DERIVS.dSigma_e(:,:,jp)*ghu'; 329 end 330 end 331 if ~isempty(indpcorr) %derivatives of ghu wrt corr parameters are zero by construction 332 for jp=1:corrparam_nbr 333 DERIVS.dOm(:,:,stderrparam_nbr+jp) = ghu*DERIVS.dSigma_e(:,:,stderrparam_nbr+jp)*ghu'; 334 DERIVS.dCorrelation_matrix(indpcorr(jp,1),indpcorr(jp,2),stderrparam_nbr+jp) = 1; 335 DERIVS.dCorrelation_matrix(indpcorr(jp,2),indpcorr(jp,1),stderrparam_nbr+jp) = 1;%symmetry 336 end 337 end 338 if ~isempty(indpmodel) %derivatives of Sigma_e wrt model parameters are zero by construction 339 for jp=1:modparam_nbr 340 DERIVS.dOm(:,:,stderrparam_nbr+corrparam_nbr+jp) = DERIVS.dghu(:,:,stderrparam_nbr+corrparam_nbr+jp)*Sigma_e*ghu' + ghu*Sigma_e*DERIVS.dghu(:,:,stderrparam_nbr+corrparam_nbr+jp)'; 341 end 342 end 343 344 %Parameter Jacobian of dynamic model derivatives (wrt selected model parameters only) 345 dYss_g = fjaco(numerical_objective_fname, modparam1, 'dynamic_model', estim_params_model, M, oo, options); 346 ind_Yss = 1:endo_nbr; 347 ind_g1 = ind_Yss(end) + (1:endo_nbr*yy0ex0_nbr); 348 DERIVS.dYss = dYss_g(ind_Yss, :); %in tensor notation, wrt selected model parameters only 349 DERIVS.dg1 = reshape(dYss_g(ind_g1,:),[endo_nbr, yy0ex0_nbr, modparam_nbr]); %in tensor notation, wrt selected model parameters only 350 if order > 1 351 ind_g2 = ind_g1(end) + (1:endo_nbr*yy0ex0_nbr^2); 352 DERIVS.dg2 = reshape(sparse(dYss_g(ind_g2,:)),[endo_nbr, yy0ex0_nbr^2*modparam_nbr]); %blockwise in matrix notation, i.e. [dg2_dp1 dg2_dp2 ...], where dg2_dpj has dimension endo_nbr by yy0ex0_nbr^2 353 end 354 if order > 2 355 ind_g3 = ind_g2(end) + (1:endo_nbr*yy0ex0_nbr^3); 356 DERIVS.dg3 = reshape(sparse(dYss_g(ind_g3,:)),[endo_nbr, yy0ex0_nbr^3*modparam_nbr]); %blockwise in matrix notation, i.e. [dg3_dp1 dg3_dp2 ...], where dg3_dpj has dimension endo_nbr by yy0ex0_nbr^3 357 end 358 359 if d2flag 360 % Hessian (wrt paramters) of steady state and first-order solution matrices ghx and Om 361 % note that hessian_sparse.m (contrary to hessian.m) does not take symmetry into account, but focuses already on unique values 362 options.order = 1; %make sure only first order 363 d2Yss_KalmanA_Om = hessian_sparse(numerical_objective_fname, xparam1, gstep, 'Kalman_Transition', estim_params, M, oo, options); 364 options.order = order; %make sure to set back 365 ind_KalmanA = ind_Yss(end) + (1:endo_nbr^2); 366 DERIVS.d2KalmanA = d2Yss_KalmanA_Om(ind_KalmanA, indp2tottot2); %only unique elements 367 DERIVS.d2Om = d2Yss_KalmanA_Om(ind_KalmanA(end)+1:end , indp2tottot2); %only unique elements 368 DERIVS.d2Yss = zeros(endo_nbr,modparam_nbr,modparam_nbr); %initialize 369 for j = 1:endo_nbr 370 DERIVS.d2Yss(j,:,:) = dyn_unvech(full(d2Yss_KalmanA_Om(j,indp2modmod2))); %Hessian for d2Yss, but without cross derivatives 371 end 372 end 373 374 return %[END OF MAIN FUNCTION]!!!!! 375end 376 377if analytic_derivation_mode == -2 378%% Numerical two-sided finite difference method to compute parameter derivatives of steady state and dynamic model, 379% i.e. dYss, dg1, dg2, dg3 as well as d2Yss, d2g1 numerically. 380% The parameter derivatives of perturbation solution matrices are computed analytically below (analytic_derivation_mode=0) 381 if order == 3 382 [~, g1, g2, g3] = feval([fname,'.dynamic'], ys(I), exo_steady_state', params, ys, 1); 383 g3 = unfold_g3(g3, yy0ex0_nbr); 384 elseif order == 2 385 [~, g1, g2] = feval([fname,'.dynamic'], ys(I), exo_steady_state', params, ys, 1); 386 elseif order == 1 387 [~, g1] = feval([fname,'.dynamic'], ys(I), exo_steady_state', params, ys, 1); 388 end 389 390 if d2flag 391 % computation of d2Yss and d2g1 392 % note that hessian_sparse does not take symmetry into account, i.e. compare hessian_sparse.m to hessian.m, but focuses already on unique values, which are duplicated below 393 options.order = 1; %d2flag requires only first order 394 d2Yss_g1 = hessian_sparse(numerical_objective_fname, modparam1, gstep, 'dynamic_model', estim_params_model, M, oo, options); % d2flag requires only first-order 395 options.order = order; %make sure to set back the order 396 d2Yss = reshape(full(d2Yss_g1(1:endo_nbr,:)), [endo_nbr modparam_nbr modparam_nbr]); %put into tensor notation 397 for j=1:endo_nbr 398 d2Yss(j,:,:) = dyn_unvech(dyn_vech(d2Yss(j,:,:))); %add duplicate values to full hessian 399 end 400 d2g1_full = d2Yss_g1(endo_nbr+1:end,:); 401 %store only nonzero unique entries and the corresponding indices of d2g1: 402 % rows: respective derivative term 403 % 1st column: equation number of the term appearing 404 % 2nd column: column number of variable in Jacobian of the dynamic model 405 % 3rd column: number of the first parameter in derivative 406 % 4th column: number of the second parameter in derivative 407 % 5th column: value of the Hessian term 408 ind_d2g1 = find(d2g1_full); 409 d2g1 = zeros(length(ind_d2g1),5); 410 for j=1:length(ind_d2g1) 411 [i1, i2] = ind2sub(size(d2g1_full),ind_d2g1(j)); 412 [ig1, ig2] = ind2sub(size(g1),i1); 413 [ip1, ip2] = ind2sub([modparam_nbr modparam_nbr],i2); 414 d2g1(j,:) = [ig1 ig2 ip1 ip2 d2g1_full(ind_d2g1(j))]; 415 end 416 clear d2g1_full d2Yss_g1; 417 end 418 419 %Parameter Jacobian of dynamic model derivatives (wrt selected model parameters only) 420 dYss_g = fjaco(numerical_objective_fname, modparam1, 'dynamic_model', estim_params_model, M, oo, options); 421 ind_Yss = 1:endo_nbr; 422 ind_g1 = ind_Yss(end) + (1:endo_nbr*yy0ex0_nbr); 423 dYss = dYss_g(ind_Yss,:); %in tensor notation, wrt selected model parameters only 424 dg1 = reshape(dYss_g(ind_g1,:),[endo_nbr,yy0ex0_nbr,modparam_nbr]); %in tensor notation 425 if order > 1 426 ind_g2 = ind_g1(end) + (1:endo_nbr*yy0ex0_nbr^2); 427 dg2 = reshape(sparse(dYss_g(ind_g2,:)),[endo_nbr, yy0ex0_nbr^2*modparam_nbr]); %blockwise in matrix notation, i.e. [dg2_dp1 dg2_dp2 ...], where dg2_dpj has dimension endo_nbr by yy0ex0_nbr^2 428 end 429 if order > 2 430 ind_g3 = ind_g2(end) + (1:endo_nbr*yy0ex0_nbr^3); 431 dg3 = reshape(sparse(dYss_g(ind_g3,:)), [endo_nbr, yy0ex0_nbr^3*modparam_nbr]); %blockwise in matrix notation, i.e. [dg3_dp1 dg3_dp2 ...], where dg3_dpj has dimension endo_nbr by yy0ex0_nbr^3 432 end 433 clear dYss_g 434 435elseif (analytic_derivation_mode == 0 || analytic_derivation_mode == 1) 436 %% Analytical computation of Jacobian and Hessian (wrt selected model parameters) of steady state, i.e. dYss and d2Yss 437 [~, g1_static] = feval([fname,'.static'], ys, exo_steady_state', params); %g1_static is [endo_nbr by endo_nbr] first-derivative (wrt all endogenous variables) of static model equations f, i.e. df/dys, in declaration order 438 try 439 rp_static = feval([fname,'.static_params_derivs'], ys, exo_steady_state', params); %rp_static is [endo_nbr by param_nbr] first-derivative (wrt all model parameters) of static model equations f, i.e. df/dparams, in declaration order 440 catch 441 error('For analytical parameter derivatives ''static_params_derivs.m'' file is needed, this can be created by putting identification(order=%d) into your mod file.',order) 442 end 443 dys = -g1_static\rp_static; %use implicit function theorem (equation 5 of Ratto and Iskrev (2012) to compute [endo_nbr by param_nbr] first-derivative (wrt all model parameters) of steady state for all endogenous variables analytically, note that dys is in declaration order 444 d2ys = zeros(endo_nbr, param_nbr, param_nbr); %initialize in tensor notation, note that d2ys is only needed for d2flag, i.e. for g1pp 445 if d2flag 446 [~, ~, g2_static] = feval([fname,'.static'], ys, exo_steady_state', params); %g2_static is [endo_nbr by endo_nbr^2] second derivative (wrt all endogenous variables) of static model equations f, i.e. d(df/dys)/dys, in declaration order 447 if order < 3 448 [~, g1, g2, g3] = feval([fname,'.dynamic'], ys(I), exo_steady_state', params, ys, 1); %note that g3 does not contain symmetric elements 449 g3 = unfold_g3(g3, yy0ex0_nbr); %add symmetric elements to g3 450 else 451 T = NaN(sum(dynamic_tmp_nbr(1:5))); 452 T = feval([fname, '.dynamic_g4_tt'], T, ys(I), exo_steady_state', params, ys, 1); 453 g1 = feval([fname, '.dynamic_g1'], T, ys(I), exo_steady_state', params, ys, 1, false); %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 454 g2 = feval([fname, '.dynamic_g2'], T, ys(I), exo_steady_state', params, ys, 1, false); %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 455 g3 = feval([fname, '.dynamic_g3'], T, ys(I), exo_steady_state', params, ys, 1, false); %note that g3 does not contain symmetric elements 456 g4 = feval([fname, '.dynamic_g4'], T, ys(I), exo_steady_state', params, ys, 1, false); %note that g4 does not contain symmetric elements 457 g3 = unfold_g3(g3, yy0ex0_nbr); %add symmetric elements to g3, %g3 is [endo_nbr by yy0ex0_nbr^3] third-derivative (wrt all dynamic variables) of dynamic model equations, i.e. (d(df/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order 458 g4 = unfold_g4(g4, yy0ex0_nbr); %add symmetric elements to g4, %g4 is [endo_nbr by yy0ex0_nbr^4] fourth-derivative (wrt all dynamic variables) of dynamic model equations, i.e. ((d(df/dyy0ex0)/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order 459 end 460 %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 461 %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 462 %g3 is [endo_nbr by yy0ex0_nbr^3] third-derivative (wrt all dynamic variables) of dynamic model equations, i.e. (d(df/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order 463 try 464 [~, g1p_static, rpp_static] = feval([fname,'.static_params_derivs'], ys, exo_steady_state', params); 465 %g1p_static is [endo_nbr by endo_nbr by param_nbr] first derivative (wrt all model parameters) of first-derivative (wrt all endogenous variables) of static model equations f, i.e. (df/dys)/dparams, in declaration order 466 %rpp_static is [#second_order_residual_terms by 4] and contains nonzero values and corresponding indices of second derivatives (wrt all model parameters) of static model equations f, i.e. d(df/dparams)/dparams, in declaration order, where 467 % column 1 contains equation number; column 2 contains first parameter; column 3 contains second parameter; column 4 contains value of derivative 468 catch 469 error('For analytical parameter derivatives ''static_params_derivs.m'' file is needed, this can be created by putting identification(order=%d) into your mod file.',order) 470 end 471 rpp_static = get_all_resid_2nd_derivs(rpp_static, endo_nbr, param_nbr); %make full matrix out of nonzero values and corresponding indices 472 %rpp_static is [endo_nbr by param_nbr by param_nbr] second derivatives (wrt all model parameters) of static model equations, i.e. d(df/dparams)/dparams, in declaration order 473 if isempty(find(g2_static)) 474 %auxiliary expression on page 8 of Ratto and Iskrev (2012) is zero, i.e. gam = 0 475 for j = 1:param_nbr 476 %using the implicit function theorem, equation 15 on page 7 of Ratto and Iskrev (2012) 477 d2ys(:,:,j) = -g1_static\rpp_static(:,:,j); 478 %d2ys is [endo_nbr by param_nbr by param_nbr] second-derivative (wrt all model parameters) of steady state of all endogenous variables, i.e. d(dys/dparams)/dparams, in declaration order 479 end 480 else 481 gam = zeros(endo_nbr,param_nbr,param_nbr); %initialize auxiliary expression on page 8 of Ratto and Iskrev (2012) 482 for j = 1:endo_nbr 483 tmp_g1p_static_dys = (squeeze(g1p_static(j,:,:))'*dys); 484 gam(j,:,:) = transpose(reshape(g2_static(j,:),[endo_nbr endo_nbr])*dys)*dys + tmp_g1p_static_dys + tmp_g1p_static_dys'; 485 end 486 for j = 1:param_nbr 487 %using the implicit function theorem, equation 15 on page 7 of Ratto and Iskrev (2012) 488 d2ys(:,:,j) = -g1_static\(rpp_static(:,:,j)+gam(:,:,j)); 489 %d2ys is [endo_nbr by param_nbr by param_nbr] second-derivative (wrt all model parameters) of steady state of all endogenous variables, i.e. d(dys/dparams)/dparams, in declaration order 490 end 491 clear g1p_static g2_static tmp_g1p_static_dys gam 492 end 493 end 494 %handling of steady state for nonstationary variables 495 if any(any(isnan(dys))) 496 [U,T] = schur(g1_static); 497 e1 = abs(ordeig(T)) < qz_criterium-1; 498 k = sum(e1); % Number of non stationary variables. 499 % Number of stationary variables: n = length(e1)-k 500 [U,T] = ordschur(U,T,e1); 501 T = T(k+1:end,k+1:end); 502 %using implicit function theorem, equation 5 of Ratto and Iskrev (2012), in declaration order 503 dys = -U(:,k+1:end)*(T\U(:,k+1:end)')*rp_static; 504 if d2flag 505 fprintf('Computation of d2ys for nonstationary variables is not yet correctly handled if g2_static is nonempty, but continue anyways...\n') 506 for j = 1:param_nbr 507 %using implicit function theorem, equation 15 of Ratto and Iskrev (2012), in declaration order 508 d2ys(:,:,j) = -U(:,k+1:end)*(T\U(:,k+1:end)')*rpp_static(:,:,j); %THIS IS NOT CORRECT, IF g2_static IS NONEMPTY. WE NEED TO ADD GAM [willi] 509 end 510 end 511 end 512 513 if d2flag 514 try 515 if order < 3 516 [~, g1p, ~, g1pp, g2p] = feval([fname,'.dynamic_params_derivs'], ys(I), exo_steady_state', params, ys, 1, dys, d2ys); 517 else 518 [~, g1p, ~, g1pp, g2p, g3p] = feval([fname,'.dynamic_params_derivs'], ys(I), exo_steady_state', params, ys, 1, dys, d2ys); 519 end 520 catch 521 error('For analytical parameter derivatives ''dynamic_params_derivs.m'' file is needed, this can be created by putting identification(order=%d) into your mod file.',order) 522 end 523 %g1pp are nonzero values and corresponding indices of second-derivatives (wrt all model parameters) of first-derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(d(df/dyy0ex0)/dparam)/dparam, rows are in declaration order, first column in declaration order 524 d2Yss = d2ys(order_var,indpmodel,indpmodel); %[endo_nbr by mod_param_nbr by mod_param_nbr], put into DR order and focus only on selected model parameters 525 else 526 if order == 1 527 try 528 [~, g1p] = feval([fname,'.dynamic_params_derivs'], ys(I), exo_steady_state', params, ys, 1, dys, d2ys); 529 %g1p is [endo_nbr by yy0ex0_nbr by param_nbr] first-derivative (wrt all model parameters) of first-derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dparam, rows are in declaration order, column in lead_lag_incidence order 530 catch 531 error('For analytical parameter derivatives ''dynamic_params_derivs.m'' file is needed, this can be created by putting identification(order=%d) into your mod file.',order) 532 end 533 [~, g1, g2 ] = feval([fname,'.dynamic'], ys(I), exo_steady_state', params, ys, 1); 534 %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 535 %g2 is [endo_nbr by yy0ex0_nbr^2] second derivatives (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 536 elseif order == 2 537 try 538 [~, g1p, ~, ~, g2p] = feval([fname,'.dynamic_params_derivs'], ys(I), exo_steady_state', params, ys, 1, dys, d2ys); 539 %g1p is [endo_nbr by yy0ex0_nbr by param_nbr] first-derivative (wrt all model parameters) of first-derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dparam, rows are in declaration order, column in lead_lag_incidence order 540 %g2p are nonzero values and corresponding indices of first-derivative (wrt all model parameters) of second-derivatives (wrt all dynamic variables) of dynamic model equations, i.e. d(d(df/dyy0ex0)/dyy0ex0)/dparam, rows are in declaration order, first and second column in declaration order 541 catch 542 error('For analytical parameter derivatives ''dynamic_params_derivs.m'' file is needed, this can be created by putting identification(order=%d) into your mod file.',order) 543 end 544 [~, g1, g2, g3] = feval([fname,'.dynamic'], ys(I), exo_steady_state', params, ys, 1); %note that g3 does not contain symmetric elements 545 g3 = unfold_g3(g3, yy0ex0_nbr); %add symmetric elements to g3 546 %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 547 %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 548 %g3 is [endo_nbr by yy0ex0_nbr^3] third-derivative (wrt all dynamic variables) of dynamic model equations, i.e. (d(df/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order 549 elseif order == 3 550 try 551 [~, g1p, ~, ~, g2p, g3p] = feval([fname,'.dynamic_params_derivs'], ys(I), exo_steady_state', params, ys, 1, dys, d2ys); 552 %g1p is [endo_nbr by yy0ex0_nbr by param_nbr] first-derivative (wrt all model parameters) of first-derivative (wrt all dynamic variables) of dynamic model equations, i.e. d(df/dyy0ex0)/dparam, rows are in declaration order, column in lead_lag_incidence order 553 %g2p are nonzero values and corresponding indices of first-derivative (wrt all model parameters) of second-derivatives (wrt all dynamic variables) of dynamic model equations, i.e. d(d(df/dyy0ex0)/dyy0ex0)/dparam, rows are in declaration order, first and second column in declaration order 554 %g3p are nonzero values and corresponding indices of first-derivative (wrt all model parameters) of third-derivatives (wrt all dynamic variables) of dynamic model equations, i.e. d(d(d(df/dyy0ex0)/dyy0ex0)/dyy0ex0)/dparam, rows are in declaration order, first, second and third column in declaration order 555 catch 556 error('For analytical parameter derivatives ''dynamic_params_derivs.m'' file is needed, this can be created by putting identification(order=%d) into your mod file.',order) 557 end 558 T = NaN(sum(dynamic_tmp_nbr(1:5))); 559 T = feval([fname, '.dynamic_g4_tt'], T, ys(I), exo_steady_state', params, ys, 1); 560 g1 = feval([fname, '.dynamic_g1'], T, ys(I), exo_steady_state', params, ys, 1, false); %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 561 g2 = feval([fname, '.dynamic_g2'], T, ys(I), exo_steady_state', params, ys, 1, false); %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 562 g3 = feval([fname, '.dynamic_g3'], T, ys(I), exo_steady_state', params, ys, 1, false); %note that g3 does not contain symmetric elements 563 g4 = feval([fname, '.dynamic_g4'], T, ys(I), exo_steady_state', params, ys, 1, false); %note that g4 does not contain symmetric elements 564 g3 = unfold_g3(g3, yy0ex0_nbr); %add symmetric elements to g3, %g3 is [endo_nbr by yy0ex0_nbr^3] third-derivative (wrt all dynamic variables) of dynamic model equations, i.e. (d(df/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order 565 g4 = unfold_g4(g4, yy0ex0_nbr); %add symmetric elements to g4, %g4 is [endo_nbr by yy0ex0_nbr^4] fourth-derivative (wrt all dynamic variables) of dynamic model equations, i.e. ((d(df/dyy0ex0)/dyy0ex0)/dyy0ex0)/dyy0ex0, rows are in declaration order, columns in lead_lag_incidence order 566 end 567 end 568 % Parameter Jacobian of steady state in different orderings, note dys is in declaration order 569 dYss = dys(order_var, indpmodel); %in DR-order, focus only on selected model parameters 570 dyy0 = dys(I,:); %in lead_lag_incidence order, Jacobian of dynamic (without exogenous) variables, focus on all model parameters 571 dyy0ex0 = sparse([dyy0; zeros(exo_nbr,param_nbr)]); %in lead_lag_incidence order, Jacobian of dynamic (with exogenous) variables, focus on all model parameters 572 573 %% Analytical computation of Jacobian (wrt selected model parameters) of first-derivative of dynamic model, i.e. dg1 574 % Note that we use the implicit function theorem (see Ratto and Iskrev (2012) formula 7): 575 % Let g1 = df/dyy0ex0 be the first derivative (wrt all dynamic variables) of the dynamic model, then 576 % dg1 denotes the first-derivative (wrt model parameters) of g1 evaluated at the steady state. 577 % Note that g1 is a function of both the model parameters p and of the steady state of all dynamic variables, which also depend on the model parameters. 578 % Hence, implicitly g1=g1(p,yy0ex0(p)) and dg1 consists of two parts: 579 % (1) g1p: direct derivative (wrt to all model parameters) given by the preprocessor 580 % (2) g2*dyy0ex0: contribution of derivative of steady state of dynamic variables (wrt all model parameters) 581 % Note that in a stochastic context ex0 and dex0 is always zero and hence can be skipped in the computations 582 dg1 = zeros(endo_nbr,yy0ex0_nbr,param_nbr); %initialize part (2) 583 for j = 1:endo_nbr 584 [II, JJ] = ind2sub([yy0ex0_nbr yy0ex0_nbr], find(g2(j,:))); %g2 is [endo_nbr by yy0ex0_nbr^2] 585 for i = 1:yy0ex0_nbr 586 is = find(II==i); 587 is = is(find(JJ(is)<=yy0_nbr)); %focus only on oo.dr.ys(I) derivatives as exogenous variables are 0 in a stochastic context 588 if ~isempty(is) 589 tmp_g2 = full(g2(j,find(g2(j,:)))); 590 dg1(j,i,:) = tmp_g2(is)*dyy0(JJ(is),:); %put into tensor notation 591 end 592 end 593 end 594 dg1 = g1p + dg1; %add part (1) 595 dg1 = dg1(:,:,indpmodel); %focus only on selected model parameters 596 597 if order>1 598 %% Analytical computation of Jacobian (wrt selected model parameters) of second derivative of dynamic model, i.e. dg2 599 % We use the implicit function theorem: 600 % Let g2 = d2f/(dyy0ex0*dyy0ex0) denote the second derivative (wrt all dynamic variables) of the dynamic model, then 601 % dg2 denotes the first-derivative (wrt all model parameters) of g2 evaluated at the steady state. 602 % Note that g2 is a function of both the model parameters p and of the steady state of all dynamic variables, which also depend on the parameters. 603 % Hence, implicitly g2=g2(p,yy0ex0(p)) and dg2 consists of two parts: 604 % (1) g2p: direct derivative wrt to all model parameters given by the preprocessor 605 % and 606 % (2) g3*dyy0ex0: contribution of derivative of steady state of dynamic variables (wrt all model parameters) 607 % Note that we focus on selected model parameters only, i.e. indpmodel 608 % Also note that we stack the parameter derivatives blockwise instead of in tensors 609 dg2 = reshape(g3,[endo_nbr*yy0ex0_nbr^2 yy0ex0_nbr])*dyy0ex0(:,indpmodel); %part (2) 610 dg2 = reshape(dg2, [endo_nbr, yy0ex0_nbr^2*modparam_nbr]); 611 for jj = 1:size(g2p,1) 612 jpos = find(indpmodel==g2p(jj,4)); 613 if jpos~=0 614 dg2(g2p(jj,1), k2yy0ex0(g2p(jj,2),g2p(jj,3))+(jpos-1)*yy0ex0_nbr^2) = dg2(g2p(jj,1), k2yy0ex0(g2p(jj,2),g2p(jj,3))+(jpos-1)*yy0ex0_nbr^2) + g2p(jj,5); %add part (1) 615 end 616 end 617 end 618 619 if order>2 620 %% Analytical computation of Jacobian (wrt selected model parameters) of third derivative of dynamic model, i.e. dg3 621 % We use the implicit function theorem: 622 % Let g3 = d3f/(dyy0ex0*dyy0ex0*dyy0ex0) denote the third derivative (wrt all dynamic variables) of the dynamic model, then 623 % dg3 denotes the first-derivative (wrt all model parameters) of g3 evaluated at the steady state. 624 % Note that g3 is a function of both the model parameters p and of the steady state of all dynamic variables, which also depend on the parameters. 625 % Hence, implicitly g3=g3(p,yy0ex0(p)) and dg3 consists of two parts: 626 % (1) g3p: direct derivative wrt to all model parameters given by the preprocessor 627 % and 628 % (2) g4*dyy0ex0: contribution of derivative of steady state of dynamic variables (wrt all model parameters) 629 % Note that we focus on selected model parameters only, i.e. indpmodel 630 % Also note that we stack the parameter derivatives blockwise instead of in tensors 631 dg3 = reshape(g4,[endo_nbr*yy0ex0_nbr^3 yy0ex0_nbr])*dyy0ex0(:,indpmodel); %part (2) 632 dg3 = reshape(dg3, [endo_nbr, yy0ex0_nbr^3*modparam_nbr]); 633 for jj = 1:size(g3p,1) 634 jpos = find(indpmodel==g3p(jj,5)); 635 if jpos~=0 636 idyyy = unique(perms([g3p(jj,2) g3p(jj,3) g3p(jj,4)]),'rows'); %note that g3p does not contain symmetric terms, so we use the perms and unique functions 637 for k = 1:size(idyyy,1) 638 dg3(g3p(jj,1), k3yy0ex0(idyyy(k,1),idyyy(k,2),idyyy(k,3))+(jpos-1)*yy0ex0_nbr^3) = dg3(g3p(jj,1), k3yy0ex0(idyyy(k,1),idyyy(k,2),idyyy(k,3))+(jpos-1)*yy0ex0_nbr^3) + g3p(jj,6); %add part (1) 639 end 640 end 641 end 642 end 643 644 if d2flag 645 %% Analytical computation of Hessian (wrt selected model parameters) of first-derivative of dynamic model, i.e. d2g1 646 % We use the implicit function theorem (above we already computed: dg1 = g1p + g2*dyy0ex0): 647 % Accordingly, d2g1, the second-derivative (wrt model parameters), consists of five parts (ignoring transposes, see Ratto and Iskrev (2012) formula 16) 648 % (1) d(g1p)/dp = g1pp 649 % (2) d(g1p)/dyy0ex0*d(yy0ex0)/dp = g2p * dyy0ex0 650 % (3) d(g2)/dp * dyy0ex0 = g2p * dyy0ex0 651 % (4) d(g2)/dyy0ex0*d(dyy0ex0)/dp * dyy0ex0 = g3 * dyy0ex0 * dyy0ex0 652 % (5) g2 * d(dyy0ex0)/dp = g2 * d2yy0ex0 653 % Note that part 2 and 3 are equivalent besides the use of transpose (see Ratto and Iskrev (2012) formula 16) 654 d2g1_full = sparse(endo_nbr*yy0ex0_nbr, param_nbr*param_nbr); %initialize 655 g3_tmp = reshape(g3,[endo_nbr*yy0ex0_nbr*yy0ex0_nbr yy0ex0_nbr]); 656 d2g1_part4_left = sparse(endo_nbr*yy0ex0_nbr*yy0ex0_nbr,param_nbr); 657 for j = 1:param_nbr 658 %compute first two terms of part 4 659 d2g1_part4_left(:,j) = g3_tmp*dyy0ex0(:,j); 660 end 661 for j=1:endo_nbr 662 %Note that in the following we skip exogenous variables as these are 0 by construction in a stochastic setting 663 d2g1_part5 = reshape(g2(j,:), [yy0ex0_nbr yy0ex0_nbr]); 664 d2g1_part5 = d2g1_part5(:,1:yy0_nbr)*reshape(d2ys(I,:,:),[yy0_nbr,param_nbr*param_nbr]); 665 for i=1:yy0ex0_nbr 666 ind_part4 = sub2ind([endo_nbr yy0ex0_nbr yy0ex0_nbr], ones(yy0ex0_nbr,1)*j ,ones(yy0ex0_nbr,1)*i, (1:yy0ex0_nbr)'); 667 d2g1_part4 = (d2g1_part4_left(ind_part4,:))'*dyy0ex0; 668 d2g1_part2_and_part3 = (get_hess_deriv(g2p,j,i,yy0ex0_nbr,param_nbr))'*dyy0ex0; 669 d2g1_part1 = get_2nd_deriv_mat(g1pp,j,i,param_nbr); 670 d2g1_tmp = d2g1_part1 + d2g1_part2_and_part3 + d2g1_part2_and_part3' + d2g1_part4 + reshape(d2g1_part5(i,:,:),[param_nbr param_nbr]); 671 d2g1_tmp = d2g1_tmp(indpmodel,indpmodel); %focus only on model parameters 672 if any(any(d2g1_tmp)) 673 ind_d2g1_tmp = find(triu(d2g1_tmp)); 674 d2g1_full(sub2ind([endo_nbr yy0ex0_nbr],j,i), ind_d2g1_tmp) = transpose(d2g1_tmp(ind_d2g1_tmp)); 675 end 676 end 677 end 678 clear d2g1_tmp d2g1_part1 d2g1_part2_and_part3 d2g1_part4 d2g1_part4_left d2g1_part5 679 %store only nonzero entries and the corresponding indices of d2g1: 680 % 1st column: equation number of the term appearing 681 % 2nd column: column number of variable in Jacobian of the dynamic model 682 % 3rd column: number of the first parameter in derivative 683 % 4th column: number of the second parameter in derivative 684 % 5th column: value of the Hessian term 685 ind_d2g1 = find(d2g1_full); 686 d2g1 = zeros(length(ind_d2g1),5); 687 for j=1:length(ind_d2g1) 688 [i1, i2] = ind2sub(size(d2g1_full),ind_d2g1(j)); 689 [ig1, ig2] = ind2sub(size(g1),i1); 690 [ip1, ip2] = ind2sub([modparam_nbr modparam_nbr],i2); 691 d2g1(j,:) = [ig1 ig2 ip1 ip2 d2g1_full(ind_d2g1(j))]; 692 end 693 clear d2g1_full; 694 end 695end 696 697%% clear variables that are not used any more 698clear dys dyy0 dyy0ex0 d2ys 699clear rp_static rpp_static 700clear g1_static g1p_static g1p g1pp 701clear g2_static g2p tmp_g2 g3_tmp 702clear ind_d2g1 ind_d2g1_tmp ind_part4 i j i1 i2 ig1 ig2 I II JJ ip1 ip2 is 703if order == 1 704 clear g2 g3 705elseif order == 2 706 clear g3 707end 708 709%% Construct Jacobian (wrt all selected parameters) of Sigma_e and Corr_e for Gaussian innovations 710dSigma_e = zeros(exo_nbr,exo_nbr,totparam_nbr); %initialize 711dCorrelation_matrix = zeros(exo_nbr,exo_nbr,totparam_nbr); %initialize 712% Compute Jacobians wrt stderr parameters (these come first) 713% note that derivatives wrt stderr parameters are zero by construction for Corr_e 714if ~isempty(indpstderr) 715 for jp = 1:stderrparam_nbr 716 dSigma_e(indpstderr(jp),indpstderr(jp),jp) = 2*sqrt(Sigma_e(indpstderr(jp),indpstderr(jp))); 717 if isdiag(Sigma_e) == 0 % if there are correlated errors add cross derivatives 718 indotherex0 = 1:exo_nbr; 719 indotherex0(indpstderr(jp)) = []; 720 for kk = indotherex0 721 dSigma_e(indpstderr(jp), kk, jp) = Correlation_matrix(indpstderr(jp),kk)*sqrt(Sigma_e(kk,kk)); 722 dSigma_e(kk, indpstderr(jp), jp) = dSigma_e(indpstderr(jp), kk, jp); %symmetry 723 end 724 end 725 end 726end 727% Compute Jacobians wrt corr parameters (these come second) 728if ~isempty(indpcorr) 729 for jp = 1:corrparam_nbr 730 dSigma_e(indpcorr(jp,1),indpcorr(jp,2),stderrparam_nbr+jp) = sqrt(Sigma_e(indpcorr(jp,1),indpcorr(jp,1)))*sqrt(Sigma_e(indpcorr(jp,2),indpcorr(jp,2))); 731 dSigma_e(indpcorr(jp,2),indpcorr(jp,1),stderrparam_nbr+jp) = dSigma_e(indpcorr(jp,1),indpcorr(jp,2),stderrparam_nbr+jp); %symmetry 732 dCorrelation_matrix(indpcorr(jp,1),indpcorr(jp,2),stderrparam_nbr+jp) = 1; 733 dCorrelation_matrix(indpcorr(jp,2),indpcorr(jp,1),stderrparam_nbr+jp) = 1;%symmetry 734 end 735end 736% note that derivatives wrt model parameters (these come third) are zero by construction for Sigma_e and Corr_e 737 738%% Analytical computation of Jacobian (wrt selected model parameters) of first-order solution matrices ghx, ghu, and of Om=ghu*Sigma_e*ghu' 739%Notation: 740% fy_ = g1(:,nonzeros(klag)) in DR order 741% fy0 = g1(:,nonzeros(kcurr)) in DR order 742% fyp = g1(:,nonzeros(klead)) in DR order 743if analytic_derivation_mode == 1 744 % The following derivations are based on Iskrev (2010) and its online appendix A. 745 % Basic idea is to make use of the implicit function theorem. 746 % Define Kalman transition matrix KalmanA = [0 ghx 0], where the first zero spans nstatic columns, and the second zero nfwrd columns 747 % At first order we have: F = GAM0*KalmanA - GAM1*KalmanA*KalmanA - GAM2 = 0, where GAM0=fy0, GAM1=-fyp, GAM2 = -fy_ 748 % Note that F is a function of parameters p and KalmanA, which is also a function of p,therefore, F = F(p,KalmanA(p))=0, and hence, 749 % dF = Fp + dF_dKalmanA*dKalmanA = 0 or dKalmanA = - Fp/dF_dKalmanA 750 751 % Some auxiliary matrices 752 I_endo = speye(endo_nbr); 753 KalmanA = [zeros(endo_nbr,nstatic) ghx zeros(endo_nbr,nfwrd)]; 754 755 % Reshape to write first dynamic derivatives in the Magnus and Neudecker style, i.e. dvec(X)/dp 756 GAM0 = zeros(endo_nbr,endo_nbr); 757 GAM0(:,kcurr~=0,:) = g1(:,nonzeros(kcurr)); 758 dGAM0 = zeros(endo_nbr,endo_nbr,modparam_nbr); 759 dGAM0(:,kcurr~=0,:) = dg1(:,nonzeros(kcurr),:); 760 dGAM0 = reshape(dGAM0, endo_nbr*endo_nbr, modparam_nbr); 761 762 GAM1 = zeros(endo_nbr,endo_nbr); 763 GAM1(:,klead~=0,:) = -g1(:,nonzeros(klead)); 764 dGAM1 = zeros(endo_nbr,endo_nbr,modparam_nbr); 765 dGAM1(:,klead~=0,:) = -dg1(:,nonzeros(klead),:); 766 dGAM1 = reshape(dGAM1, endo_nbr*endo_nbr, modparam_nbr); 767 768 dGAM2 = zeros(endo_nbr,endo_nbr,modparam_nbr); 769 dGAM2(:,klag~=0,:) = -dg1(:,nonzeros(klag),:); 770 dGAM2 = reshape(dGAM2, endo_nbr*endo_nbr, modparam_nbr); 771 772 GAM3 = -g1(:,yy0_nbr+1:end); 773 dGAM3 = reshape(-dg1(:,yy0_nbr+1:end,:), endo_nbr*exo_nbr, modparam_nbr); 774 775 % Compute dKalmanA via implicit function 776 dF_dKalmanAghx = kron(I_endo,GAM0) - kron(KalmanA',GAM1) - kron(I_endo,GAM1*KalmanA); %equation 31 in Appendix A of Iskrev (2010) 777 Fp = kron(KalmanA',I_endo)*dGAM0 - kron( (KalmanA')^2,I_endo)*dGAM1 - dGAM2; %equation 32 in Appendix A of Iskrev (2010) 778 dKalmanA = -dF_dKalmanAghx\Fp; 779 780 % Compute dBB from expressions 33 in Iskrev (2010) Appendix A 781 MM = GAM0-GAM1*KalmanA; %this corresponds to matrix M in Ratto and Iskrev (2012, page 6) 782 invMM = MM\eye(endo_nbr); 783 dghu = - kron( (invMM*GAM3)' , invMM ) * ( dGAM0 - kron( KalmanA' , I_endo ) * dGAM1 - kron( I_endo , GAM1 ) * dKalmanA ) + kron( speye(exo_nbr), invMM ) * dGAM3; 784 785 % Add derivatives for stderr and corr parameters, which are zero by construction 786 dKalmanA = [zeros(endo_nbr^2, stderrparam_nbr+corrparam_nbr) dKalmanA]; 787 dghu = [zeros(endo_nbr*exo_nbr, stderrparam_nbr+corrparam_nbr) dghu]; 788 789 % Compute dOm = dvec(ghu*Sigma_e*ghu') from expressions 34 in Iskrev (2010) Appendix A 790 dOm = kron(I_endo,ghu*Sigma_e)*(commutation(endo_nbr, exo_nbr)*dghu)... 791 + kron(ghu,ghu)*reshape(dSigma_e, exo_nbr^2, totparam_nbr) + kron(ghu*Sigma_e,I_endo)*dghu; 792 793 % Put into tensor notation 794 dKalmanA = reshape(dKalmanA, endo_nbr, endo_nbr, totparam_nbr); 795 dghx = dKalmanA(:, nstatic+(1:nspred), stderrparam_nbr+corrparam_nbr+1:end); %get rid of zeros and focus on modparams only 796 dghu = reshape(dghu, endo_nbr, exo_nbr, totparam_nbr); 797 dghu = dghu(:,:,stderrparam_nbr+corrparam_nbr+1:end); %focus only on modparams 798 dOm = reshape(dOm, endo_nbr, endo_nbr, totparam_nbr); 799 clear dF_dKalmanAghx Fp dGAM0 dGAM1 dGAM2 dGAM3 MM invMM I_endo 800 801elseif (analytic_derivation_mode == 0 || analytic_derivation_mode == -2) 802 % Here we make use of more efficient generalized sylvester equations 803 % Notation: ghx_ = ghx(idx_states,:), ghx0 = ghx(kcurr~=0,:), ghxp = ghx(klead~=0,:) 804 % Note that at first-order we have the following expressions, which are (numerically) zero: 805 % * for ghx: g1*zx = fyp*ghxp*ghx_ + fy0*ghx0 + fy_ = A*ghx + fy_ = 0 806 % Taking the differential yields a generalized sylvester equation to get dghx: A*dghx + B*dghx*ghx_ = -dg1*zx 807 % * for ghu: g1*zu = A*ghu + fu = 0 808 % Taking the differential yields an invertible equation to get dghu: A*dghu = -(dfu + dA*ghu) 809 % INITIALIZATIONS 810 % Note that A and B are the perturbation matrices (NOT the Kalman transition matrices)! 811 A = zeros(endo_nbr,endo_nbr); 812 A(:,kcurr~=0) = g1(:,nonzeros(kcurr)); 813 A(:,idx_states) = A(:,idx_states) + g1(:,nonzeros(klead))*ghx(klead~=0,:); 814 B = zeros(endo_nbr,endo_nbr); 815 B(:,nstatic+npred+1:end) = g1(:,nonzeros(klead)); 816 zx = [eye(nspred); 817 ghx(kcurr~=0,:); 818 ghx(klead~=0,:)*ghx(idx_states,:); 819 zeros(exo_nbr,nspred)]; 820 dRHSx = zeros(endo_nbr,nspred,modparam_nbr); 821 for jp=1:modparam_nbr 822 dRHSx(:,:,jp) = -dg1(:,kyy0,jp)*zx(1:yy0_nbr,:); 823 end 824 %use iterated generalized sylvester equation to compute dghx 825 dghx = sylvester3(A,B,ghx(idx_states,:),dRHSx); 826 flag = 1; icount = 0; 827 while flag && icount < 4 828 [dghx, flag] = sylvester3a(dghx,A,B,ghx(idx_states,:),dRHSx); 829 icount = icount+1; 830 end 831 832 %Compute dOm, dghu, dA, dB 833 dOm = zeros(endo_nbr,endo_nbr,totparam_nbr); %as Om=ghu*Sigma_e*ghu', we need to use totparam_nbr, because there is also a contribution from stderr and corr parameters, which we compute after modparams 834 dghu = zeros(endo_nbr,exo_nbr,modparam_nbr); 835 dA = zeros(endo_nbr,endo_nbr,modparam_nbr); %dA is also needed at higher orders 836 dA(:,kcurr~=0,:) = dg1(:,nonzeros(kcurr),:); 837 invA = inv(A); %also needed at higher orders 838 for jp=1:modparam_nbr 839 dA(:,idx_states,jp) = dA(:,idx_states,jp) + dg1(:,nonzeros(klead),jp)*ghx(klead~=0,:) + g1(:,nonzeros(klead))*dghx(klead~=0,:,jp); 840 dghu(:,:,jp) = -invA*( dg1(:,yy0_nbr+1:end,jp) + dA(:,:,jp)*ghu); 841 dOm(:,:,stderrparam_nbr+corrparam_nbr+jp) = dghu(:,:,jp)*Sigma_e*ghu' + ghu*Sigma_e*dghu(:,:,jp)'; 842 end 843 %add stderr and corr derivatives to Om=ghu*Sigma_e*ghu' 844 if ~isempty(indpstderr) 845 for jp = 1:stderrparam_nbr 846 dOm(:,:,jp) = ghu*dSigma_e(:,:,jp)*ghu'; 847 end 848 end 849 if ~isempty(indpcorr) 850 for jp = 1:corrparam_nbr 851 dOm(:,:,stderrparam_nbr+jp) = ghu*dSigma_e(:,:,stderrparam_nbr+jp)*ghu'; 852 end 853 end 854end 855 856%% Analytical computation of Jacobian (wrt selected model parameters) of second-order solution matrices ghxx,ghxu,ghuu and ghs2 857if order > 1 858 % Notation: ghxx_ = ghxx(idx_states,:), ghxx0 = ghxx(kcurr~=0,:), ghxxp = ghxx(klead~=0,:) 859 % and similar for ghxu, ghuu and ghs2 860 % Note that at second-order we have the following expressions, which are (numerically) zero: 861 % * for ghxx: A*ghxx + B*ghxx*kron(ghx_,ghx_) + g2*kron(zx,zx) = 0 862 % Taking the differential yields a generalized sylvester equation to get dghxx: A*dghxx + B*dghxx*kron(ghx_,ghx_) = RHSxx 863 % * for ghxu: A*ghxu + B*ghxx*kron(ghx_,ghu_) + g2*kron(zx,zu) = 0 864 % Taking the differential yields an invertible equation to get dghxu: A*dghxu = RHSxu 865 % * for ghuu: A*ghuu + B*ghxx*kron(ghu_,ghu_) + g2*kron(zu,zu) = 0 866 % Taking the differential yields an invertible equation to get dghuu: A*dghuu = RHSuu 867 % * for ghs2: Ahs2*ghs2 + (gg2_{++}*kron(ghup,ghup) + fyp*ghuup)*vec(Sigma_e) = 0 868 % Taking the differential yields an invertible equation to get dghs2: S*dghs2 = RHSs2 869 % * due to certainty equivalence and zero mean shocks, we note that ghxs and ghus are zero, and thus not computed 870 dB = zeros(endo_nbr,endo_nbr,modparam_nbr); %this matrix is also needed at higher orders 871 dB(:,nstatic+npred+1:end,:) = dg1(:,nonzeros(klead),:); 872 S = A + B; %needed for dghs2, but also at higher orders 873 dS = dA + dB; 874 invS = inv(S); 875 G_x_x = kron(ghx(idx_states,:),ghx(idx_states,:)); 876 dG_x_x = zeros(size(G_x_x,1),size(G_x_x,2),modparam_nbr); 877 dzx = zeros(size(zx,1),size(zx,2),modparam_nbr); 878 dRHSghxx = zeros(endo_nbr,nspred^2,modparam_nbr); 879 for jp=1:modparam_nbr 880 dzx(:,:,jp) = [zeros(nspred,nspred); 881 dghx(kcurr~=0,:,jp); 882 dghx(klead~=0,:,jp)*ghx(idx_states,:) + ghx(klead~=0,:)*dghx(idx_states,:,jp); 883 zeros(exo_nbr,nspred)]; 884 [dRHS_1, err] = sparse_hessian_times_B_kronecker_C(dg2(:,k2yy0ex0(kyy0,kyy0)+(jp-1)*yy0ex0_nbr^2),zx(1:yy0_nbr,:),threads_BC); mexErrCheck('sparse_hessian_times_B_kronecker_C', err); 885 [dRHS_2, err] = sparse_hessian_times_B_kronecker_C(g2(:,k2yy0ex0(kyy0,kyy0)),dzx(1:yy0_nbr,:,jp),zx(1:yy0_nbr,:),threads_BC); mexErrCheck('sparse_hessian_times_B_kronecker_C', err); 886 [dRHS_3, err] = sparse_hessian_times_B_kronecker_C(g2(:,k2yy0ex0(kyy0,kyy0)),zx(1:yy0_nbr,:),dzx(1:yy0_nbr,:,jp),threads_BC); mexErrCheck('sparse_hessian_times_B_kronecker_C', err); 887 dG_x_x(:,:,jp) = kron(dghx(idx_states,:,jp),ghx(idx_states,:)) + kron(ghx(idx_states,:),dghx(idx_states,:,jp)); 888 dRHSghxx(:,:,jp) = -( (dRHS_1+dRHS_2+dRHS_3) + dA(:,:,jp)*ghxx + dB(:,:,jp)*ghxx*G_x_x + B*ghxx*dG_x_x(:,:,jp) ); 889 end 890 %use iterated generalized sylvester equation to compute dghxx 891 dghxx = sylvester3(A,B,G_x_x,dRHSghxx); 892 flag = 1; icount = 0; 893 while flag && icount < 4 894 [dghxx, flag] = sylvester3a(dghxx,A,B,G_x_x,dRHSghxx); 895 icount = icount+1; 896 end 897 zu = [zeros(nspred,exo_nbr); 898 ghu(kcurr~=0,:); 899 ghx(klead~=0,:)*ghu(idx_states,:); 900 eye(exo_nbr)]; 901 [abcOutxu,err] = A_times_B_kronecker_C(ghxx,ghx(idx_states,:),ghu(idx_states,:)); mexErrCheck('A_times_B_kronecker_C', err); %auxiliary expressions for dghxu 902 [abcOutuu, err] = A_times_B_kronecker_C(ghxx,ghu(idx_states,:)); mexErrCheck('A_times_B_kronecker_C', err); %auxiliary expressions for dghuu 903 [RHSs2, err] = sparse_hessian_times_B_kronecker_C(g2(:,k2yy0ex0(nonzeros(klead),nonzeros(klead))), ghu(klead~=0,:),threads_BC); mexErrCheck('sparse_hessian_times_B_kronecker_C', err); 904 RHSs2 = RHSs2 + g1(:,nonzeros(klead))*ghuu(klead~=0,:); 905 dzu = zeros(size(zu,1),size(zu,2),modparam_nbr); 906 dghxu = zeros(endo_nbr,nspred*exo_nbr,modparam_nbr); 907 dghuu = zeros(endo_nbr,exo_nbr*exo_nbr,modparam_nbr); 908 dghs2 = zeros(endo_nbr,totparam_nbr); %note that for modparam we ignore the contribution by dSigma_e and add it after the computations for stderr and corr 909 for jp=1:modparam_nbr 910 dzu(:,:,jp) = [zeros(nspred,exo_nbr); 911 dghu(kcurr~=0,:,jp); 912 dghx(klead~=0,:,jp)*ghu(idx_states,:) + ghx(klead~=0,:)*dghu(idx_states,:,jp); 913 zeros(exo_nbr,exo_nbr)]; 914 %compute dghxu 915 [dRHS_1, err] = sparse_hessian_times_B_kronecker_C(dg2(:,k2yy0ex0(kyy0ex0,kyy0ex0)+(jp-1)*yy0ex0_nbr^2),zx,zu,threads_BC); mexErrCheck('sparse_hessian_times_B_kronecker_C', err); 916 [dRHS_2, err] = sparse_hessian_times_B_kronecker_C(g2(:,k2yy0ex0(kyy0ex0,kyy0ex0)),dzx(:,:,jp),zu,threads_BC); mexErrCheck('sparse_hessian_times_B_kronecker_C', err); 917 [dRHS_3, err] = sparse_hessian_times_B_kronecker_C(g2(:,k2yy0ex0(kyy0ex0,kyy0ex0)),zx,dzu(:,:,jp),threads_BC); mexErrCheck('sparse_hessian_times_B_kronecker_C', err); 918 [dabcOut_1,err] = A_times_B_kronecker_C(dghxx(:,:,jp),ghx(idx_states,:),ghu(idx_states,:)); mexErrCheck('A_times_B_kronecker_C', err); 919 [dabcOut_2,err] = A_times_B_kronecker_C(ghxx,dghx(idx_states,:,jp),ghu(idx_states,:)); mexErrCheck('A_times_B_kronecker_C', err); 920 [dabcOut_3,err] = A_times_B_kronecker_C(ghxx,ghx(idx_states,:),dghu(idx_states,:,jp)); mexErrCheck('A_times_B_kronecker_C', err); 921 dghxu(:,:,jp) = -invA * ( dA(:,:,jp)*ghxu + (dRHS_1+dRHS_2+dRHS_3) + dB(:,:,jp)*abcOutxu + B*(dabcOut_1+dabcOut_2+dabcOut_3) ); 922 %compute dghuu 923 [dRHS_1, err] = sparse_hessian_times_B_kronecker_C(dg2(:,k2yy0ex0(kyy0ex0,kyy0ex0)+(jp-1)*yy0ex0_nbr^2),zu,threads_BC); mexErrCheck('sparse_hessian_times_B_kronecker_C', err); 924 [dRHS_2, err] = sparse_hessian_times_B_kronecker_C(g2(:,k2yy0ex0(kyy0ex0,kyy0ex0)),dzu(:,:,jp),zu,threads_BC); mexErrCheck('sparse_hessian_times_B_kronecker_C', err); 925 [dRHS_3, err] = sparse_hessian_times_B_kronecker_C(g2(:,k2yy0ex0(kyy0ex0,kyy0ex0)),zu,dzu(:,:,jp),threads_BC); mexErrCheck('sparse_hessian_times_B_kronecker_C', err); 926 [dabcOut_1, err] = A_times_B_kronecker_C(dghxx(:,:,jp),ghu(idx_states,:)); mexErrCheck('A_times_B_kronecker_C', err); 927 [dabcOut_2, err] = A_times_B_kronecker_C(ghxx,dghu(idx_states,:,jp),ghu(idx_states,:)); mexErrCheck('A_times_B_kronecker_C', err); 928 [dabcOut_3, err] = A_times_B_kronecker_C(ghxx,ghu(idx_states,:),dghu(idx_states,:,jp)); mexErrCheck('A_times_B_kronecker_C', err); 929 dghuu(:,:,jp) = -invA * ( dA(:,:,jp)*ghuu + (dRHS_1+dRHS_2+dRHS_3) + dB(:,:,jp)*abcOutuu + B*(dabcOut_1+dabcOut_2+dabcOut_3) ); 930 %compute dghs2 931 [dRHS_1, err] = sparse_hessian_times_B_kronecker_C(dg2(:,k2yy0ex0(nonzeros(klead),nonzeros(klead))+(jp-1)*yy0ex0_nbr^2), ghu(klead~=0,:),threads_BC); mexErrCheck('sparse_hessian_times_B_kronecker_C', err); 932 [dRHS_2, err] = sparse_hessian_times_B_kronecker_C(g2(:,k2yy0ex0(nonzeros(klead),nonzeros(klead))), dghu(klead~=0,:,jp), ghu(klead~=0,:),threads_BC); mexErrCheck('sparse_hessian_times_B_kronecker_C', err); 933 [dRHS_3, err] = sparse_hessian_times_B_kronecker_C(g2(:,k2yy0ex0(nonzeros(klead),nonzeros(klead))), ghu(klead~=0,:), dghu(klead~=0,:,jp),threads_BC); mexErrCheck('sparse_hessian_times_B_kronecker_C', err); 934 dghs2(:,stderrparam_nbr+corrparam_nbr+jp) = -invS * ( dS(:,:,jp)*ghs2 + ( (dRHS_1+dRHS_2+dRHS_3) + dg1(:,nonzeros(klead),jp)*ghuu(klead~=0,:) + g1(:,nonzeros(klead))*dghuu(klead~=0,:,jp) )*M.Sigma_e(:) ); 935 end 936 %add contributions by dSigma_e to dghs2 937 if ~isempty(indpstderr) 938 for jp = 1:stderrparam_nbr 939 dghs2(:,jp) = -invS * RHSs2*vec(dSigma_e(:,:,jp)); 940 end 941 end 942 if ~isempty(indpcorr) 943 for jp = 1:corrparam_nbr 944 dghs2(:,stderrparam_nbr+jp) = -invS * RHSs2*vec(dSigma_e(:,:,stderrparam_nbr+jp)); 945 end 946 end 947end 948 949if order > 2 950 % NOTE: The computations can be improved significantly, particularly for larger models [to do: @wmutschl] 951 % Notation: ghxxx_ = ghxxx(idx_states,:), ghxxx0 = ghxxx(kcurr~=0,:), ghxxxp = ghxxx(klead~=0,:) 952 % and similar for ghxxu, ghxuu, ghuuu, ghxss, ghuss 953 954 % Note that at third-order we have the following expressions, which are (numerically) zero, given suitable tensor-unfolding permuation matrices P: 955 % * for ghxxx: A*ghxxx + B*ghxxx*kron(ghx_,kron(ghx_,ghx_)) + g3*kron(kron(zx,zx),zx) + g2*kron(zx,zxx)*P_x_xx + fyp*ghxxp*kron(ghx_,ghxx_)*P_x_xx = 0 956 % Taking the differential yields a generalized sylvester equation to get dghxxx: A*dghxxx + B*dghxxx*kron(ghx_,kron(ghx_,ghx_)) = RHSxxx 957 % * for ghxxu: A*ghxxu + B*ghxxx*kron(ghx_,kron(ghx_,ghu_)) + gg3*kron(zx,kron(zx,zu)) + gg2*(kron(zx,zxu)*P_x_xu + kron(zxx,zu)) + fyp*ghxxp*(kron(ghx_,ghxu_)*P_x_xu + kron(ghxx_,ghu_)) = 0 958 % Taking the differential yields an invertible equation to get dghxxu: A*dghxxu = RHSxxu 959 % * for ghxuu: A*ghxuu + B*ghxxx*kron(ghx_,kron(ghu_,ghu_)) + gg3*kron(zx,kron(zu,zu)) + gg2*(kron(zxu,zu)*P_xu_u + kron(zx,zuu)) + fyp*ghxxp*(kron(ghxu_,ghu_)*Pxu_u + kron(ghx_,ghuu_)) = 0 960 % Taking the differential yields an invertible equation to get dghxuu: A*dghxuu = RHSxuu 961 % * for ghuuu: A*ghuuu + B*ghxxx*kron(ghu_,kron(ghu_,ghu_)) + gg3*kron(kron(zu,zu),zu) + gg2*kron(zu,zuu)*P_u_uu + fyp*ghxxp*kron(ghu_,ghuu_)*P_u_uu = 0 962 % Taking the differential yields an invertible equation to get dghuuu: A*dghuuu = RHSuuu 963 % * for ghxss: A*ghxss + B*ghxss*ghx_ + fyp*ghxxp*kron(ghx_,ghss_) + gg2*kron(zx,zss) + Fxupup*kron(Ix,Sigma_e(:)) = 0 964 % Taking the differential yields a generalized sylvester equation to get dghxss: A*dghxss + B*dghxss*ghx_ = RHSxss 965 % * for ghuss: A*ghuss + B*ghxss*ghu_ + gg2*kron(zu,zss) + fyp*ghxxp*kron(ghu_,ghss_) + Fuupup*kron(Iu,Sigma_e(:)) = 0 966 % Taking the differential yields an invertible equation to get dghuss: A*dghuss = RHSuss 967 % * due to certainty equivalence and Gaussian shocks, we note that ghxxs, ghxus, ghuus, and ghsss are zero and thus not computed 968 969 % permutation matrices 970 id_xxx = reshape(1:nspred^3,1,nspred,nspred,nspred); 971 id_uux = reshape(1:nspred*exo_nbr^2,1,exo_nbr,exo_nbr,nspred); 972 id_uxx = reshape(1:nspred^2*exo_nbr,1,exo_nbr,nspred,nspred); 973 id_uuu = reshape(1:exo_nbr^3,1,exo_nbr,exo_nbr,exo_nbr); 974 I_xxx = speye(nspred^3); 975 I_xxu = speye(nspred^2*exo_nbr); 976 I_xuu = speye(nspred*exo_nbr^2); 977 I_uuu = speye(exo_nbr^3); 978 P_x_xx = I_xxx(:,ipermute(id_xxx,[1,3,4,2])) + I_xxx(:,ipermute(id_xxx,[1,2,4,3])) + I_xxx(:,ipermute(id_xxx,[1,2,3,4])); 979 P_x_xu = I_xxu(:,ipermute(id_uxx,[1,2,3,4])) + I_xxu(:,ipermute(id_uxx,[1,2,4,3])); 980 P_xu_u = I_xuu(:,ipermute(id_uux,[1,2,3,4])) + I_xuu(:,ipermute(id_uux,[1,3,2,4])); 981 P_u_uu = I_uuu(:,ipermute(id_uuu,[1,3,4,2])) + I_uuu(:,ipermute(id_uuu,[1,2,4,3])) + I_uuu(:,ipermute(id_uuu,[1,2,3,4])); 982 P_uu_u = I_uuu(:,ipermute(id_uuu,[1,2,3,4])) + I_uuu(:,ipermute(id_uuu,[1,3,4,2])); 983 984 zxx = [spalloc(nspred,nspred^2,0); 985 ghxx(kcurr~=0,:); 986 ghxx(klead~=0,:)*G_x_x + ghx(klead~=0,:)*ghxx(idx_states,:); 987 spalloc(exo_nbr,nspred^2,0)]; 988 G_x_x_x = kron(ghx(idx_states,:), G_x_x); 989 G_x_xx = kron(ghx(idx_states,:),ghxx(idx_states,:)); 990 Z_x_x = kron(zx,zx); 991 Z_x_x_x = kron(zx,Z_x_x); 992 Z_x_xx = kron(zx,zxx); 993 fyp_ghxxp = sparse(g1(:,nonzeros(klead))*ghxx(klead~=0,:)); 994 B_ghxxx = B*ghxxx; 995 dzxx = zeros(size(zxx,1),size(zxx,2),modparam_nbr); 996 dfyp_ghxxp = zeros(size(fyp_ghxxp,1),size(fyp_ghxxp,2),modparam_nbr); 997 998 dRHSghxxx = zeros(endo_nbr,nspred^3,modparam_nbr); 999 for jp=1:modparam_nbr 1000 dzxx(:,:,jp) = [zeros(nspred,nspred^2); 1001 dghxx(kcurr~=0,:,jp); 1002 dghxx(klead~=0,:,jp)*G_x_x + ghxx(klead~=0,:)*dG_x_x(:,:,jp) + dghx(klead~=0,:,jp)*ghxx(idx_states,:) + ghx(klead~=0,:)*dghxx(idx_states,:,jp); 1003 zeros(exo_nbr,nspred^2)]; 1004 dG_x_x_x = kron(dghx(idx_states,:,jp),G_x_x) + kron(ghx(idx_states,:),dG_x_x(:,:,jp)); 1005 dG_x_xx = kron(dghx(idx_states,:,jp),ghxx(idx_states,:)) + kron(ghx(idx_states,:),dghxx(idx_states,:,jp)); 1006 dZ_x_x = kron(dzx(:,:,jp), zx) + kron(zx, dzx(:,:,jp)); 1007 dZ_x_x_x = kron(dzx(:,:,jp), Z_x_x) + kron(zx, dZ_x_x); 1008 dZ_x_xx = kron(dzx(:,:,jp), zxx) + kron(zx, dzxx(:,:,jp)); 1009 dfyp_ghxxp(:,:,jp) = dg1(:,nonzeros(klead),jp)*ghxx(klead~=0,:) + g1(:,nonzeros(klead))*dghxx(klead~=0,:,jp); 1010 dRHSghxxx(:,:,jp) = dA(:,:,jp)*ghxxx + dB(:,:,jp)*ghxxx*G_x_x_x + B_ghxxx*dG_x_x_x; 1011 dRHSghxxx(:,:,jp) = dRHSghxxx(:,:,jp) + dg3(:,k3yy0ex0(kyy0ex0,kyy0ex0,kyy0ex0)+(jp-1)*yy0ex0_nbr^3)*Z_x_x_x + g3(:,k3yy0ex0(kyy0ex0,kyy0ex0,kyy0ex0))*dZ_x_x_x; 1012 dRHSghxxx(:,:,jp) = dRHSghxxx(:,:,jp) + dg2(:,k2yy0ex0(kyy0ex0,kyy0ex0)+(jp-1)*yy0ex0_nbr^2)*Z_x_xx*P_x_xx + g2(:,k2yy0ex0(kyy0ex0,kyy0ex0))*dZ_x_xx*P_x_xx; 1013 dRHSghxxx(:,:,jp) = dRHSghxxx(:,:,jp) + dfyp_ghxxp(:,:,jp)*G_x_xx*P_x_xx + fyp_ghxxp*dG_x_xx*P_x_xx; 1014 end 1015 dRHSghxxx = -dRHSghxxx; 1016 %use iterated generalized sylvester equation to compute dghxxx 1017 dghxxx = sylvester3(A,B,G_x_x_x,dRHSghxxx); 1018 flag = 1; icount = 0; 1019 while flag && icount < 4 1020 [dghxxx, flag] = sylvester3a(dghxxx,A,B,G_x_x_x,dRHSghxxx); 1021 icount = icount+1; 1022 end 1023 1024 %Auxiliary expressions for dghxxu, dghxuu, dghuuu, dghxss, dghuss 1025 G_x_u = kron(ghx(idx_states,:),ghu(idx_states,:)); 1026 G_u_u = kron(ghu(idx_states,:),ghu(idx_states,:)); 1027 zxu = [zeros(nspred,nspred*exo_nbr); 1028 ghxu(kcurr~=0,:); 1029 ghxx(klead~=0,:)*G_x_u + ghx(klead~=0,:)*ghxu(idx_states,:); 1030 zeros(exo_nbr,exo_nbr*nspred)]; 1031 zuu = [zeros(nspred,exo_nbr^2); 1032 ghuu(kcurr~=0,:); 1033 ghxx(klead~=0,:)*G_u_u + ghx(klead~=0,:)*ghuu(idx_states,:); 1034 zeros(exo_nbr,exo_nbr^2)]; 1035 Z_x_u = kron(zx,zu); 1036 Z_u_u = kron(zu,zu); 1037 Z_x_xu = kron(zx,zxu); 1038 Z_xx_u = kron(zxx,zu); 1039 Z_xu_u = kron(zxu,zu); 1040 Z_x_uu = kron(zx,zuu); 1041 Z_u_uu = kron(zu,zuu); 1042 Z_x_x_u = kron(Z_x_x,zu); 1043 Z_x_u_u = kron(Z_x_u,zu); 1044 Z_u_u_u = kron(Z_u_u,zu); 1045 G_x_xu = kron(ghx(idx_states,:),ghxu(idx_states,:)); 1046 G_xx_u = kron(ghxx(idx_states,:),ghu(idx_states,:)); 1047 G_xu_u = kron(ghxu(idx_states,:),ghu(idx_states,:)); 1048 G_x_uu = kron(ghx(idx_states,:),ghuu(idx_states,:)); 1049 G_u_uu = kron(ghu(idx_states,:),ghuu(idx_states,:)); 1050 G_x_x_u = kron(G_x_x,ghu(idx_states,:)); 1051 G_x_u_u = kron(G_x_u,ghu(idx_states,:)); 1052 G_u_u_u = kron(G_u_u,ghu(idx_states,:)); 1053 aux_ZP_x_xu_Z_xx_u = Z_x_xu*P_x_xu + Z_xx_u; 1054 aux_ZP_xu_u_Z_x_uu = Z_xu_u*P_xu_u + Z_x_uu; 1055 aux_GP_x_xu_G_xx_u = G_x_xu*P_x_xu + G_xx_u; 1056 aux_GP_xu_u_G_x_uu = G_xu_u*P_xu_u + G_x_uu; 1057 dghxxu = zeros(endo_nbr,nspred^2*exo_nbr,modparam_nbr); 1058 dghxuu = zeros(endo_nbr,nspred*exo_nbr^2,modparam_nbr); 1059 dghuuu = zeros(endo_nbr,exo_nbr^3,modparam_nbr); 1060 1061 %stuff for ghxss 1062 zup = [zeros(nspred,exo_nbr); 1063 zeros(length(nonzeros(kcurr)),exo_nbr); 1064 ghu(klead~=0,:); 1065 zeros(exo_nbr,exo_nbr)]; 1066 zss = [zeros(nspred,1); 1067 ghs2(kcurr~=0,:); 1068 ghs2(klead~=0,:) + ghx(klead~=0,:)*ghs2(idx_states,:); 1069 zeros(exo_nbr,1)]; 1070 zxup = [zeros(nspred,nspred*exo_nbr); 1071 zeros(length(nonzeros(kcurr)),nspred*exo_nbr); 1072 ghxu(klead~=0,:)*kron(ghx(idx_states,:),eye(exo_nbr)); 1073 zeros(exo_nbr,nspred*exo_nbr)]; 1074 zupup = [zeros(nspred,exo_nbr^2); 1075 zeros(length(nonzeros(kcurr)),exo_nbr^2); 1076 ghuu(klead~=0,:); 1077 zeros(exo_nbr,exo_nbr^2)]; 1078 G_x_ss = kron(ghx(idx_states,:),ghs2(idx_states,:)); 1079 Z_x_ss = kron(zx,zss); 1080 Z_up_up = kron(zup,zup); 1081 Z_xup_up = kron(zxup,zup); 1082 Z_x_upup = kron(zx,zupup); 1083 Z_x_up_up = kron(zx,Z_up_up); 1084 aux_ZP_xup_up_Z_x_upup = Z_xup_up*P_xu_u + Z_x_upup; 1085 Fxupup = g3(:,k3yy0ex0(kyy0ex0,kyy0ex0,kyy0ex0))*Z_x_up_up + g2(:,k2yy0ex0(kyy0ex0,kyy0ex0))*aux_ZP_xup_up_Z_x_upup + g1(:,nonzeros(klead))*ghxuu(klead~=0,:)*kron(ghx(idx_states,:),eye(exo_nbr^2)); 1086 Ix_vecSig_e = kron(speye(nspred),Sigma_e(:)); 1087 dRHSxss = zeros(endo_nbr,nspred,totparam_nbr); 1088 1089 %stuff for ghuss 1090 zuup = [zeros(nspred,exo_nbr^2); 1091 zeros(length(nonzeros(kcurr)),exo_nbr^2); 1092 ghxu(klead~=0,:)*kron(ghu(idx_states,:),eye(exo_nbr)); 1093 zeros(exo_nbr,exo_nbr^2)]; 1094 G_u_ss = kron(ghu(idx_states,:),ghs2(idx_states,:)); 1095 Z_u_ss = kron(zu,zss); 1096 Z_u_upup = kron(zu,zupup); 1097 Z_uup_up = kron(zuup,zup); 1098 Z_u_up_up = kron(zu,Z_up_up); 1099 aux_ZP_uup_up_Z_u_upup = Z_uup_up*P_uu_u + Z_u_upup; 1100 Fuupup = g3(:,k3yy0ex0(kyy0ex0,kyy0ex0,kyy0ex0))*Z_u_up_up + g2(:,k2yy0ex0(kyy0ex0,kyy0ex0))*aux_ZP_uup_up_Z_u_upup + g1(:,nonzeros(klead))*ghxuu(klead~=0,:)*kron(ghu(idx_states,:),eye(exo_nbr^2)); 1101 Iu_vecSig_e = kron(speye(exo_nbr),Sigma_e(:)); 1102 dRHSuss = zeros(endo_nbr,exo_nbr,totparam_nbr); 1103 1104 for jp=1:modparam_nbr 1105 dG_x_u = kron(dghx(idx_states,:,jp), ghu(idx_states,:)) + kron(ghx(idx_states,:), dghu(idx_states,:,jp)); 1106 dG_u_u = kron(dghu(idx_states,:,jp), ghu(idx_states,:)) + kron(ghu(idx_states,:), dghu(idx_states,:,jp)); 1107 dzxu = [zeros(nspred,nspred*exo_nbr); 1108 dghxu(kcurr~=0,:,jp); 1109 dghxx(klead~=0,:,jp)*G_x_u + ghxx(klead~=0,:)*dG_x_u + dghx(klead~=0,:,jp)*ghxu(idx_states,:) + ghx(klead~=0,:)*dghxu(idx_states,:,jp); 1110 zeros(exo_nbr,nspred*exo_nbr)]; 1111 dzuu = [zeros(nspred,exo_nbr^2); 1112 dghuu(kcurr~=0,:,jp); 1113 dghxx(klead~=0,:,jp)*G_u_u + ghxx(klead~=0,:)*dG_u_u + dghx(klead~=0,:,jp)*ghuu(idx_states,:) + ghx(klead~=0,:)*dghuu(idx_states,:,jp); 1114 zeros(exo_nbr,exo_nbr^2)]; 1115 dG_x_xu = kron(dghx(idx_states,:,jp),ghxu(idx_states,:)) + kron(ghx(idx_states,:),dghxu(idx_states,:,jp)); 1116 dG_x_uu = kron(dghx(idx_states,:,jp),ghuu(idx_states,:)) + kron(ghx(idx_states,:),dghuu(idx_states,:,jp)); 1117 dG_u_uu = kron(dghu(idx_states,:,jp),ghuu(idx_states,:)) + kron(ghu(idx_states,:),dghuu(idx_states,:,jp)); 1118 dG_xx_u = kron(dghxx(idx_states,:,jp),ghu(idx_states,:)) + kron(ghxx(idx_states,:),dghu(idx_states,:,jp)); 1119 dG_xu_u = kron(dghxu(idx_states,:,jp),ghu(idx_states,:)) + kron(ghxu(idx_states,:),dghu(idx_states,:,jp)); 1120 dG_x_x_u = kron(G_x_x,dghu(idx_states,:,jp)) + kron(dG_x_x(:,:,jp),ghu(idx_states,:)); 1121 dG_x_u_u = kron(G_x_u,dghu(idx_states,:,jp)) + kron(dG_x_u,ghu(idx_states,:)); 1122 dG_u_u_u = kron(G_u_u,dghu(idx_states,:,jp)) + kron(dG_u_u,ghu(idx_states,:)); 1123 dZ_x_u = kron(dzx(:,:,jp),zu) + kron(zx,dzu(:,:,jp)); 1124 dZ_u_u = kron(dzu(:,:,jp),zu) + kron(zu,dzu(:,:,jp)); 1125 dZ_x_x_u = kron(dzx(:,:,jp), Z_x_u) + kron(zx, dZ_x_u); 1126 dZ_x_u_u = kron(dZ_x_u, zu) + kron(Z_x_u, dzu(:,:,jp)); 1127 dZ_u_u_u = kron(dZ_u_u, zu) + kron(Z_u_u, dzu(:,:,jp)); 1128 dZ_xx_u = kron(dzxx(:,:,jp), zu) + kron(zxx, dzu(:,:,jp)); 1129 dZ_xu_u = kron(dzxu, zu) + kron(zxu, dzu(:,:,jp)); 1130 dZ_x_xu = kron(dzx(:,:,jp), zxu) + kron(zx, dzxu); 1131 dZ_x_uu = kron(dzx(:,:,jp), zuu) + kron(zx, dzuu); 1132 dZ_u_uu = kron(dzu(:,:,jp), zuu) + kron(zu, dzuu); 1133 dB_ghxxx = dB(:,:,jp)*ghxxx + B*dghxxx(:,:,jp); 1134 %Compute dghxxu 1135 dRHS = dA(:,:,jp)*ghxxu + dB_ghxxx*G_x_x_u + B_ghxxx*dG_x_x_u; 1136 dRHS = dRHS + dg3(:,k3yy0ex0(kyy0ex0,kyy0ex0,kyy0ex0)+(jp-1)*yy0ex0_nbr^3)*Z_x_x_u + g3(:,k3yy0ex0(kyy0ex0,kyy0ex0,kyy0ex0))*dZ_x_x_u; 1137 dRHS = dRHS + dg2(:,k2yy0ex0(kyy0ex0,kyy0ex0)+(jp-1)*yy0ex0_nbr^2)*aux_ZP_x_xu_Z_xx_u + g2(:,k2yy0ex0(kyy0ex0,kyy0ex0))*( dZ_x_xu*P_x_xu + dZ_xx_u ); 1138 dRHS = dRHS + dfyp_ghxxp(:,:,jp)*aux_GP_x_xu_G_xx_u + fyp_ghxxp*( dG_x_xu*P_x_xu + dG_xx_u ); 1139 dghxxu(:,:,jp) = invA* (-dRHS); 1140 %Compute dghxuu 1141 dRHS = dA(:,:,jp)*ghxuu + dB_ghxxx*G_x_u_u + B_ghxxx*dG_x_u_u; 1142 dRHS = dRHS + dg3(:,k3yy0ex0(kyy0ex0,kyy0ex0,kyy0ex0)+(jp-1)*yy0ex0_nbr^3)*Z_x_u_u + g3(:,k3yy0ex0(kyy0ex0,kyy0ex0,kyy0ex0))*dZ_x_u_u; 1143 dRHS = dRHS + dg2(:,k2yy0ex0(kyy0ex0,kyy0ex0)+(jp-1)*yy0ex0_nbr^2)*aux_ZP_xu_u_Z_x_uu + g2(:,k2yy0ex0(kyy0ex0,kyy0ex0))*( dZ_xu_u*P_xu_u + dZ_x_uu ); 1144 dRHS = dRHS + dfyp_ghxxp(:,:,jp)*aux_GP_xu_u_G_x_uu + fyp_ghxxp*( dG_xu_u*P_xu_u + dG_x_uu ); 1145 dghxuu(:,:,jp) = invA* (-dRHS); 1146 %Compute dghuuu 1147 dRHS = dA(:,:,jp)*ghuuu + dB_ghxxx*G_u_u_u + B_ghxxx*dG_u_u_u; 1148 dRHS = dRHS + dg3(:,k3yy0ex0(kyy0ex0,kyy0ex0,kyy0ex0)+(jp-1)*yy0ex0_nbr^3)*Z_u_u_u + g3(:,k3yy0ex0(kyy0ex0,kyy0ex0,kyy0ex0))*dZ_u_u_u; 1149 dRHS = dRHS + dg2(:,k2yy0ex0(kyy0ex0,kyy0ex0)+(jp-1)*yy0ex0_nbr^2)*Z_u_uu*P_u_uu + g2(:,k2yy0ex0(kyy0ex0,kyy0ex0))*dZ_u_uu*P_u_uu; 1150 dRHS = dRHS + dfyp_ghxxp(:,:,jp)*G_u_uu*P_u_uu + fyp_ghxxp*dG_u_uu*P_u_uu; 1151 dghuuu(:,:,jp) = invA* (-dRHS); 1152 %Compute dRHSxss 1153 dzup = [zeros(nspred,exo_nbr); 1154 zeros(length(nonzeros(kcurr)),exo_nbr); 1155 dghu(klead~=0,:,jp); 1156 zeros(exo_nbr,exo_nbr)]; 1157 dzss = [zeros(nspred,1); 1158 dghs2(kcurr~=0,stderrparam_nbr+corrparam_nbr+jp); 1159 dghs2(klead~=0,stderrparam_nbr+corrparam_nbr+jp) + dghx(klead~=0,:,jp)*ghs2(idx_states,:) + ghx(klead~=0,:)*dghs2(idx_states,stderrparam_nbr+corrparam_nbr+jp); 1160 zeros(exo_nbr,1)]; 1161 dzxup = [zeros(nspred,nspred*exo_nbr); 1162 zeros(length(nonzeros(kcurr)),nspred*exo_nbr); 1163 dghxu(klead~=0,:,jp)*kron(ghx(idx_states,:),eye(exo_nbr)) + ghxu(klead~=0,:)*kron(dghx(idx_states,:,jp),eye(exo_nbr)); 1164 zeros(exo_nbr,nspred*exo_nbr)]; 1165 dzupup = [zeros(nspred,exo_nbr^2); 1166 zeros(length(nonzeros(kcurr)),exo_nbr^2); 1167 dghuu(klead~=0,:,jp); 1168 zeros(exo_nbr,exo_nbr^2)]; 1169 dG_x_ss = kron(dghx(idx_states,:,jp),ghs2(idx_states,:)) + kron(ghx(idx_states,:),dghs2(idx_states,stderrparam_nbr+corrparam_nbr+jp)); 1170 dZ_x_ss = kron(dzx(:,:,jp),zss) + kron(zx,dzss); 1171 dZ_up_up = kron(dzup,zup) + kron(zup,dzup); 1172 dZ_xup_up = kron(dzxup,zup) + kron(zxup,dzup); 1173 dZ_x_upup = kron(dzx(:,:,jp),zupup) + kron(zx,dzupup); 1174 dZ_x_up_up = kron(dzx(:,:,jp),Z_up_up) + kron(zx,dZ_up_up); 1175 daux_ZP_xup_up_Z_x_upup = dZ_xup_up*P_xu_u + dZ_x_upup; 1176 dFxupup = dg3(:,k3yy0ex0(kyy0ex0,kyy0ex0,kyy0ex0)+(jp-1)*yy0ex0_nbr^3)*Z_x_up_up + g3(:,k3yy0ex0(kyy0ex0,kyy0ex0,kyy0ex0))*dZ_x_up_up... 1177 + dg2(:,k2yy0ex0(kyy0ex0,kyy0ex0)+(jp-1)*yy0ex0_nbr^2)*aux_ZP_xup_up_Z_x_upup + g2(:,k2yy0ex0(kyy0ex0,kyy0ex0))*daux_ZP_xup_up_Z_x_upup... 1178 + dg1(:,nonzeros(klead),jp)*ghxuu(klead~=0,:)*kron(ghx(idx_states,:),eye(exo_nbr^2)) + g1(:,nonzeros(klead))*dghxuu(klead~=0,:,jp)*kron(ghx(idx_states,:),eye(exo_nbr^2)) + g1(:,nonzeros(klead))*ghxuu(klead~=0,:)*kron(dghx(idx_states,:,jp),eye(exo_nbr^2)); 1179 dRHSxss(:,:,stderrparam_nbr+corrparam_nbr+jp) = dA(:,:,jp)*ghxss + dB(:,:,jp)*ghxss*ghx(idx_states,:) + B*ghxss*dghx(idx_states,:,jp); 1180 dRHSxss(:,:,stderrparam_nbr+corrparam_nbr+jp) = dRHSxss(:,:,stderrparam_nbr+corrparam_nbr+jp) + dfyp_ghxxp(:,:,jp)*G_x_ss + fyp_ghxxp*dG_x_ss; 1181 dRHSxss(:,:,stderrparam_nbr+corrparam_nbr+jp) = dRHSxss(:,:,stderrparam_nbr+corrparam_nbr+jp) + dg2(:,k2yy0ex0(kyy0ex0,kyy0ex0)+(jp-1)*yy0ex0_nbr^2)*Z_x_ss + g2(:,k2yy0ex0(kyy0ex0,kyy0ex0))*dZ_x_ss; 1182 dRHSxss(:,:,stderrparam_nbr+corrparam_nbr+jp) = dRHSxss(:,:,stderrparam_nbr+corrparam_nbr+jp) + dFxupup*Ix_vecSig_e; %missing contribution by dSigma_e 1183 %Compute dRHSuss 1184 dzuup = [zeros(nspred,exo_nbr^2); 1185 zeros(length(nonzeros(kcurr)),exo_nbr^2); 1186 dghxu(klead~=0,:,jp)*kron(ghu(idx_states,:),eye(exo_nbr)) + ghxu(klead~=0,:)*kron(dghu(idx_states,:,jp),eye(exo_nbr)); 1187 zeros(exo_nbr,exo_nbr^2)]; 1188 dG_u_ss = kron(dghu(idx_states,:,jp),ghs2(idx_states,:)) + kron(ghu(idx_states,:),dghs2(idx_states,stderrparam_nbr+corrparam_nbr+jp)); 1189 dZ_u_ss = kron(dzu(:,:,jp),zss) + kron(zu,dzss); 1190 dZ_u_upup = kron(dzu(:,:,jp),zupup) + kron(zu,dzupup); 1191 dZ_uup_up = kron(dzuup,zup) + kron(zuup,dzup); 1192 dZ_u_up_up = kron(dzu(:,:,jp),Z_up_up) + kron(zu,dZ_up_up); 1193 daux_ZP_uup_up_Z_u_upup = dZ_uup_up*P_uu_u + dZ_u_upup; 1194 dFuupup = dg3(:,k3yy0ex0(kyy0ex0,kyy0ex0,kyy0ex0)+(jp-1)*yy0ex0_nbr^3)*Z_u_up_up + g3(:,k3yy0ex0(kyy0ex0,kyy0ex0,kyy0ex0))*dZ_u_up_up... 1195 + dg2(:,k2yy0ex0(kyy0ex0,kyy0ex0)+(jp-1)*yy0ex0_nbr^2)*aux_ZP_uup_up_Z_u_upup + g2(:,k2yy0ex0(kyy0ex0,kyy0ex0))*daux_ZP_uup_up_Z_u_upup... 1196 + dg1(:,nonzeros(klead),jp)*ghxuu(klead~=0,:)*kron(ghu(idx_states,:),eye(exo_nbr^2)) + g1(:,nonzeros(klead))*dghxuu(klead~=0,:,jp)*kron(ghu(idx_states,:),eye(exo_nbr^2)) + g1(:,nonzeros(klead))*ghxuu(klead~=0,:)*kron(dghu(idx_states,:,jp),eye(exo_nbr^2)); 1197 dRHSuss(:,:,stderrparam_nbr+corrparam_nbr+jp) = dA(:,:,jp)*ghuss + dB(:,:,jp)*ghxss*ghu(idx_states,:) + B*ghxss*dghu(idx_states,:,jp); %missing dghxss 1198 dRHSuss(:,:,stderrparam_nbr+corrparam_nbr+jp) = dRHSuss(:,:,stderrparam_nbr+corrparam_nbr+jp) + dfyp_ghxxp(:,:,jp)*G_u_ss + fyp_ghxxp*dG_u_ss; 1199 dRHSuss(:,:,stderrparam_nbr+corrparam_nbr+jp) = dRHSuss(:,:,stderrparam_nbr+corrparam_nbr+jp) + dg2(:,k2yy0ex0(kyy0ex0,kyy0ex0)+(jp-1)*yy0ex0_nbr^2)*Z_u_ss + g2(:,k2yy0ex0(kyy0ex0,kyy0ex0))*dZ_u_ss; 1200 dRHSuss(:,:,stderrparam_nbr+corrparam_nbr+jp) = dRHSuss(:,:,stderrparam_nbr+corrparam_nbr+jp) + dFuupup*Iu_vecSig_e; %contribution by dSigma_e only for stderr and corr params 1201 end 1202 %Add contribution for stderr and corr params to dRHSxss and dRHSuss 1203 if ~isempty(indpstderr) 1204 for jp = 1:stderrparam_nbr 1205 dzss = [zeros(nspred,1); 1206 dghs2(kcurr~=0,jp); 1207 dghs2(klead~=0,jp) + ghx(klead~=0,:)*dghs2(idx_states,jp); 1208 zeros(exo_nbr,1)]; 1209 dG_x_ss = kron(ghx(idx_states,:),dghs2(idx_states,jp)); 1210 dZ_x_ss = kron(zx,dzss); 1211 dRHSxss(:,:,jp) = Fxupup*kron(speye(nspred),vec(dSigma_e(:,:,jp))); 1212 dRHSxss(:,:,jp) = dRHSxss(:,:,jp) + fyp_ghxxp*dG_x_ss; 1213 dRHSxss(:,:,jp) = dRHSxss(:,:,jp) + g2(:,k2yy0ex0(kyy0ex0,kyy0ex0))*dZ_x_ss; 1214 1215 dG_u_ss = kron(ghu(idx_states,:),dghs2(idx_states,jp)); 1216 dZ_u_ss = kron(zu,dzss); 1217 dRHSuss(:,:,jp) = Fuupup*kron(speye(exo_nbr),vec(dSigma_e(:,:,jp))); 1218 dRHSuss(:,:,jp) = dRHSuss(:,:,jp) + fyp_ghxxp*dG_u_ss; 1219 dRHSuss(:,:,jp) = dRHSuss(:,:,jp) + g2(:,k2yy0ex0(kyy0ex0,kyy0ex0))*dZ_u_ss; 1220 end 1221 end 1222 if ~isempty(indpcorr) 1223 for jp = (stderrparam_nbr+1):(stderrparam_nbr+corrparam_nbr) 1224 dzss = [zeros(nspred,1); 1225 dghs2(kcurr~=0,jp); 1226 dghs2(klead~=0,jp) + ghx(klead~=0,:)*dghs2(idx_states,jp); 1227 zeros(exo_nbr,1)]; 1228 dG_x_ss = kron(ghx(idx_states,:),dghs2(idx_states,jp)); 1229 dZ_x_ss = kron(zx,dzss); 1230 dRHSxss(:,:,jp) = Fxupup*kron(speye(nspred),vec(dSigma_e(:,:,jp))); 1231 dRHSxss(:,:,jp) = dRHSxss(:,:,jp) + fyp_ghxxp*dG_x_ss; 1232 dRHSxss(:,:,jp) = dRHSxss(:,:,jp) + g2(:,k2yy0ex0(kyy0ex0,kyy0ex0))*dZ_x_ss; 1233 1234 dG_u_ss = kron(ghu(idx_states,:),dghs2(idx_states,jp)); 1235 dZ_u_ss = kron(zu,dzss); 1236 dRHSuss(:,:,jp) = Fuupup*kron(speye(exo_nbr),vec(dSigma_e(:,:,jp))); 1237 dRHSuss(:,:,jp) = dRHSuss(:,:,jp) + fyp_ghxxp*dG_u_ss; 1238 dRHSuss(:,:,jp) = dRHSuss(:,:,jp) + g2(:,k2yy0ex0(kyy0ex0,kyy0ex0))*dZ_u_ss; 1239 end 1240 end 1241 dRHSxss = -dRHSxss; 1242 %use iterated generalized sylvester equation to compute dghxss 1243 dghxss = sylvester3(A,B,ghx(idx_states,:),dRHSxss); 1244 flag = 1; icount = 0; 1245 while flag && icount < 4 1246 [dghxss, flag] = sylvester3a(dghxss,A,B,ghx(idx_states,:),dRHSxss); 1247 icount = icount+1; 1248 end 1249 %Add contribution by dghxss to dRHSuss and compute it 1250 dghuss = zeros(endo_nbr,exo_nbr,totparam_nbr); 1251 for jp = 1:totparam_nbr 1252 dRHS = dRHSuss(:,:,jp) + B*dghxss(:,:,jp)*ghu(idx_states,:); 1253 dghuss(:,:,jp) = invA* (-dRHS); 1254 end 1255end 1256 1257%% Store into structure 1258DERIVS.dg1 = dg1; 1259DERIVS.dSigma_e = dSigma_e; 1260DERIVS.dCorrelation_matrix = dCorrelation_matrix; 1261DERIVS.dYss = dYss; 1262DERIVS.dghx = cat(3,zeros(endo_nbr,nspred,stderrparam_nbr+corrparam_nbr), dghx); 1263DERIVS.dghu = cat(3,zeros(endo_nbr,exo_nbr,stderrparam_nbr+corrparam_nbr), dghu); 1264DERIVS.dOm = dOm; 1265if order > 1 1266 DERIVS.dg2 = dg2; 1267 DERIVS.dghxx = cat(3,zeros(endo_nbr,nspred^2,stderrparam_nbr+corrparam_nbr), dghxx); 1268 DERIVS.dghxu = cat(3,zeros(endo_nbr,nspred*exo_nbr,stderrparam_nbr+corrparam_nbr), dghxu); 1269 DERIVS.dghuu = cat(3,zeros(endo_nbr,exo_nbr^2,stderrparam_nbr+corrparam_nbr), dghuu); 1270 DERIVS.dghs2 = dghs2; 1271end 1272if order > 2 1273 DERIVS.dg3 = dg3; 1274 DERIVS.dghxxx = cat(3,zeros(endo_nbr,nspred^3,stderrparam_nbr+corrparam_nbr), dghxxx); 1275 DERIVS.dghxxu = cat(3,zeros(endo_nbr,nspred^2*exo_nbr,stderrparam_nbr+corrparam_nbr), dghxxu); 1276 DERIVS.dghxuu = cat(3,zeros(endo_nbr,nspred*exo_nbr^2,stderrparam_nbr+corrparam_nbr), dghxuu); 1277 DERIVS.dghuuu = cat(3,zeros(endo_nbr,exo_nbr^3,stderrparam_nbr+corrparam_nbr), dghuuu); 1278 DERIVS.dghxss = dghxss; 1279 DERIVS.dghuss = dghuss; 1280end 1281 1282%% Construct Hessian (wrt all selected parameters) of ghx, and Om=ghu*Sigma_e*ghu' 1283if d2flag 1284 % Construct Hessian (wrt all selected parameters) of Sigma_e 1285 % note that we only need to focus on (stderr x stderr), (stderr x corr), (corr x stderr) parameters, because derivatives wrt all other second-cross parameters are zero by construction 1286 d2Sigma_e = zeros(exo_nbr,exo_nbr,totparam_nbr^2); %initialize full matrix, even though we'll reduce it later to unique upper triangular values 1287 % Compute Hessian of Sigma_e wrt (stderr x stderr) parameters 1288 if ~isempty(indp2stderrstderr) 1289 for jp = 1:stderrparam_nbr 1290 for ip = 1:jp 1291 if jp == ip %same stderr parameters 1292 d2Sigma_e(indpstderr(jp),indpstderr(jp),indp2stderrstderr(ip,jp)) = 2; 1293 else %different stderr parameters 1294 if isdiag(Sigma_e) == 0 % if there are correlated errors 1295 d2Sigma_e(indpstderr(jp),indpstderr(ip),indp2stderrstderr(ip,jp)) = Correlation_matrix(indpstderr(jp),indpstderr(ip)); 1296 d2Sigma_e(indpstderr(ip),indpstderr(jp),indp2stderrstderr(ip,jp)) = Correlation_matrix(indpstderr(jp),indpstderr(ip)); %symmetry 1297 end 1298 end 1299 end 1300 end 1301 end 1302 % Compute Hessian of Sigma_e wrt (stderr x corr) parameters 1303 if ~isempty(indp2stderrcorr) 1304 for jp = 1:stderrparam_nbr 1305 for ip = 1:corrparam_nbr 1306 if indpstderr(jp) == indpcorr(ip,1) %if stderr is equal to first index of corr parameter, then derivative is equal to stderr of second index 1307 d2Sigma_e(indpstderr(jp),indpcorr(ip,2),indp2stderrcorr(jp,ip)) = sqrt(Sigma_e(indpcorr(ip,2),indpcorr(ip,2))); 1308 d2Sigma_e(indpcorr(ip,2),indpstderr(jp),indp2stderrcorr(jp,ip)) = sqrt(Sigma_e(indpcorr(ip,2),indpcorr(ip,2))); % symmetry 1309 end 1310 if indpstderr(jp) == indpcorr(ip,2) %if stderr is equal to second index of corr parameter, then derivative is equal to stderr of first index 1311 d2Sigma_e(indpstderr(jp),indpcorr(ip,1),indp2stderrcorr(jp,ip)) = sqrt(Sigma_e(indpcorr(ip,1),indpcorr(ip,1))); 1312 d2Sigma_e(indpcorr(ip,1),indpstderr(jp),indp2stderrcorr(jp,ip)) = sqrt(Sigma_e(indpcorr(ip,1),indpcorr(ip,1))); % symmetry 1313 end 1314 end 1315 end 1316 end 1317 d2Sigma_e = d2Sigma_e(:,:,indp2tottot2); %focus on upper triangular hessian values only 1318 1319 % Construct nonzero derivatives wrt to t+1, i.e. GAM1=-f_{y^+} in Villemot (2011) 1320 GAM1 = zeros(endo_nbr,endo_nbr); 1321 GAM1(:,klead~=0,:) = -g1(:,nonzeros(klead)); 1322 dGAM1 = zeros(endo_nbr,endo_nbr,modparam_nbr); 1323 dGAM1(:,klead~=0,:) = -dg1(:,nonzeros(klead),:); 1324 indind = ismember(d2g1(:,2),nonzeros(klead)); 1325 tmp = d2g1(indind,:); 1326 tmp(:,end)=-tmp(:,end); 1327 d2GAM1 = tmp; 1328 indklead = find(klead~=0); 1329 for j = 1:size(tmp,1) 1330 inxinx = (nonzeros(klead)==tmp(j,2)); 1331 d2GAM1(j,2) = indklead(inxinx); 1332 end 1333 1334 % Construct nonzero derivatives wrt to t, i.e. GAM0=f_{y^0} in Villemot (2011) 1335 GAM0 = zeros(endo_nbr,endo_nbr); 1336 GAM0(:,kcurr~=0,:) = g1(:,nonzeros(kcurr)); 1337 dGAM0 = zeros(endo_nbr,endo_nbr,modparam_nbr); 1338 dGAM0(:,kcurr~=0,:) = dg1(:,nonzeros(kcurr),:); 1339 indind = ismember(d2g1(:,2),nonzeros(kcurr)); 1340 tmp = d2g1(indind,:); 1341 d2GAM0 = tmp; 1342 indkcurr = find(kcurr~=0); 1343 for j = 1:size(tmp,1) 1344 inxinx = (nonzeros(kcurr)==tmp(j,2)); 1345 d2GAM0(j,2) = indkcurr(inxinx); 1346 end 1347 1348 % Construct nonzero derivatives wrt to t-1, i.e. GAM2=-f_{y^-} in Villemot (2011) 1349 % GAM2 = zeros(endo_nbr,endo_nbr); 1350 % GAM2(:,klag~=0) = -g1(:,nonzeros(klag)); 1351 % dGAM2 = zeros(endo_nbr,endo_nbr,modparam_nbr); 1352 % dGAM2(:,klag~=0) = -dg1(:,nonzeros(klag),:); 1353 indind = ismember(d2g1(:,2),nonzeros(klag)); 1354 tmp = d2g1(indind,:); 1355 tmp(:,end) = -tmp(:,end); 1356 d2GAM2 = tmp; 1357 indklag = find(klag~=0); 1358 for j = 1:size(tmp,1) 1359 inxinx = (nonzeros(klag)==tmp(j,2)); 1360 d2GAM2(j,2) = indklag(inxinx); 1361 end 1362 1363 % Construct nonzero derivatives wrt to u_t, i.e. GAM3=-f_{u} in Villemot (2011) 1364 % GAM3 = -g1(:,yy0_nbr+1:end); 1365 % dGAM3 = -dg1(:,yy0_nbr+1:end,:); 1366 cols_ex = yy0_nbr+(1:yy0ex0_nbr); 1367 indind = ismember(d2g1(:,2),cols_ex); 1368 tmp = d2g1(indind,:); 1369 tmp(:,end) = -tmp(:,end); 1370 d2GAM3 = tmp; 1371 for j = 1:size(tmp,1) 1372 inxinx = find(cols_ex==tmp(j,2)); 1373 d2GAM3(j,2) = inxinx; 1374 end 1375 1376 clear d2g1 tmp 1377 1378 % Compute Hessian (wrt selected params) of ghx using generalized sylvester equations, see equations 17 and 18 in Ratto and Iskrev (2012) 1379 % solves MM*d2KalmanA+N*d2KalmanA*P = QQ where d2KalmanA are second order derivatives (wrt model parameters) of KalmanA 1380 QQ = zeros(endo_nbr,endo_nbr,floor(sqrt(modparam_nbr2))); 1381 jcount=0; 1382 cumjcount=0; 1383 jinx = []; 1384 x2x=sparse(endo_nbr*endo_nbr,modparam_nbr2); 1385 dKalmanA = zeros(endo_nbr,endo_nbr,modparam_nbr); 1386 dKalmanA(:,idx_states,:) = dghx; 1387 MM = (GAM0-GAM1*KalmanA); 1388 invMM = inv(MM); 1389 for i=1:modparam_nbr 1390 for j=1:i 1391 elem1 = (get_2nd_deriv(d2GAM0,endo_nbr,endo_nbr,j,i)-get_2nd_deriv(d2GAM1,endo_nbr,endo_nbr,j,i)*KalmanA); 1392 elem1 = get_2nd_deriv(d2GAM2,endo_nbr,endo_nbr,j,i)-elem1*KalmanA; 1393 elemj0 = dGAM0(:,:,j)-dGAM1(:,:,j)*KalmanA; 1394 elemi0 = dGAM0(:,:,i)-dGAM1(:,:,i)*KalmanA; 1395 elem2 = -elemj0*dKalmanA(:,:,i)-elemi0*dKalmanA(:,:,j); 1396 elem2 = elem2 + ( dGAM1(:,:,j)*dKalmanA(:,:,i) + dGAM1(:,:,i)*dKalmanA(:,:,j) )*KalmanA; 1397 elem2 = elem2 + GAM1*( dKalmanA(:,:,i)*dKalmanA(:,:,j) + dKalmanA(:,:,j)*dKalmanA(:,:,i)); 1398 jcount=jcount+1; 1399 jinx = [jinx; [j i]]; 1400 QQ(:,:,jcount) = elem1+elem2; 1401 if jcount==floor(sqrt(modparam_nbr2)) || (j*i)==modparam_nbr^2 1402 if (j*i)==modparam_nbr^2 1403 QQ = QQ(:,:,1:jcount); 1404 end 1405 xx2=sylvester3(MM,-GAM1,KalmanA,QQ); 1406 flag=1; 1407 icount=0; 1408 while flag && icount<4 1409 [xx2, flag]=sylvester3a(xx2,MM,-GAM1,KalmanA,QQ); 1410 icount = icount + 1; 1411 end 1412 x2x(:,cumjcount+1:cumjcount+jcount)=reshape(xx2,[endo_nbr*endo_nbr jcount]); 1413 cumjcount=cumjcount+jcount; 1414 jcount = 0; 1415 jinx = []; 1416 end 1417 end 1418 end 1419 clear xx2; 1420 jcount = 0; 1421 icount = 0; 1422 cumjcount = 0; 1423 MAX_DIM_MAT = 100000000; 1424 ncol = max(1,floor(MAX_DIM_MAT/(8*endo_nbr*(endo_nbr+1)/2))); 1425 ncol = min(ncol, totparam_nbr2); 1426 d2KalmanA = sparse(endo_nbr*endo_nbr,totparam_nbr2); 1427 d2Om = sparse(endo_nbr*(endo_nbr+1)/2,totparam_nbr2); 1428 d2KalmanA_tmp = zeros(endo_nbr*endo_nbr,ncol); 1429 d2Om_tmp = zeros(endo_nbr*(endo_nbr+1)/2,ncol); 1430 tmpDir = CheckPath('tmp_derivs',dname); 1431 offset = stderrparam_nbr+corrparam_nbr; 1432 % d2B = zeros(m,n,tot_param_nbr,tot_param_nbr); 1433 for j=1:totparam_nbr 1434 for i=1:j 1435 jcount=jcount+1; 1436 if j<=offset %stderr and corr parameters 1437 y = KalmanB*d2Sigma_e(:,:,jcount)*KalmanB'; 1438 d2Om_tmp(:,jcount) = dyn_vech(y); 1439 else %model parameters 1440 jind = j-offset; 1441 iind = i-offset; 1442 if i<=offset 1443 y = dghu(:,:,jind)*dSigma_e(:,:,i)*KalmanB'+KalmanB*dSigma_e(:,:,i)*dghu(:,:,jind)'; 1444 % y(abs(y)<1.e-8)=0; 1445 d2Om_tmp(:,jcount) = dyn_vech(y); 1446 else 1447 icount=icount+1; 1448 dKalmanAj = reshape(x2x(:,icount),[endo_nbr endo_nbr]); 1449 % x = get_2nd_deriv(x2x,m,m,iind,jind);%xx2(:,:,jcount); 1450 elem1 = (get_2nd_deriv(d2GAM0,endo_nbr,endo_nbr,iind,jind)-get_2nd_deriv(d2GAM1,endo_nbr,endo_nbr,iind,jind)*KalmanA); 1451 elem1 = elem1 -( dGAM1(:,:,jind)*dKalmanA(:,:,iind) + dGAM1(:,:,iind)*dKalmanA(:,:,jind) ); 1452 elemj0 = dGAM0(:,:,jind)-dGAM1(:,:,jind)*KalmanA-GAM1*dKalmanA(:,:,jind); 1453 elemi0 = dGAM0(:,:,iind)-dGAM1(:,:,iind)*KalmanA-GAM1*dKalmanA(:,:,iind); 1454 elem0 = elemj0*dghu(:,:,iind)+elemi0*dghu(:,:,jind); 1455 y = invMM * (get_2nd_deriv(d2GAM3,endo_nbr,exo_nbr,iind,jind)-elem0-(elem1-GAM1*dKalmanAj)*KalmanB); 1456 % d2B(:,:,j+length(indexo),i+length(indexo)) = y; 1457 % d2B(:,:,i+length(indexo),j+length(indexo)) = y; 1458 y = y*Sigma_e*KalmanB'+KalmanB*Sigma_e*y'+ ... 1459 dghu(:,:,jind)*Sigma_e*dghu(:,:,iind)'+dghu(:,:,iind)*Sigma_e*dghu(:,:,jind)'; 1460 % x(abs(x)<1.e-8)=0; 1461 d2KalmanA_tmp(:,jcount) = vec(dKalmanAj); 1462 % y(abs(y)<1.e-8)=0; 1463 d2Om_tmp(:,jcount) = dyn_vech(y); 1464 end 1465 end 1466 if jcount==ncol || i*j==totparam_nbr^2 1467 d2KalmanA(:,cumjcount+1:cumjcount+jcount) = d2KalmanA_tmp(:,1:jcount); 1468 % d2KalmanA(:,:,j+length(indexo),i+length(indexo)) = x; 1469 % d2KalmanA(:,:,i+length(indexo),j+length(indexo)) = x; 1470 d2Om(:,cumjcount+1:cumjcount+jcount) = d2Om_tmp(:,1:jcount); 1471 % d2Om(:,:,j+length(indexo),i+length(indexo)) = y; 1472 % d2Om(:,:,i+length(indexo),j+length(indexo)) = y; 1473 save([tmpDir filesep 'd2KalmanA_' int2str(cumjcount+1) '_' int2str(cumjcount+jcount) '.mat'],'d2KalmanA') 1474 save([tmpDir filesep 'd2Om_' int2str(cumjcount+1) '_' int2str(cumjcount+jcount) '.mat'],'d2Om') 1475 cumjcount = cumjcount+jcount; 1476 jcount=0; 1477 % d2KalmanA = sparse(m1*m1,tot_param_nbr*(tot_param_nbr+1)/2); 1478 % d2Om = sparse(m1*(m1+1)/2,tot_param_nbr*(tot_param_nbr+1)/2); 1479 d2KalmanA_tmp = zeros(endo_nbr*endo_nbr,ncol); 1480 d2Om_tmp = zeros(endo_nbr*(endo_nbr+1)/2,ncol); 1481 end 1482 end 1483 end 1484 1485 %Store into structure 1486 DERIVS.d2Yss = d2Yss; 1487 DERIVS.d2KalmanA = d2KalmanA; 1488 DERIVS.d2Om = d2Om; 1489end 1490 1491return 1492 1493%% AUXILIARY FUNCTIONS %% 1494%%%%%%%%%%%%%%%%%%%%%%%%% 1495 1496function g22 = get_2nd_deriv(gpp,m,n,i,j) 1497% inputs: 1498% - gpp: [#second_order_Jacobian_terms by 5] double Hessian matrix (wrt parameters) of a matrix 1499% rows: respective derivative term 1500% 1st column: equation number of the term appearing 1501% 2nd column: column number of variable in Jacobian 1502% 3rd column: number of the first parameter in derivative 1503% 4th column: number of the second parameter in derivative 1504% 5th column: value of the Hessian term 1505% - m: scalar number of equations 1506% - n: scalar number of variables 1507% - i: scalar number for which first parameter 1508% - j: scalar number for which second parameter 1509 1510g22=zeros(m,n); 1511is=find(gpp(:,3)==i); 1512is=is(find(gpp(is,4)==j)); 1513 1514if ~isempty(is) 1515 g22(sub2ind([m,n],gpp(is,1),gpp(is,2)))=gpp(is,5)'; 1516end 1517return 1518 1519function g22 = get_2nd_deriv_mat(gpp,i,j,npar) 1520% inputs: 1521% - gpp: [#second_order_Jacobian_terms by 5] double Hessian matrix of (wrt parameters) of dynamic Jacobian 1522% rows: respective derivative term 1523% 1st column: equation number of the term appearing 1524% 2nd column: column number of variable in Jacobian of the dynamic model 1525% 3rd column: number of the first parameter in derivative 1526% 4th column: number of the second parameter in derivative 1527% 5th column: value of the Hessian term 1528% - i: scalar number for which model equation 1529% - j: scalar number for which variable in Jacobian of dynamic model 1530% - npar: scalar Number of model parameters, i.e. equals param_nbr 1531% 1532% output: 1533% g22: [npar by npar] Hessian matrix (wrt parameters) of Jacobian of dynamic model for equation i 1534% rows: first parameter in Hessian 1535% columns: second paramater in Hessian 1536 1537g22=zeros(npar,npar); 1538is=find(gpp(:,1)==i); 1539is=is(find(gpp(is,2)==j)); 1540 1541if ~isempty(is) 1542 g22(sub2ind([npar,npar],gpp(is,3),gpp(is,4)))=gpp(is,5)'; 1543end 1544return 1545 1546function r22 = get_all_resid_2nd_derivs(rpp,m,npar) 1547% inputs: 1548% - rpp: [#second_order_residual_terms by 4] double Hessian matrix (wrt paramters) of model equations 1549% rows: respective derivative term 1550% 1st column: equation number of the term appearing 1551% 2nd column: number of the first parameter in derivative 1552% 3rd column: number of the second parameter in derivative 1553% 4th column: value of the Hessian term 1554% - m: scalar Number of residuals (or model equations), i.e. equals endo_nbr 1555% - npar: scalar Number of model parameters, i.e. equals param_nbr 1556% 1557% output: 1558% r22: [endo_nbr by param_nbr by param_nbr] Hessian matrix of model equations with respect to parameters 1559% rows: equations in order of declaration 1560% 1st columns: first parameter number in derivative 1561% 2nd columns: second parameter in derivative 1562 1563r22=zeros(m,npar,npar); 1564 1565for is=1:size(rpp,1) 1566 % Keep symmetry in hessian, hence 2 and 3 as well as 3 and 2, i.e. d2f/(dp1 dp2) = d2f/(dp2 dp1) 1567 r22(rpp(is,1),rpp(is,2),rpp(is,3))=rpp(is,4); 1568 r22(rpp(is,1),rpp(is,3),rpp(is,2))=rpp(is,4); 1569end 1570 1571return 1572 1573function h2 = get_hess_deriv(hp,i,j,m,npar) 1574% inputs: 1575% - hp: [#first_order_Hessian_terms by 5] double Jacobian matrix (wrt paramters) of dynamic Hessian 1576% rows: respective derivative term 1577% 1st column: equation number of the term appearing 1578% 2nd column: column number of first variable in Hessian of the dynamic model 1579% 3rd column: column number of second variable in Hessian of the dynamic model 1580% 4th column: number of the parameter in derivative 1581% 5th column: value of the Hessian term 1582% - i: scalar number for which model equation 1583% - j: scalar number for which first variable in Hessian of dynamic model variable 1584% - m: scalar Number of dynamic model variables + exogenous vars, i.e. yy0_nbr + exo_nbr 1585% - npar: scalar Number of model parameters, i.e. equals param_nbr 1586% 1587% output: 1588% h2: [(yy0_nbr + exo_nbr) by param_nbr] Jacobian matrix (wrt parameters) of dynamic Hessian 1589% rows: second dynamic or exogenous variables in Hessian of specific model equation of the dynamic model 1590% columns: parameters 1591 1592h2=zeros(m,npar); 1593is1=find(hp(:,1)==i); 1594is=is1(find(hp(is1,2)==j)); 1595 1596if ~isempty(is) 1597 h2(sub2ind([m,npar],hp(is,3),hp(is,4)))=hp(is,5)'; 1598end 1599 1600return 1601