960 lines
54 KiB
Matlab
960 lines
54 KiB
Matlab
function [dA, dB, dSigma_e, dOm, dYss, dg1, d2A, d2Om, d2Yss] = get_first_order_solution_params_deriv(A, B, estim_params, M, oo, options, kronflag, indpmodel, indpstderr, indpcorr, indvar)
|
|
%[dA, dB, dSigma_e, dOm, dYss, dg1, d2A, d2Om, d2Yss] = get_first_order_solution_params_deriv(A, B, estim_params, M, oo, options, kronflag, indpmodel, indpstderr, indpcorr, indvar)
|
|
% previously getH.m
|
|
% -------------------------------------------------------------------------
|
|
% Computes first and second derivatives (with respect to parameters) of
|
|
% (1) reduced-form solution (dA, dB, dSigma_e, dOm, d2A, d2Om)
|
|
% (1) steady-state (dYss, d2Yss)
|
|
% (3) Jacobian (wrt to dynamic variables) of dynamic model (dg1)
|
|
% Note that the order in the parameter Jacobians is the following:
|
|
% first stderr parameters, second corr parameters, third model parameters
|
|
% =========================================================================
|
|
% INPUTS
|
|
% A: [endo_nbr by endo_nbr] Transition matrix from Kalman filter
|
|
% for all endogenous declared variables, in DR order
|
|
% B: [endo_nbr by exo_nbr] Transition matrix from Kalman filter
|
|
% mapping shocks today to endogenous variables today, in DR order
|
|
% estim_params: [structure] storing the estimation information
|
|
% M: [structure] storing the model information
|
|
% oo: [structure] storing the reduced-form solution results
|
|
% options: [structure] storing the options
|
|
% kronflag: [scalar] method to compute Jacobians (equal to analytic_derivation_mode in options_ident). Default:0
|
|
% * 0: efficient sylvester equation method to compute
|
|
% analytical derivatives as in Ratto & Iskrev (2011)
|
|
% * 1: kronecker products method to compute analytical
|
|
% derivatives as in Iskrev (2010)
|
|
% * -1: numerical two-sided finite difference method to
|
|
% compute numerical derivatives of all output arguments
|
|
% using function identification_numerical_objective.m
|
|
% (previously thet2tau.m)
|
|
% * -2: numerical two-sided finite difference method to
|
|
% compute numerically dYss, dg1, d2Yss and d2g1, the other
|
|
% output arguments are computed analytically as in kronflag=0
|
|
% indpmodel: [modparam_nbr by 1] index of estimated parameters in M_.params;
|
|
% corresponds to model parameters (no stderr and no corr)
|
|
% in estimated_params block; if estimated_params block is
|
|
% not available, then all model parameters are selected
|
|
% indpstderr: [stderrparam_nbr by 1] index of estimated standard errors,
|
|
% i.e. for all exogenous variables where "stderr" is given
|
|
% in the estimated_params block; if estimated_params block
|
|
% is not available, then all stderr parameters are selected
|
|
% indpcorr: [corrparam_nbr by 2] matrix of estimated correlations,
|
|
% i.e. for all exogenous variables where "corr" is given
|
|
% in the estimated_params block; if estimated_params block
|
|
% is not available, then no corr parameters are selected
|
|
% indvar: [var_nbr by 1] index of considered (or observed) variables
|
|
% -------------------------------------------------------------------------
|
|
% OUTPUTS
|
|
% dA: [var_nbr by var_nbr by totparam_nbr] in DR order
|
|
% Jacobian (wrt to all parameters) of transition matrix A
|
|
% dB: [var_nbr by exo_nbr by totparam_nbr] in DR order
|
|
% Jacobian (wrt to all parameters) of transition matrix B
|
|
% dSigma_e: [exo_nbr by exo_nbr by totparam_nbr] in declaration order
|
|
% Jacobian (wrt to all paramters) of M_.Sigma_e
|
|
% dOm: [var_nbr by var_nbr by totparam_nbr] in DR order
|
|
% Jacobian (wrt to all paramters) of Om = (B*M_.Sigma_e*B')
|
|
% dYss: [var_nbr by modparam_nbr] in DR order
|
|
% Jacobian (wrt model parameters only) of steady state
|
|
% dg1: [endo_nbr by (dynamicvar_nbr + exo_nbr) by modparam_nbr] in DR order
|
|
% Jacobian (wrt to model parameters only) of Jacobian of dynamic model
|
|
% d2A: [var_nbr*var_nbr by totparam_nbr*(totparam_nbr+1)/2] in DR order
|
|
% Unique entries of Hessian (wrt all parameters) of transition matrix A
|
|
% d2Om: [var_nbr*(var_nbr+1)/2 by totparam_nbr*(totparam_nbr+1)/2] in DR order
|
|
% Unique entries of Hessian (wrt all parameters) of Omega
|
|
% d2Yss: [var_nbr by modparam_nbr by modparam_nbr] in DR order
|
|
% Unique entries of Hessian (wrt model parameters only) of steady state
|
|
% -------------------------------------------------------------------------
|
|
% This function is called by
|
|
% * dsge_likelihood.m
|
|
% * get_identification_jacobians.m (previously getJJ.m)
|
|
% -------------------------------------------------------------------------
|
|
% This function calls
|
|
% * [fname,'.dynamic']
|
|
% * [fname,'.dynamic_params_derivs']
|
|
% * [fname,'.static']
|
|
% * [fname,'.static_params_derivs']
|
|
% * commutation
|
|
% * dyn_vech
|
|
% * dyn_unvech
|
|
% * fjaco
|
|
% * get_2nd_deriv (embedded)
|
|
% * get_2nd_deriv_mat(embedded)
|
|
% * get_all_parameters
|
|
% * get_all_resid_2nd_derivs (embedded)
|
|
% * get_hess_deriv (embedded)
|
|
% * hessian_sparse
|
|
% * sylvester3
|
|
% * sylvester3a
|
|
% * identification_numerical_objective.m (previously thet2tau.m)
|
|
% =========================================================================
|
|
% Copyright (C) 2010-2019 Dynare Team
|
|
%
|
|
% This file is part of Dynare.
|
|
%
|
|
% Dynare is free software: you can redistribute it and/or modify
|
|
% it under the terms of the GNU General Public License as published by
|
|
% the Free Software Foundation, either version 3 of the License, or
|
|
% (at your option) any later version.
|
|
%
|
|
% Dynare is distributed in the hope that it will be useful,
|
|
% but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
% GNU General Public License for more details.
|
|
%
|
|
% You should have received a copy of the GNU General Public License
|
|
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
|
% =========================================================================
|
|
fname = M.fname;
|
|
dname = M.dname;
|
|
maximum_exo_lag = M.maximum_exo_lag;
|
|
maximum_exo_lead = M.maximum_exo_lead;
|
|
maximum_endo_lag = M.maximum_endo_lag;
|
|
maximum_endo_lead = M.maximum_endo_lead;
|
|
lead_lag_incidence = M.lead_lag_incidence;
|
|
[I,~] = find(lead_lag_incidence'); %I is used to select nonzero columns of the Jacobian of endogenous variables in dynamic model files
|
|
|
|
ys = oo.dr.ys; %steady state of endogenous variables in declaration order
|
|
yy0 = oo.dr.ys(I); %steady state of dynamic (endogenous and auxiliary variables) in DR order
|
|
ex0 = oo.exo_steady_state'; %steady state of exogenous variables in declaration order
|
|
params0 = M.params; %values at which to evaluate dynamic, static and param_derivs files
|
|
Sigma_e0 = M.Sigma_e; %covariance matrix of exogenous shocks
|
|
Corr_e0 = M.Correlation_matrix; %correlation matrix of exogenous shocks
|
|
stderr_e0 = sqrt(diag(Sigma_e0)); %standard errors of exogenous shocks
|
|
|
|
param_nbr = M.param_nbr; %number of all declared model parameters in mod file
|
|
modparam_nbr = length(indpmodel); %number of model parameters to be used
|
|
stderrparam_nbr = length(indpstderr); %number of stderr parameters to be used
|
|
corrparam_nbr = size(indpcorr,1); %number of stderr parameters to be used
|
|
totparam_nbr = modparam_nbr + stderrparam_nbr + corrparam_nbr; %total number of parameters to be used
|
|
|
|
if nargout > 6
|
|
modparam_nbr2 = modparam_nbr*(modparam_nbr+1)/2; %number of unique entries of model parameters only in second-order derivative matrix
|
|
totparam_nbr2 = totparam_nbr*(totparam_nbr+1)/2; %number of unique entries of all parameters in second-order derivative matrix
|
|
%get indices of elements in second derivatives of parameters
|
|
indp2tottot = reshape(1:totparam_nbr^2,totparam_nbr,totparam_nbr);
|
|
indp2stderrstderr = indp2tottot(1:stderrparam_nbr , 1:stderrparam_nbr);
|
|
indp2stderrcorr = indp2tottot(1:stderrparam_nbr , stderrparam_nbr+1:stderrparam_nbr+corrparam_nbr);
|
|
indp2modmod = indp2tottot(stderrparam_nbr+corrparam_nbr+1:stderrparam_nbr+corrparam_nbr+modparam_nbr , stderrparam_nbr+corrparam_nbr+1:stderrparam_nbr+corrparam_nbr+modparam_nbr);
|
|
if totparam_nbr ~=1
|
|
indp2tottot2 = dyn_vech(indp2tottot); %index of unique second-order derivatives
|
|
else
|
|
indp2tottot2 = indp2tottot;
|
|
end
|
|
if modparam_nbr ~= 1
|
|
indp2modmod2 = dyn_vech(indp2modmod); %get rid of cross derivatives
|
|
else
|
|
indp2modmod2 = indp2modmod;
|
|
end
|
|
end
|
|
endo_nbr = size(A,1); %number of all declared endogenous variables
|
|
var_nbr = length(indvar); %number of considered variables
|
|
exo_nbr = size(B,2); %number of exogenous shocks in model
|
|
|
|
if kronflag == -1
|
|
% numerical two-sided finite difference method using function identification_numerical_objective.m (previously thet2tau.m) for Jacobian (wrt parameters) of A, B, Sig, Om, Yss, and g1
|
|
para0 = get_all_parameters(estim_params, M); %get all selected parameters in estimated_params block, stderr and corr come first, then model parameters
|
|
if isempty(para0)
|
|
%if there is no estimated_params block, consider all stderr and all model parameters, but no corr parameters
|
|
para0 = [stderr_e0', params0'];
|
|
end
|
|
%Jacobians (wrt paramters) of steady state, solution matrices A and B, as well as Sigma_e for ALL variables [outputflag=0]
|
|
dYssABSige = fjaco('identification_numerical_objective', para0, 0, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvar);
|
|
M.params = params0; %make sure values are set back
|
|
M.Sigma_e = Sigma_e0; %make sure values are set back
|
|
M.Correlation_matrix = Corr_e0 ; %make sure values are set back
|
|
% get Jacobians for Yss, A, and B from dYssABSige
|
|
indYss = 1:var_nbr;
|
|
indA = (var_nbr+1):(var_nbr+var_nbr^2);
|
|
indB = (var_nbr+var_nbr^2+1):(var_nbr+var_nbr^2+var_nbr*exo_nbr);
|
|
indSigma_e = (var_nbr+var_nbr^2+var_nbr*exo_nbr+1):(var_nbr+var_nbr^2+var_nbr*exo_nbr+exo_nbr*(exo_nbr+1)/2);
|
|
dYss = dYssABSige(indYss , stderrparam_nbr+corrparam_nbr+1:end); %in tensor notation only wrt model parameters
|
|
dA = reshape(dYssABSige(indA , :) , [var_nbr var_nbr totparam_nbr]); %in tensor notation
|
|
dB = reshape(dYssABSige(indB , :) , [var_nbr exo_nbr totparam_nbr]); %in tensor notation
|
|
|
|
dOm = zeros(var_nbr,var_nbr,totparam_nbr); %initialize in tensor notation
|
|
dSigma_e = zeros(exo_nbr,exo_nbr,totparam_nbr); %initialize in tensor notation
|
|
% get Jacobians of Sigma_e and Om wrt stderr parameters
|
|
if ~isempty(indpstderr)
|
|
for jp=1:stderrparam_nbr
|
|
dSigma_e(:,:,jp) = dyn_unvech(dYssABSige(indSigma_e , jp));
|
|
dOm(:,:,jp) = B*dSigma_e(:,:,jp)*B'; %note that derivatives of B wrt stderr parameters are zero by construction
|
|
end
|
|
end
|
|
% get Jacobians of Sigma_e and Om wrt corr parameters
|
|
if ~isempty(indpcorr)
|
|
for jp=1:corrparam_nbr
|
|
dSigma_e(:,:,stderrparam_nbr+jp) = dyn_unvech(dYssABSige(indSigma_e , stderrparam_nbr+jp));
|
|
dOm(:,:,stderrparam_nbr+jp) = B*dSigma_e(:,:,stderrparam_nbr+jp)*B'; %note that derivatives of B wrt corr parameters are zero by construction
|
|
end
|
|
end
|
|
% get Jacobian of Om wrt model parameters
|
|
if ~isempty(indpmodel)
|
|
for jp=1:modparam_nbr
|
|
dOm(:,:,stderrparam_nbr+corrparam_nbr+jp) = dB(:,:,stderrparam_nbr+corrparam_nbr+jp)*Sigma_e0*B' + B*Sigma_e0*dB(:,:,stderrparam_nbr+corrparam_nbr+jp)'; %note that derivatives of Sigma_e wrt model parameters are zero by construction
|
|
end
|
|
end
|
|
|
|
%Jacobian (wrt model parameters ONLY) of steady state and of Jacobian of all dynamic model equations [outputflag=-1]
|
|
dYssg1 = fjaco('identification_numerical_objective', params0(indpmodel), -1, estim_params, M, oo, options, indpmodel, [], [], (1:endo_nbr)');
|
|
M.params = params0; %make sure values are set back
|
|
M.Sigma_e = Sigma_e0; %make sure values are set back
|
|
M.Correlation_matrix = Corr_e0 ; %make sure values are set back
|
|
dg1 = reshape(dYssg1(endo_nbr+1:end,:),[endo_nbr, length(yy0)+length(ex0), modparam_nbr]); %get rid of steady state and in tensor notation
|
|
|
|
if nargout > 6
|
|
%Hessian (wrt paramters) of steady state, solution matrices A and Om [outputflag=-2]
|
|
% 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
|
|
d2YssAOm = hessian_sparse('identification_numerical_objective', para0', options.gstep, -2, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvar);
|
|
M.params = params0; %make sure values are set back
|
|
M.Sigma_e = Sigma_e0; %make sure values are set back
|
|
M.Correlation_matrix = Corr_e0 ; %make sure values are set back
|
|
|
|
d2A = d2YssAOm(indA , indp2tottot2); %only unique elements
|
|
d2Om = d2YssAOm(indA(end)+1:end , indp2tottot2); %only unique elements
|
|
d2Yss = zeros(var_nbr,modparam_nbr,modparam_nbr); %initialize
|
|
for j = 1:var_nbr
|
|
d2Yss(j,:,:) = dyn_unvech(full(d2YssAOm(j,indp2modmod2))); %full Hessian for d2Yss, note that here we duplicate unique values for model parameters
|
|
end
|
|
clear d2YssAOm
|
|
end
|
|
|
|
return %[END OF MAIN FUNCTION]!!!!!
|
|
end
|
|
|
|
if kronflag == -2
|
|
% numerical two-sided finite difference method to compute numerically
|
|
% dYss, dg1, d2Yss and d2g1, the rest is computed analytically (kronflag=0) below
|
|
modpara0 = params0(indpmodel); %focus only on model parameters for dYss, d2Yss and dg1
|
|
[~, g1] = feval([fname,'.dynamic'], yy0, ex0, params0, ys, 1);
|
|
%g1 is [endo_nbr by (dynamicvar_nbr+exo_nbr)] first derivative (wrt all endogenous, exogenous and auxiliary variables) of dynamic model equations, i.e. df/d[yy0;ex0], in DR order
|
|
if nargout > 6
|
|
% computation of d2Yss and d2g1, i.e. second derivative (wrt. parameters) of Jacobian (wrt endogenous and auxilary variables) of dynamic model [outputflag = -1]
|
|
% 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
|
|
d2Yssg1 = hessian_sparse('identification_numerical_objective', modpara0, options.gstep, -1, estim_params, M, oo, options, indpmodel, [], [], (1:endo_nbr)');
|
|
M.params = params0; %make sure values are set back
|
|
M.Sigma_e = Sigma_e0; %make sure values are set back
|
|
M.Correlation_matrix = Corr_e0 ; %make sure values are set back
|
|
|
|
d2Yss = reshape(full(d2Yssg1(1:endo_nbr,:)), [endo_nbr modparam_nbr modparam_nbr]); %put into tensor notation
|
|
for j=1:endo_nbr
|
|
d2Yss(j,:,:) = dyn_unvech(dyn_vech(d2Yss(j,:,:))); %add duplicate values to full hessian
|
|
end
|
|
d2g1_full = d2Yssg1(endo_nbr+1:end,:);
|
|
%store only nonzero unique entries and the corresponding indices of d2g1:
|
|
% rows: respective derivative term
|
|
% 1st column: equation number of the term appearing
|
|
% 2nd column: column number of variable in Jacobian of the dynamic model
|
|
% 3rd column: number of the first parameter in derivative
|
|
% 4th column: number of the second parameter in derivative
|
|
% 5th column: value of the Hessian term
|
|
ind_d2g1 = find(d2g1_full);
|
|
d2g1 = zeros(length(ind_d2g1),5);
|
|
for j=1:length(ind_d2g1)
|
|
[i1, i2] = ind2sub(size(d2g1_full),ind_d2g1(j));
|
|
[ig1, ig2] = ind2sub(size(g1),i1);
|
|
[ip1, ip2] = ind2sub([modparam_nbr modparam_nbr],i2);
|
|
d2g1(j,:) = [ig1 ig2 ip1 ip2 d2g1_full(ind_d2g1(j))];
|
|
end
|
|
clear d2g1_full;
|
|
end
|
|
%Jacobian (wrt parameters) of steady state and Jacobian of dynamic model equations [outputflag=-1]
|
|
dg1 = fjaco('identification_numerical_objective', modpara0, -1, estim_params, M, oo, options, indpmodel, [], [], (1:endo_nbr)');
|
|
M.params = params0; %make sure values are set back
|
|
M.Sigma_e = Sigma_e0; %make sure values are set back
|
|
M.Correlation_matrix = Corr_e0 ; %make sure values are set back
|
|
dYss = dg1(1:endo_nbr , :);
|
|
dg1 = reshape(dg1(endo_nbr+1 : end , :),[endo_nbr, length(yy0)+length(ex0), modparam_nbr]); %get rid of steady state
|
|
elseif (kronflag == 0 || kronflag == 1)
|
|
% Analytical method to compute dYss, dg1, d2Yss and d2g1
|
|
[~, g1_static] = feval([fname,'.static'], ys, ex0, params0);
|
|
%g1_static is [endo_nbr by endo_nbr] first-derivative (wrt variables) of static model equations f, i.e. df/dys, in declaration order
|
|
rp_static = feval([fname,'.static_params_derivs'], ys, repmat(ex0, maximum_exo_lag+maximum_exo_lead+1), params0);
|
|
%rp_static is [endo_nbr by param_nbr] first-derivative (wrt parameters) of static model equations f, i.e. df/dparams, in declaration order
|
|
dys = -g1_static\rp_static;
|
|
%use implicit function theorem (equation 5 of Ratto and Iskrev (2011) to compute [endo_nbr by param_nbr] first-derivative (wrt parameters) of steady state analytically, note that dys is in declaration order
|
|
d2ys = zeros(length(ys), param_nbr, param_nbr); %initialize in tensor notation
|
|
if nargout > 6
|
|
[~, ~, g2_static] = feval([fname,'.static'], ys, ex0, params0);
|
|
%g2_static is [endo_nbr by endo_nbr^2] second derivative (wrt variables) of static model equations f, i.e. d(df/dys)/dys, in declaration order
|
|
[~, g1, g2, g3] = feval([fname,'.dynamic'], yy0, ex0, params0, ys, 1);
|
|
%g1 is [endo_nbr by (dynamicvar_nbr+exo_nbr)] first derivative (wrt all endogenous, exogenous and auxiliary variables) of dynamic model equations, i.e. df/d[yy0;ex0], in DR order
|
|
%g2 is [endo_nbr by (dynamicvar_nbr+exo_nbr)^2] second derivative (wrt all endogenous, exogenous and auxiliary variables) of dynamic model equations, i.e. d(df/d[yy0;ex0])/d[yy0;ex0], in DR order
|
|
%g3 is [endo_nbr by (dynamicvar_nbr+exo_nbr)^2] third-derivative (wrt all endogenous, exogenous and auxiliary variables) of dynamic model equations, i.e. d(df/d[yy0;ex0])/d[yy0;ex0], in DR order
|
|
[~, gp_static, rpp_static] = feval([fname,'.static_params_derivs'], ys, ex0, params0);
|
|
%gp_static is [endo_nbr by endo_nbr by param_nbr] first derivative (wrt parameters) of first-derivative (wrt variables) of static model equations f, i.e. (df/dys)/dparams, in declaration order
|
|
%rpp_static are nonzero values and corresponding indices of second derivative (wrt parameters) of static model equations f, i.e. d(df/dparams)/dparams, in declaration order
|
|
rpp_static = get_all_resid_2nd_derivs(rpp_static, length(ys), param_nbr); %make full matrix out of nonzero values and corresponding indices
|
|
%rpp_static is [endo_nbr by param_nbr by param_nbr] second derivative (wrt parameters) of static model equations, i.e. d(df/dparams)/dparams, in declaration order
|
|
if isempty(find(g2_static))
|
|
%auxiliary expression on page 8 of Ratto and Iskrev (2011) is zero, i.e. gam = 0
|
|
for j = 1:param_nbr
|
|
%using the implicit function theorem, equation 15 on page 7 of Ratto and Iskrev (2011)
|
|
d2ys(:,:,j) = -g1_static\rpp_static(:,:,j);
|
|
%d2ys is [endo_nbr by param_nbr by param_nbr] second-derivative (wrt parameters) of steady state, i.e. d(dys/dparams)/dparams, in declaration order
|
|
end
|
|
else
|
|
gam = rpp_static*0; %initialize auxiliary expression on page 8 of Ratto and Iskrev (2011)
|
|
for j = 1:endo_nbr
|
|
tmp_gp_static_dys = (squeeze(gp_static(j,:,:))'*dys);
|
|
gam(j,:,:) = transpose(reshape(g2_static(j,:),[endo_nbr endo_nbr])*dys)*dys + tmp_gp_static_dys + tmp_gp_static_dys';
|
|
end
|
|
for j = 1:param_nbr
|
|
%using the implicit function theorem, equation 15 on page 7 of Ratto and Iskrev (2011)
|
|
d2ys(:,:,j) = -g1_static\(rpp_static(:,:,j)+gam(:,:,j));
|
|
%d2ys is [endo_nbr by param_nbr by param_nbr] second-derivative (wrt parameters) of steady state, i.e. d(dys/dparams)/dparams, in declaration order
|
|
end
|
|
clear gp_static g2_static tmp_gp_static_dys gam
|
|
end
|
|
end
|
|
%handling of steady state for nonstationary variables
|
|
if any(any(isnan(dys)))
|
|
[U,T] = schur(g1_static);
|
|
qz_criterium = options.qz_criterium;
|
|
e1 = abs(ordeig(T)) < qz_criterium-1;
|
|
k = sum(e1); % Number of non stationary variables.
|
|
% Number of stationary variables: n = length(e1)-k
|
|
[U,T] = ordschur(U,T,e1);
|
|
T = T(k+1:end,k+1:end);
|
|
%using implicit function theorem, equation 5 of Ratto and Iskrev (2011), in declaration order
|
|
dys = -U(:,k+1:end)*(T\U(:,k+1:end)')*rp_static;
|
|
if nargout > 6
|
|
disp('Computation of d2ys for nonstationary variables is not yet correctly handled if g2_static is nonempty, but continue anyways...')
|
|
for j = 1:param_nbr
|
|
%using implicit function theorem, equation 15 of Ratto and Iskrev (2011), in declaration order
|
|
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]
|
|
end
|
|
end
|
|
end
|
|
if nargout > 6
|
|
[~, gp, ~, gpp, hp] = feval([fname,'.dynamic_params_derivs'], yy0, ex0, params0, ys, 1, dys, d2ys);
|
|
%gp is [endo_nbr by (dynamicvar_nbr + exo_nbr) by param_nbr] first-derivative (wrt parameters) of first-derivative (wrt all endogenous, auxiliary and exogenous variables) of dynamic model equations, i.e. d(df/dvars)/dparam, in DR order
|
|
%gpp are nonzero values and corresponding indices of second-derivative (wrt parameters) of first-derivative (wrt all endogenous, auxiliary and exogenous variables) of dynamic model equations, i.e. d(d(df/dvars)/dparam)/dparam, in DR order
|
|
%hp are nonzero values and corresponding indices of first-derivative (wrt parameters) of second-derivative (wrt all endogenous, auxiliary and exogenous variables) of dynamic model equations, i.e. d(d(df/dvars)/dvars)/dparam, in DR order
|
|
d2Yss = d2ys(oo.dr.order_var,indpmodel,indpmodel);
|
|
%[endo_nbr by mod_param_nbr by mod_param_nbr], i.e. put into DR order and focus only on model parameters
|
|
else
|
|
[~, gp] = feval([fname,'.dynamic_params_derivs'], yy0, repmat(ex0, [maximum_exo_lag+maximum_exo_lead+1,1]), params0, ys, 1, dys, d2ys);
|
|
%gp is [endo_nbr by (dynamicvar_nbr + exo_nbr) by param_nbr] first-derivative (wrt parameters) of first-derivative (wrt all endogenous, auxiliary and exogenous variables) of dynamic model equations, i.e. d(df/dvars)/dparam, in DR order
|
|
[~, g1, g2 ] = feval([fname,'.dynamic'], yy0, repmat(ex0, [maximum_exo_lag+maximum_exo_lead+1,1]), params0, ys, 1);
|
|
%g1 is [endo_nbr by (dynamicvar_nbr+exo_nbr)] first derivative (wrt all endogenous, exogenous and auxiliary variables) of dynamic model equations, i.e. df/d[yy0;ex0], in DR order
|
|
%g2 is [endo_nbr by (dynamicvar_nbr+exo_nbr)^2] second derivative (wrt all endogenous, exogenous and auxiliary variables) of dynamic model equations, i.e. d(df/d[yy0;ex0])/d[yy0;ex0], in DR order
|
|
end
|
|
yy0ex0_nbr = sqrt(size(g2,2)); % number of dynamic variables + exogenous variables (length(yy0)+length(ex0))
|
|
dYss = dys(oo.dr.order_var, indpmodel); %focus only on model parameters, note dys is in declaration order, dYss is in DR-order
|
|
dyy0 = dys(I,:);
|
|
yy0_nbr = max(max(lead_lag_incidence)); % retrieve the number of states excluding columns for shocks
|
|
% Computation of dg1, i.e. first derivative (wrt. parameters) of Jacobian (wrt endogenous and auxilary variables) of dynamic model using the implicit function theorem
|
|
% Let g1 denote the Jacobian of dynamic model equations, i.e. g1 = df/d[yy0ex0], evaluated at the steady state
|
|
% Let dg1 denote the first-derivative (wrt parameters) of g1 evaluated at the steady state
|
|
% Note that g1 is a function of both the parameters and of the steady state, which also depends on the parameters.
|
|
% Hence, implicitly g1=g1(p,yy0ex0(p)) and dg1 consists of two parts (see Ratto and Iskrev (2011) formula 7):
|
|
% (1) direct derivative wrt to parameters given by the preprocessor, i.e. gp
|
|
% and
|
|
% (2) contribution of derivative of steady state (wrt parameters), i.e. g2*dyy0
|
|
% Note that in a stochastic context ex0 is always zero and hence can be skipped in the computations
|
|
dg1_part2 = gp*0; %initialize part 2, it has dimension [endo_nbr by (dynamicvar_nbr+exo_nbr) by param_nbr]
|
|
for j = 1:endo_nbr
|
|
[II, JJ] = ind2sub([yy0ex0_nbr yy0ex0_nbr], find(g2(j,:)));
|
|
%g2 is [endo_nbr by (dynamicvar_nbr+exo_nbr)^2]
|
|
for i = 1:yy0ex0_nbr
|
|
is = find(II==i);
|
|
is = is(find(JJ(is)<=yy0_nbr)); %focus only on yy0 derivatives as ex0 variables are 0 in a stochastic context
|
|
if ~isempty(is)
|
|
tmp_g2 = full(g2(j,find(g2(j,:))));
|
|
dg1_part2(j,i,:) = tmp_g2(is)*dyy0(JJ(is),:); %put into tensor notation
|
|
end
|
|
end
|
|
end
|
|
dg1 = gp + dg1_part2; %dg is sum of two parts due to implicit function theorem
|
|
dg1 = dg1(:,:,indpmodel); %focus only on model parameters
|
|
|
|
if nargout > 6
|
|
% Computation of d2g1, i.e. second derivative (wrt. parameters) of Jacobian (wrt endogenous and auxilary variables) of dynamic model using the implicit function theorem
|
|
% Let g1 denote the Jacobian of dynamic model equations, i.e. g1 = df/d[yy0ex0], evaluated at the steady state
|
|
% Let d2g1 denote the second-derivative (wrt parameters) of g1
|
|
% Note that g1 is a function of both the parameters and of the steady state, which also depends on the parameters.
|
|
% Hence, implicitly g1=g1(p,yy0ex0(p)) and the first derivative is given by dg1 = gp + g2*dyy0ex0 (see above)
|
|
% Accordingly, d2g1, the second-derivative (wrt parameters), consists of five parts (ignoring transposes, see Ratto and Iskrev (2011) formula 16)
|
|
% (1) d(gp)/dp = gpp
|
|
% (2) d(gp)/dyy0ex0*d(yy0ex0)/dp = hp * dyy0ex0
|
|
% (3) d(g2)/dp * dyy0ex0 = hp * dyy0ex0
|
|
% (4) d(g2)/dyy0ex0*d(dyy0ex0)/dp * dyy0ex0 = g3 * dyy0ex0 * dyy0ex0
|
|
% (5) g2 * d(dyy0ex0)/dp = g2 * d2yy0ex0
|
|
% Note that part 2 and 3 are equivalent besides the use of transpose (see Ratto and Iskrev (2011) formula 16)
|
|
d2g1_full = sparse(endo_nbr*yy0ex0_nbr, param_nbr*param_nbr); %initialize
|
|
dyy0ex0 = sparse([dyy0; zeros(yy0ex0_nbr-yy0_nbr,param_nbr)]); %Jacobian (wrt model parameters) of steady state of dynamic (endogenous and auxiliary) and exogenous variables
|
|
|
|
g3 = unfold_g3(g3, yy0ex0_nbr);
|
|
g3_tmp = reshape(g3,[endo_nbr*yy0ex0_nbr*yy0ex0_nbr yy0ex0_nbr]);
|
|
d2g1_part4_left = sparse(endo_nbr*yy0ex0_nbr*yy0ex0_nbr,param_nbr);
|
|
for j = 1:param_nbr
|
|
%compute first two terms of part 4
|
|
d2g1_part4_left(:,j) = g3_tmp*dyy0ex0(:,j);
|
|
end
|
|
|
|
for j=1:endo_nbr
|
|
%Note that in the following we focus only on dynamic variables as exogenous variables are 0 by construction in a stochastic setting
|
|
d2g1_part5 = reshape(g2(j,:), [yy0ex0_nbr yy0ex0_nbr]);
|
|
d2g1_part5 = d2g1_part5(:,1:yy0_nbr)*reshape(d2ys(I,:,:),[yy0_nbr,param_nbr*param_nbr]);
|
|
for i=1:yy0ex0_nbr
|
|
ind_part4 = sub2ind([endo_nbr yy0ex0_nbr yy0ex0_nbr], ones(yy0ex0_nbr,1)*j ,ones(yy0ex0_nbr,1)*i, (1:yy0ex0_nbr)');
|
|
d2g1_part4 = (d2g1_part4_left(ind_part4,:))'*dyy0ex0;
|
|
d2g1_part2_and_part3 = (get_hess_deriv(hp,j,i,yy0ex0_nbr,param_nbr))'*dyy0ex0;
|
|
d2g1_part1 = get_2nd_deriv_mat(gpp,j,i,param_nbr);
|
|
d2g1_tmp = d2g1_part1 + d2g1_part2_and_part3 + d2g1_part2_and_part3' + d2g1_part4 + reshape(d2g1_part5(i,:,:),[param_nbr param_nbr]);
|
|
d2g1_tmp = d2g1_tmp(indpmodel,indpmodel); %focus only on model parameters
|
|
if any(any(d2g1_tmp))
|
|
ind_d2g1_tmp = find(triu(d2g1_tmp));
|
|
d2g1_full(sub2ind([endo_nbr yy0ex0_nbr],j,i), ind_d2g1_tmp) = transpose(d2g1_tmp(ind_d2g1_tmp));
|
|
end
|
|
end
|
|
end
|
|
clear d2g1_tmp d2g1_part1 d2g1_part2_and_part3 d2g1_part4 d2g1_part4_left d2g1_part5
|
|
%store only nonzero entries and the corresponding indices of d2g1:
|
|
% rows: respective derivative term
|
|
% 1st column: equation number of the term appearing
|
|
% 2nd column: column number of variable in Jacobian of the dynamic model
|
|
% 3rd column: number of the first parameter in derivative
|
|
% 4th column: number of the second parameter in derivative
|
|
% 5th column: value of the Hessian term
|
|
ind_d2g1 = find(d2g1_full);
|
|
d2g1 = zeros(length(ind_d2g1),5);
|
|
for j=1:length(ind_d2g1)
|
|
[i1, i2] = ind2sub(size(d2g1_full),ind_d2g1(j));
|
|
[ig1, ig2] = ind2sub(size(g1),i1);
|
|
[ip1, ip2] = ind2sub([modparam_nbr modparam_nbr],i2);
|
|
d2g1(j,:) = [ig1 ig2 ip1 ip2 d2g1_full(ind_d2g1(j))];
|
|
end
|
|
clear d2g1_full;
|
|
end
|
|
end
|
|
% clear variables that are not used any more
|
|
clear rp_static g1_static
|
|
clear ys dys dyy0 dyy0ex0
|
|
clear dg1_part2 tmp_g2
|
|
clear g2 gp rpp_static g2_static gp_static d2ys
|
|
clear hp g3 g3_tmp gpp
|
|
clear ind_d2g1 ind_d2g1_tmp ind_part4 i j i1 i2 ig1 ig2 I II JJ ip1 ip2 is
|
|
|
|
% Construct nonzero derivatives wrt to t+1, t, and t-1 variables using kstate
|
|
klen = maximum_endo_lag + maximum_endo_lead + 1; %total length
|
|
k11 = lead_lag_incidence(find([1:klen] ~= maximum_endo_lag+1),:);
|
|
g1nonzero = g1(:,nonzeros(k11'));
|
|
dg1nonzero = dg1(:,nonzeros(k11'),:);
|
|
if nargout > 6
|
|
indind = ismember(d2g1(:,2),nonzeros(k11'));
|
|
tmp = d2g1(indind,:);
|
|
d2g1nonzero = tmp;
|
|
for j = 1:size(tmp,1)
|
|
inxinx = find(nonzeros(k11')==tmp(j,2));
|
|
d2g1nonzero(j,2) = inxinx;
|
|
end
|
|
end
|
|
kstate = oo.dr.kstate;
|
|
|
|
% Construct nonzero derivatives wrt to t+1, i.e. GAM1=-f_{y^+} in Villemot (2011)
|
|
GAM1 = zeros(endo_nbr,endo_nbr);
|
|
dGAM1 = zeros(endo_nbr,endo_nbr,modparam_nbr);
|
|
k1 = find(kstate(:,2) == maximum_endo_lag+2 & kstate(:,3));
|
|
GAM1(:, kstate(k1,1)) = -g1nonzero(:,kstate(k1,3));
|
|
dGAM1(:, kstate(k1,1), :) = -dg1nonzero(:,kstate(k1,3),:);
|
|
if nargout > 6
|
|
indind = ismember(d2g1nonzero(:,2),kstate(k1,3));
|
|
tmp = d2g1nonzero(indind,:);
|
|
tmp(:,end)=-tmp(:,end);
|
|
d2GAM1 = tmp;
|
|
for j = 1:size(tmp,1)
|
|
inxinx = (kstate(k1,3)==tmp(j,2));
|
|
d2GAM1(j,2) = kstate(k1(inxinx),1);
|
|
end
|
|
end
|
|
|
|
% Construct nonzero derivatives wrt to t, i.e. GAM0=f_{y^0} in Villemot (2011)
|
|
[~,cols_b,cols_j] = find(lead_lag_incidence(maximum_endo_lag+1, oo.dr.order_var));
|
|
GAM0 = zeros(endo_nbr,endo_nbr);
|
|
dGAM0 = zeros(endo_nbr,endo_nbr,modparam_nbr);
|
|
GAM0(:,cols_b) = g1(:,cols_j);
|
|
dGAM0(:,cols_b,:) = dg1(:,cols_j,:);
|
|
if nargout > 6
|
|
indind = ismember(d2g1(:,2),cols_j);
|
|
tmp = d2g1(indind,:);
|
|
d2GAM0 = tmp;
|
|
for j = 1:size(tmp,1)
|
|
inxinx = (cols_j==tmp(j,2));
|
|
d2GAM0(j,2) = cols_b(inxinx);
|
|
end
|
|
end
|
|
|
|
% Construct nonzero derivatives wrt to t-1, i.e. GAM2=-f_{y^-} in Villemot (2011)
|
|
k2 = find(kstate(:,2) == maximum_endo_lag+1 & kstate(:,4));
|
|
GAM2 = zeros(endo_nbr,endo_nbr);
|
|
dGAM2 = zeros(endo_nbr,endo_nbr,modparam_nbr);
|
|
GAM2(:, kstate(k2,1)) = -g1nonzero(:,kstate(k2,4));
|
|
dGAM2(:, kstate(k2,1), :) = -dg1nonzero(:,kstate(k2,4),:);
|
|
if nargout > 6
|
|
indind = ismember(d2g1nonzero(:,2),kstate(k2,4));
|
|
tmp = d2g1nonzero(indind,:);
|
|
tmp(:,end) = -tmp(:,end);
|
|
d2GAM2 = tmp;
|
|
for j = 1:size(tmp,1)
|
|
inxinx = (kstate(k2,4)==tmp(j,2));
|
|
d2GAM2(j,2) = kstate(k2(inxinx),1);
|
|
end
|
|
end
|
|
|
|
% Construct nonzero derivatives wrt to u_t, i.e. GAM3=-f_{u} in Villemot (2011)
|
|
GAM3 = -g1(:,length(yy0)+1:end);
|
|
dGAM3 = -dg1(:,length(yy0)+1:end,:);
|
|
if nargout > 6
|
|
cols_ex = [length(yy0)+1:size(g1,2)];
|
|
indind = ismember(d2g1(:,2),cols_ex);
|
|
tmp = d2g1(indind,:);
|
|
tmp(:,end) = -tmp(:,end);
|
|
d2GAM3 = tmp;
|
|
for j = 1:size(tmp,1)
|
|
inxinx = find(cols_ex==tmp(j,2));
|
|
d2GAM3(j,2) = inxinx;
|
|
end
|
|
clear d2g1 d2g1nonzero tmp
|
|
end
|
|
clear cols_b cols_ex cols_j k1 k11 k2 klen kstate
|
|
clear g1nonzero dg1nonzero g1 yy0
|
|
|
|
%% Construct first derivative of Sigma_e
|
|
dSigma_e = zeros(exo_nbr,exo_nbr,totparam_nbr); %initialize
|
|
% note that derivatives wrt model parameters are zero by construction
|
|
% Compute first derivative of Sigma_e wrt stderr parameters (these come first)
|
|
if ~isempty(indpstderr)
|
|
for jp = 1:stderrparam_nbr
|
|
dSigma_e(indpstderr(jp),indpstderr(jp),jp) = 2*stderr_e0(indpstderr(jp));
|
|
if isdiag(Sigma_e0) == 0 % if there are correlated errors add cross derivatives
|
|
indotherex0 = 1:exo_nbr;
|
|
indotherex0(indpstderr(jp)) = [];
|
|
for kk = indotherex0
|
|
dSigma_e(indpstderr(jp), kk, jp) = Corr_e0(indpstderr(jp),kk)*stderr_e0(kk);
|
|
dSigma_e(kk, indpstderr(jp), jp) = dSigma_e(indpstderr(jp), kk, jp); %symmetry
|
|
end
|
|
end
|
|
end
|
|
end
|
|
% Compute first derivative of Sigma_e wrt corr parameters (these come second)
|
|
if ~isempty(indpcorr)
|
|
for jp = 1:corrparam_nbr
|
|
dSigma_e(indpcorr(jp,1),indpcorr(jp,2),stderrparam_nbr+jp) = stderr_e0(indpcorr(jp,1))*stderr_e0(indpcorr(jp,2));
|
|
dSigma_e(indpcorr(jp,2),indpcorr(jp,1),stderrparam_nbr+jp) = dSigma_e(indpcorr(jp,1),indpcorr(jp,2),stderrparam_nbr+jp); %symmetry
|
|
end
|
|
end
|
|
|
|
%% Construct second derivative of Sigma_e
|
|
if nargout > 6
|
|
% note that derivatives wrt (mod x mod) and (corr x corr) parameters
|
|
% are zero by construction; hence we only need to focus on (stderr x stderr), and (stderr x corr)
|
|
d2Sigma_e = zeros(exo_nbr,exo_nbr,totparam_nbr^2); %initialize full matrix, even though we'll reduce it later on to unique upper triangular values
|
|
% Compute upper triangular values of Hessian of Sigma_e wrt (stderr x stderr) parameters
|
|
if ~isempty(indp2stderrstderr)
|
|
for jp = 1:stderrparam_nbr
|
|
for ip = 1:jp
|
|
if jp == ip %same stderr parameters
|
|
d2Sigma_e(indpstderr(jp),indpstderr(jp),indp2stderrstderr(ip,jp)) = 2;
|
|
else %different stderr parameters
|
|
if isdiag(Sigma_e0) == 0 % if there are correlated errors
|
|
d2Sigma_e(indpstderr(jp),indpstderr(ip),indp2stderrstderr(ip,jp)) = Corr_e0(indpstderr(jp),indpstderr(ip));
|
|
d2Sigma_e(indpstderr(ip),indpstderr(jp),indp2stderrstderr(ip,jp)) = Corr_e0(indpstderr(jp),indpstderr(ip)); %symmetry
|
|
end
|
|
end
|
|
end
|
|
end
|
|
end
|
|
% Compute upper triangular values of Hessian of Sigma_e wrt (stderr x corr) parameters
|
|
if ~isempty(indp2stderrcorr)
|
|
for jp = 1:stderrparam_nbr
|
|
for ip = 1:corrparam_nbr
|
|
if indpstderr(jp) == indpcorr(ip,1) %if stderr equal to first index of corr parameter, derivative is equal to stderr corresponding to second index
|
|
d2Sigma_e(indpstderr(jp),indpcorr(ip,2),indp2stderrcorr(jp,ip)) = stderr_e0(indpcorr(ip,2));
|
|
d2Sigma_e(indpcorr(ip,2),indpstderr(jp),indp2stderrcorr(jp,ip)) = stderr_e0(indpcorr(ip,2)); % symmetry
|
|
end
|
|
if indpstderr(jp) == indpcorr(ip,2) %if stderr equal to second index of corr parameter, derivative is equal to stderr corresponding to first index
|
|
d2Sigma_e(indpstderr(jp),indpcorr(ip,1),indp2stderrcorr(jp,ip)) = stderr_e0(indpcorr(ip,1));
|
|
d2Sigma_e(indpcorr(ip,1),indpstderr(jp),indp2stderrcorr(jp,ip)) = stderr_e0(indpcorr(ip,1)); % symmetry
|
|
end
|
|
end
|
|
end
|
|
end
|
|
d2Sigma_e = d2Sigma_e(:,:,indp2tottot2); %focus on upper triangular hessian values
|
|
end
|
|
|
|
|
|
if kronflag == 1
|
|
% The following derivations are based on Iskrev (2010) and its online appendix A.
|
|
% Basic idea is to make use of the implicit function theorem.
|
|
% Let F = GAM0*A - GAM1*A*A - GAM2 = 0
|
|
% Note that F is a function of parameters p and A, which is also a
|
|
% function of p,therefore, F = F(p,A(p)), and hence,
|
|
% dF = Fp + dF_dA*dA or dA = - Fp/dF_dA
|
|
|
|
% Some auxiliary matrices
|
|
I_endo = speye(endo_nbr);
|
|
I_exo = speye(exo_nbr);
|
|
|
|
% Reshape to write derivatives in the Magnus and Neudecker style, i.e. dvec(X)/dp
|
|
dGAM0 = reshape(dGAM0, endo_nbr^2, modparam_nbr);
|
|
dGAM1 = reshape(dGAM1, endo_nbr^2, modparam_nbr);
|
|
dGAM2 = reshape(dGAM2, endo_nbr^2, modparam_nbr);
|
|
dGAM3 = reshape(dGAM3, endo_nbr*exo_nbr, modparam_nbr);
|
|
dSigma_e = reshape(dSigma_e, exo_nbr^2, totparam_nbr);
|
|
|
|
% Compute dA via implicit function
|
|
dF_dA = kron(I_endo,GAM0) - kron(A',GAM1) - kron(I_endo,GAM1*A); %equation 31 in Appendix A of Iskrev (2010)
|
|
Fp = kron(A',I_endo)*dGAM0 - kron( (A')^2,I_endo)*dGAM1 - dGAM2; %equation 32 in Appendix A of Iskrev (2010)
|
|
dA = -dF_dA\Fp;
|
|
|
|
% Compute dB from expressions 33 in Iskrev (2010) Appendix A
|
|
MM = GAM0-GAM1*A; %this corresponds to matrix M in Ratto and Iskrev (2011, page 6) and will be used if nargout > 6 below
|
|
invMM = MM\eye(endo_nbr);
|
|
dB = - kron( (invMM*GAM3)' , invMM ) * ( dGAM0 - kron( A' , I_endo ) * dGAM1 - kron( I_endo , GAM1 ) * dA ) + kron( I_exo, invMM ) * dGAM3 ;
|
|
dBt = commutation(endo_nbr, exo_nbr)*dB; %transose of derivative using the commutation matrix
|
|
|
|
% Add derivatives for stderr and corr parameters, which are zero by construction
|
|
dA = [zeros(endo_nbr^2, stderrparam_nbr+corrparam_nbr) dA];
|
|
dB = [zeros(endo_nbr*exo_nbr, stderrparam_nbr+corrparam_nbr) dB];
|
|
dBt = [zeros(endo_nbr*exo_nbr, stderrparam_nbr+corrparam_nbr) dBt];
|
|
|
|
% Compute dOm = dvec(B*Sig*B') from expressions 34 in Iskrev (2010) Appendix A
|
|
dOm = kron(I_endo,B*Sigma_e0)*dBt + kron(B,B)*dSigma_e + kron(B*Sigma_e0,I_endo)*dB;
|
|
|
|
% Put into tensor notation
|
|
dA = reshape(dA, endo_nbr, endo_nbr, totparam_nbr);
|
|
dB = reshape(dB, endo_nbr, exo_nbr, totparam_nbr);
|
|
dOm = reshape(dOm, endo_nbr, endo_nbr, totparam_nbr);
|
|
dSigma_e = reshape(dSigma_e, exo_nbr, exo_nbr, totparam_nbr);
|
|
if nargout > 6
|
|
% Put back into tensor notation as these will be reused later
|
|
dGAM0 = reshape(dGAM0, endo_nbr, endo_nbr, modparam_nbr);
|
|
dGAM1 = reshape(dGAM1, endo_nbr, endo_nbr, modparam_nbr);
|
|
dGAM2 = reshape(dGAM2, endo_nbr, endo_nbr, modparam_nbr);
|
|
dGAM3 = reshape(dGAM3, endo_nbr, exo_nbr, modparam_nbr);
|
|
dAA = dA(:, :, stderrparam_nbr+corrparam_nbr+1:end); %this corresponds to matrix dA in Ratto and Iskrev (2011, page 6), i.e. derivative of A with respect to model parameters only in tensor notation
|
|
dBB = dB(:, :, stderrparam_nbr+corrparam_nbr+1:end); %dBB is for all endogenous variables, whereas dB is only for selected variables
|
|
N = -GAM1; %this corresponds to matrix N in Ratto and Iskrev (2011, page 6)
|
|
P = A; %this corresponds to matrix P in Ratto and Iskrev (2011, page 6)
|
|
end
|
|
|
|
% Focus only on selected variables
|
|
dYss = dYss(indvar,:);
|
|
dA = dA(indvar,indvar,:);
|
|
dB = dB(indvar,:,:);
|
|
dOm = dOm(indvar,indvar,:);
|
|
|
|
elseif (kronflag == 0 || kronflag == -2)
|
|
% generalized sylvester equation solves MM*dAA+N*dAA*P=Q from Ratto and Iskrev (2011) equation 11 where
|
|
% dAA is derivative of A with respect to model parameters only in tensor notation
|
|
MM = (GAM0-GAM1*A);
|
|
N = -GAM1;
|
|
P = A;
|
|
Q_rightpart = zeros(endo_nbr,endo_nbr,modparam_nbr); %initialize
|
|
Q = Q_rightpart; %initialize and compute matrix Q in Ratto and Iskrev (2011, page 6)
|
|
for j = 1:modparam_nbr
|
|
Q_rightpart(:,:,j) = (dGAM0(:,:,j)-dGAM1(:,:,j)*A);
|
|
Q(:,:,j) = dGAM2(:,:,j)-Q_rightpart(:,:,j)*A;
|
|
end
|
|
%use iterated generalized sylvester equation to compute dAA
|
|
dAA = sylvester3(MM,N,P,Q);
|
|
flag = 1; icount = 0;
|
|
while flag && icount < 4
|
|
[dAA, flag] = sylvester3a(dAA,MM,N,P,Q);
|
|
icount = icount+1;
|
|
end
|
|
|
|
%stderr parameters come first, then corr parameters, model parameters come last
|
|
%note that stderr and corr derivatives are:
|
|
% - zero by construction for A and B
|
|
% - depend only on dSig for Om
|
|
dOm = zeros(var_nbr, var_nbr, totparam_nbr);
|
|
dA = zeros(var_nbr, var_nbr, totparam_nbr);
|
|
dB = zeros(var_nbr, exo_nbr, totparam_nbr);
|
|
if nargout > 6
|
|
dBB = zeros(endo_nbr, exo_nbr, modparam_nbr); %dBB is always for all endogenous variables, whereas dB is only for selected variables
|
|
end
|
|
|
|
%compute derivative of Om=B*Sig*B' that depends on Sig (other part is added later)
|
|
if ~isempty(indpstderr)
|
|
for j = 1:stderrparam_nbr
|
|
BSigjBt = B*dSigma_e(:,:,j)*B';
|
|
dOm(:,:,j) = BSigjBt(indvar,indvar);
|
|
end
|
|
end
|
|
if ~isempty(indpcorr)
|
|
for j = 1:corrparam_nbr
|
|
BSigjBt = B*dSigma_e(:,:,stderrparam_nbr+j)*B';
|
|
dOm(:,:,stderrparam_nbr+j) = BSigjBt(indvar,indvar);
|
|
end
|
|
end
|
|
|
|
%compute derivative of B and the part of Om=B*Sig*B' that depends on B (other part is computed above)
|
|
invMM = inv(MM);
|
|
for j = 1:modparam_nbr
|
|
dAAj = dAA(:,:,j);
|
|
dBj = invMM * ( dGAM3(:,:,j) - (Q_rightpart(:,:,j) -GAM1*dAAj ) * B ); %equation 14 in Ratto and Iskrev (2011), except in the paper there is a typo as the last B is missing
|
|
dOmj = dBj*Sigma_e0*B'+B*Sigma_e0*dBj';
|
|
%store derivatives in tensor notation
|
|
dA(:, :, stderrparam_nbr+corrparam_nbr+j) = dAAj(indvar,indvar);
|
|
dB(:, :, stderrparam_nbr+corrparam_nbr+j) = dBj(indvar,:);
|
|
dOm(:, :, stderrparam_nbr+corrparam_nbr+j) = dOmj(indvar,indvar);
|
|
if nargout > 6
|
|
dBB(:, :, j) = dBj;
|
|
end
|
|
end
|
|
dYss = dYss(indvar,:); % Focus only on relevant variables
|
|
end
|
|
|
|
%% Compute second-order derivatives (wrt params) of solution matrices using generalized sylvester equations, see equations 17 and 18 in Ratto and Iskrev (2011)
|
|
if nargout > 6
|
|
% solves MM*d2AA+N*d2AA*P = QQ where d2AA are second order derivatives (wrt model parameters) of A
|
|
d2Yss = d2Yss(indvar,:,:);
|
|
QQ = zeros(endo_nbr,endo_nbr,floor(sqrt(modparam_nbr2)));
|
|
jcount=0;
|
|
cumjcount=0;
|
|
jinx = [];
|
|
x2x=sparse(endo_nbr*endo_nbr,modparam_nbr2);
|
|
for i=1:modparam_nbr
|
|
for j=1:i
|
|
elem1 = (get_2nd_deriv(d2GAM0,endo_nbr,endo_nbr,j,i)-get_2nd_deriv(d2GAM1,endo_nbr,endo_nbr,j,i)*A);
|
|
elem1 = get_2nd_deriv(d2GAM2,endo_nbr,endo_nbr,j,i)-elem1*A;
|
|
elemj0 = dGAM0(:,:,j)-dGAM1(:,:,j)*A;
|
|
elemi0 = dGAM0(:,:,i)-dGAM1(:,:,i)*A;
|
|
elem2 = -elemj0*dAA(:,:,i)-elemi0*dAA(:,:,j);
|
|
elem2 = elem2 + ( dGAM1(:,:,j)*dAA(:,:,i) + dGAM1(:,:,i)*dAA(:,:,j) )*A;
|
|
elem2 = elem2 + GAM1*( dAA(:,:,i)*dAA(:,:,j) + dAA(:,:,j)*dAA(:,:,i));
|
|
jcount=jcount+1;
|
|
jinx = [jinx; [j i]];
|
|
QQ(:,:,jcount) = elem1+elem2;
|
|
if jcount==floor(sqrt(modparam_nbr2)) || (j*i)==modparam_nbr^2
|
|
if (j*i)==modparam_nbr^2
|
|
QQ = QQ(:,:,1:jcount);
|
|
end
|
|
xx2=sylvester3(MM,N,P,QQ);
|
|
flag=1;
|
|
icount=0;
|
|
while flag && icount<4
|
|
[xx2, flag]=sylvester3a(xx2,MM,N,P,QQ);
|
|
icount = icount + 1;
|
|
end
|
|
x2x(:,cumjcount+1:cumjcount+jcount)=reshape(xx2,[endo_nbr*endo_nbr jcount]);
|
|
cumjcount=cumjcount+jcount;
|
|
jcount = 0;
|
|
jinx = [];
|
|
end
|
|
end
|
|
end
|
|
clear d xx2;
|
|
jcount = 0;
|
|
icount = 0;
|
|
cumjcount = 0;
|
|
MAX_DIM_MAT = 100000000;
|
|
ncol = max(1,floor(MAX_DIM_MAT/(8*var_nbr*(var_nbr+1)/2)));
|
|
ncol = min(ncol, totparam_nbr2);
|
|
d2A = sparse(var_nbr*var_nbr,totparam_nbr2);
|
|
d2Om = sparse(var_nbr*(var_nbr+1)/2,totparam_nbr2);
|
|
d2A_tmp = zeros(var_nbr*var_nbr,ncol);
|
|
d2Om_tmp = zeros(var_nbr*(var_nbr+1)/2,ncol);
|
|
tmpDir = CheckPath('tmp_derivs',dname);
|
|
offset = stderrparam_nbr+corrparam_nbr;
|
|
% d2B = zeros(m,n,tot_param_nbr,tot_param_nbr);
|
|
for j=1:totparam_nbr
|
|
for i=1:j
|
|
jcount=jcount+1;
|
|
if j<=offset %stderr and corr parameters
|
|
y = B*d2Sigma_e(:,:,jcount)*B';
|
|
d2Om_tmp(:,jcount) = dyn_vech(y(indvar,indvar));
|
|
else %model parameters
|
|
jind = j-offset;
|
|
iind = i-offset;
|
|
if i<=offset
|
|
y = dBB(:,:,jind)*dSigma_e(:,:,i)*B'+B*dSigma_e(:,:,i)*dBB(:,:,jind)';
|
|
% y(abs(y)<1.e-8)=0;
|
|
d2Om_tmp(:,jcount) = dyn_vech(y(indvar,indvar));
|
|
else
|
|
icount=icount+1;
|
|
dAAj = reshape(x2x(:,icount),[endo_nbr endo_nbr]);
|
|
% x = get_2nd_deriv(x2x,m,m,iind,jind);%xx2(:,:,jcount);
|
|
elem1 = (get_2nd_deriv(d2GAM0,endo_nbr,endo_nbr,iind,jind)-get_2nd_deriv(d2GAM1,endo_nbr,endo_nbr,iind,jind)*A);
|
|
elem1 = elem1 -( dGAM1(:,:,jind)*dAA(:,:,iind) + dGAM1(:,:,iind)*dAA(:,:,jind) );
|
|
elemj0 = dGAM0(:,:,jind)-dGAM1(:,:,jind)*A-GAM1*dAA(:,:,jind);
|
|
elemi0 = dGAM0(:,:,iind)-dGAM1(:,:,iind)*A-GAM1*dAA(:,:,iind);
|
|
elem0 = elemj0*dBB(:,:,iind)+elemi0*dBB(:,:,jind);
|
|
y = invMM * (get_2nd_deriv(d2GAM3,endo_nbr,exo_nbr,iind,jind)-elem0-(elem1-GAM1*dAAj)*B);
|
|
% d2B(:,:,j+length(indexo),i+length(indexo)) = y;
|
|
% d2B(:,:,i+length(indexo),j+length(indexo)) = y;
|
|
y = y*Sigma_e0*B'+B*Sigma_e0*y'+ ...
|
|
dBB(:,:,jind)*Sigma_e0*dBB(:,:,iind)'+dBB(:,:,iind)*Sigma_e0*dBB(:,:,jind)';
|
|
% x(abs(x)<1.e-8)=0;
|
|
d2A_tmp(:,jcount) = vec(dAAj(indvar,indvar));
|
|
% y(abs(y)<1.e-8)=0;
|
|
d2Om_tmp(:,jcount) = dyn_vech(y(indvar,indvar));
|
|
end
|
|
end
|
|
if jcount==ncol || i*j==totparam_nbr^2
|
|
d2A(:,cumjcount+1:cumjcount+jcount) = d2A_tmp(:,1:jcount);
|
|
% d2A(:,:,j+length(indexo),i+length(indexo)) = x;
|
|
% d2A(:,:,i+length(indexo),j+length(indexo)) = x;
|
|
d2Om(:,cumjcount+1:cumjcount+jcount) = d2Om_tmp(:,1:jcount);
|
|
% d2Om(:,:,j+length(indexo),i+length(indexo)) = y;
|
|
% d2Om(:,:,i+length(indexo),j+length(indexo)) = y;
|
|
save([tmpDir filesep 'd2A_' int2str(cumjcount+1) '_' int2str(cumjcount+jcount) '.mat'],'d2A')
|
|
save([tmpDir filesep 'd2Om_' int2str(cumjcount+1) '_' int2str(cumjcount+jcount) '.mat'],'d2Om')
|
|
cumjcount = cumjcount+jcount;
|
|
jcount=0;
|
|
% d2A = sparse(m1*m1,tot_param_nbr*(tot_param_nbr+1)/2);
|
|
% d2Om = sparse(m1*(m1+1)/2,tot_param_nbr*(tot_param_nbr+1)/2);
|
|
d2A_tmp = zeros(var_nbr*var_nbr,ncol);
|
|
d2Om_tmp = zeros(var_nbr*(var_nbr+1)/2,ncol);
|
|
end
|
|
end
|
|
end
|
|
end
|
|
|
|
return
|
|
|
|
function g22 = get_2nd_deriv(gpp,m,n,i,j)
|
|
% inputs:
|
|
% - gpp: [#second_order_Jacobian_terms by 5] double Hessian matrix (wrt parameters) of a matrix
|
|
% rows: respective derivative term
|
|
% 1st column: equation number of the term appearing
|
|
% 2nd column: column number of variable in Jacobian
|
|
% 3rd column: number of the first parameter in derivative
|
|
% 4th column: number of the second parameter in derivative
|
|
% 5th column: value of the Hessian term
|
|
% - m: scalar number of equations
|
|
% - n: scalar number of variables
|
|
% - i: scalar number for which first parameter
|
|
% - j: scalar number for which second parameter
|
|
|
|
g22=zeros(m,n);
|
|
is=find(gpp(:,3)==i);
|
|
is=is(find(gpp(is,4)==j));
|
|
|
|
if ~isempty(is)
|
|
g22(sub2ind([m,n],gpp(is,1),gpp(is,2)))=gpp(is,5)';
|
|
end
|
|
return
|
|
|
|
function g22 = get_2nd_deriv_mat(gpp,i,j,npar)
|
|
% inputs:
|
|
% - gpp: [#second_order_Jacobian_terms by 5] double Hessian matrix of (wrt parameters) of dynamic Jacobian
|
|
% rows: respective derivative term
|
|
% 1st column: equation number of the term appearing
|
|
% 2nd column: column number of variable in Jacobian of the dynamic model
|
|
% 3rd column: number of the first parameter in derivative
|
|
% 4th column: number of the second parameter in derivative
|
|
% 5th column: value of the Hessian term
|
|
% - i: scalar number for which model equation
|
|
% - j: scalar number for which variable in Jacobian of dynamic model
|
|
% - npar: scalar Number of model parameters, i.e. equals M_.param_nbr
|
|
%
|
|
% output:
|
|
% g22: [npar by npar] Hessian matrix (wrt parameters) of Jacobian of dynamic model for equation i
|
|
% rows: first parameter in Hessian
|
|
% columns: second paramater in Hessian
|
|
|
|
g22=zeros(npar,npar);
|
|
is=find(gpp(:,1)==i);
|
|
is=is(find(gpp(is,2)==j));
|
|
|
|
if ~isempty(is)
|
|
g22(sub2ind([npar,npar],gpp(is,3),gpp(is,4)))=gpp(is,5)';
|
|
end
|
|
return
|
|
|
|
function g22 = get_all_2nd_derivs(gpp,m,n,npar,fsparse)
|
|
|
|
if nargin==4 || isempty(fsparse)
|
|
fsparse=0;
|
|
end
|
|
if fsparse
|
|
g22=sparse(m*n,npar*npar);
|
|
else
|
|
g22=zeros(m,n,npar,npar);
|
|
end
|
|
% c=ones(npar,npar);
|
|
% c=triu(c);
|
|
% ic=find(c);
|
|
|
|
for is=1:length(gpp)
|
|
% d=zeros(npar,npar);
|
|
% d(gpp(is,3),gpp(is,4))=1;
|
|
% indx = find(ic==find(d));
|
|
if fsparse
|
|
g22(sub2ind([m,n],gpp(is,1),gpp(is,2)),sub2ind([npar,npar],gpp(is,3),gpp(is,4)))=gpp(is,5);
|
|
else
|
|
g22(gpp(is,1),gpp(is,2),gpp(is,3),gpp(is,4))=gpp(is,5);
|
|
end
|
|
end
|
|
|
|
return
|
|
|
|
function r22 = get_all_resid_2nd_derivs(rpp,m,npar)
|
|
% inputs:
|
|
% - rpp: [#second_order_residual_terms by 4] double Hessian matrix (wrt paramters) of model equations
|
|
% rows: respective derivative term
|
|
% 1st column: equation number of the term appearing
|
|
% 2nd column: number of the first parameter in derivative
|
|
% 3rd column: number of the second parameter in derivative
|
|
% 4th column: value of the Hessian term
|
|
% - m: scalar Number of residuals (or model equations), i.e. equals endo_nbr
|
|
% - npar: scalar Number of model parameters, i.e. equals param_nbr
|
|
%
|
|
% output:
|
|
% r22: [endo_nbr by param_nbr by param_nbr] Hessian matrix of model equations with respect to parameters
|
|
% rows: equations in order of declaration
|
|
% 1st columns: first parameter number in derivative
|
|
% 2nd columns: second parameter in derivative
|
|
|
|
r22=zeros(m,npar,npar);
|
|
|
|
for is=1:length(rpp)
|
|
% Keep symmetry in hessian, hence 2 and 3 as well as 3 and 2, i.e. d2f/(dp1 dp2) = d2f/(dp2 dp1)
|
|
r22(rpp(is,1),rpp(is,2),rpp(is,3))=rpp(is,4);
|
|
r22(rpp(is,1),rpp(is,3),rpp(is,2))=rpp(is,4);
|
|
end
|
|
|
|
return
|
|
|
|
function h2 = get_all_hess_derivs(hp,r,m,npar)
|
|
|
|
h2=zeros(r,m,m,npar);
|
|
|
|
for is=1:length(hp)
|
|
h2(hp(is,1),hp(is,2),hp(is,3),hp(is,4))=hp(is,5);
|
|
end
|
|
|
|
return
|
|
|
|
function h2 = get_hess_deriv(hp,i,j,m,npar)
|
|
% inputs:
|
|
% - hp: [#first_order_Hessian_terms by 5] double Jacobian matrix (wrt paramters) of dynamic Hessian
|
|
% rows: respective derivative term
|
|
% 1st column: equation number of the term appearing
|
|
% 2nd column: column number of first variable in Hessian of the dynamic model
|
|
% 3rd column: column number of second variable in Hessian of the dynamic model
|
|
% 4th column: number of the parameter in derivative
|
|
% 5th column: value of the Hessian term
|
|
% - i: scalar number for which model equation
|
|
% - j: scalar number for which first variable in Hessian of dynamic model variable
|
|
% - m: scalar Number of dynamic model variables + exogenous vars, i.e. dynamicvar_nbr + exo_nbr
|
|
% - npar: scalar Number of model parameters, i.e. equals M_.param_nbr
|
|
%
|
|
% output:
|
|
% h2: [(dynamicvar_nbr + exo_nbr) by M_.param_nbr] Jacobian matrix (wrt parameters) of dynamic Hessian
|
|
% rows: second dynamic or exogenous variables in Hessian of specific model equation of the dynamic model
|
|
% columns: parameters
|
|
|
|
h2=zeros(m,npar);
|
|
is1=find(hp(:,1)==i);
|
|
is=is1(find(hp(is1,2)==j));
|
|
|
|
if ~isempty(is)
|
|
h2(sub2ind([m,npar],hp(is,3),hp(is,4)))=hp(is,5)';
|
|
end
|
|
|
|
return
|