stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'perturbation_params_perivs' into 'master' 

Added parameter derivatives of perturbation solution up to 3 order

Closes #1595

See merge request Dynare/dynare!1683
time-shift
Sébastien Villemot 2019-12-17 18:17:09 +00:00
parent d5f90c8376 5a8c206760
commit d7f9f83a0f
18 changed files with 2856 additions and 1083 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
matlab/dsge_likelihood.m
View File

@ -111,7 +111,7 @@ function [fval,info,exit_flag,DLIK,Hess,SteadyState,trend_coeff,Model,DynareOpti
%! @sp 2
%! @strong{This function calls:}
%! @sp 1
%! @ref{dynare_resolve}, @ref{lyapunov_symm}, @ref{lyapunov_solver}, @ref{compute_Pinf_Pstar}, @ref{kalman_filter_d}, @ref{missing_observations_kalman_filter_d}, @ref{univariate_kalman_filter_d}, @ref{kalman_steady_state}, @ref{get_first_order_solution_params_deriv}, @ref{kalman_filter}, @ref{score}, @ref{AHessian}, @ref{missing_observations_kalman_filter}, @ref{univariate_kalman_filter}, @ref{priordens}
%! @ref{dynare_resolve}, @ref{lyapunov_symm}, @ref{lyapunov_solver}, @ref{compute_Pinf_Pstar}, @ref{kalman_filter_d}, @ref{missing_observations_kalman_filter_d}, @ref{univariate_kalman_filter_d}, @ref{kalman_steady_state}, @ref{get_perturbation_params_deriv}, @ref{kalman_filter}, @ref{score}, @ref{AHessian}, @ref{missing_observations_kalman_filter}, @ref{univariate_kalman_filter}, @ref{priordens}
%! @end deftypefn
%@eod:
@ -522,11 +522,29 @@ if analytic_derivation
else
indparam=[];
end
if full_Hess
[DT, ~, ~, DOm, DYss, ~, D2T, D2Om, D2Yss] = get_first_order_solution_params_deriv(A, B, EstimatedParameters, Model,DynareResults,DynareOptions,kron_flag,indparam,indexo,[],iv);
else
[DT, ~, ~, DOm, DYss] = get_first_order_solution_params_deriv(A, B, EstimatedParameters, Model,DynareResults,DynareOptions,kron_flag,indparam,indexo,[],iv);
old_order = DynareOptions.order;
if DynareOptions.order > 1%not sure whether this check is necessary
DynareOptions.order = 1; fprintf('Reset order to 1 for analytical parameter derivatives.\n');
end
old_analytic_derivation_mode = DynareOptions.analytic_derivation_mode;
DynareOptions.analytic_derivation_mode = kron_flag;
if full_Hess
DERIVS = get_perturbation_params_derivs(Model, DynareOptions, EstimatedParameters, DynareResults, indparam, indexo, [], true);
indD2T = reshape(1:Model.endo_nbr^2, Model.endo_nbr, Model.endo_nbr);
indD2Om = dyn_unvech(1:Model.endo_nbr*(Model.endo_nbr+1)/2);
D2T = DERIVS.d2KalmanA(indD2T(iv,iv),:);
D2Om = DERIVS.d2Om(dyn_vech(indD2Om(iv,iv)),:);
D2Yss = DERIVS.d2Yss(iv,:,:);
else
DERIVS = get_perturbation_params_derivs(Model, DynareOptions, EstimatedParameters, DynareResults, indparam, indexo, [], false);
end
DT = zeros(Model.endo_nbr, Model.endo_nbr, size(DERIVS.dghx,3));
DT(:,Model.nstatic+(1:Model.nspred),:) = DERIVS.dghx;
DT = DT(iv,iv,:);
DOm = DERIVS.dOm(iv,iv,:);
DYss = DERIVS.dYss(iv,:);
DynareOptions.order = old_order; %make sure order is reset (not sure if necessary)
DynareOptions.analytic_derivation_mode = old_analytic_derivation_mode;%make sure analytic_derivation_mode is reset (not sure if necessary)
else
DT = derivatives_info.DT(iv,iv,:);
DOm = derivatives_info.DOm(iv,iv,:);

7
matlab/dynare_identification.m
View File

@ -284,7 +284,6 @@ options_ = set_default_option(options_,'datafile','');
options_.mode_compute = 0;
options_.plot_priors = 0;
options_.smoother = 1;
options_.options_ident = options_ident; %store identification options into global options
[dataset_, dataset_info, xparam1, hh, M_, options_, oo_, estim_params_, bayestopt_, bounds] = dynare_estimation_init(M_.endo_names, fname, 1, M_, options_, oo_, estim_params_, bayestopt_);
%outputting dataset_ is needed for Octave
@ -294,6 +293,7 @@ options_ident = set_default_option(options_ident,'analytic_derivation_mode', opt
% 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 Jacobians 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 Jacobians are then computed analytically as in options 0
options_.analytic_derivation_mode = options_ident.analytic_derivation_mode; %overwrite setting
% initialize persistent variables in prior_draw
if prior_exist
@ -385,12 +385,13 @@ options_ident = set_default_option(options_ident,'max_dim_subsets_groups',min([4
% In identification_checks_via_subsets.m, when checks_via_subsets=1, this
% option sets the maximum dimension of groups of parameters for which
% the corresponding rank criteria is checked
options_.options_ident = options_ident; %store identification options into global options
MaxNumberOfBytes = options_.MaxNumberOfBytes; %threshold when to save from memory to files
store_options_ident = options_ident;
iload = options_ident.load_ident_files;
SampleSize = options_ident.prior_mc;
if iload <=0
%% Perform new identification analysis, i.e. do not load previous analysis
if prior_exist

33
matlab/fjaco.m
View File

@ -10,7 +10,7 @@ function fjac = fjaco(f,x,varargin)
% OUTPUT
% fjac : finite difference Jacobian
%
% Copyright (C) 2010-2017 Dynare Team
% Copyright (C) 2010-2017,2019 Dynare Team
%
% This file is part of Dynare.
%
@ -29,7 +29,10 @@ function fjac = fjaco(f,x,varargin)
ff=feval(f,x,varargin{:});
tol = eps.^(1/3);
tol = eps.^(1/3); %some default value
if strcmp(func2str(f),'get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification_numerical_objective')
tol= varargin{5}.dynatol.x;
end
h = tol.*max(abs(x),1);
xh1=x+h; xh0=x-h;
h=xh1-xh0;
@ -37,8 +40,34 @@ fjac = NaN(length(ff),length(x));
for j=1:length(x)
xx = x;
xx(j) = xh1(j); f1=feval(f,xx,varargin{:});
if isempty(f1) && (strcmp(func2str(f),'get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification_numerical_objective') )
[~,info]=feval(f,xx,varargin{:});
disp_info_error_identification_perturbation(info,j);
end
xx(j) = xh0(j); f0=feval(f,xx,varargin{:});
if isempty(f0) && (strcmp(func2str(f),'get_perturbation_params_derivs_numerical_objective') || strcmp(func2str(f),'identification_numerical_objective') )
[~,info]=feval(f,xx,varargin{:});
disp_info_error_identification_perturbation(info,j)
end
fjac(:,j) = (f1-f0)/h(j);
end
feval(f,x,varargin{:});
%Auxiliary functions
function disp_info_error_identification_perturbation(info,j)
% there are errors in the solution algorithm
probl_par = get_the_name(j,varargin{5}.TeX,varargin{3},varargin{2},varargin{5});
skipline()
message = get_error_message(info,0,varargin{5});
fprintf('Parameter error in numerical two-sided difference method:\n')
fprintf('Cannot solve the model for %s (info = %d, %s)\n', probl_par, info(1), message);
fprintf('Possible solutions:\n')
fprintf(' -- check your mod file, calibration and steady state computations carefully\n');
fprintf(' -- use analytic derivatives, i.e. set analytic_derivation_mode=0\n');
fprintf(' -- use an estimated_params block without %s or change its value\n', probl_par);
fprintf(' -- change numerical tolerance level in fjaco.m (you can tune ''options_.dynatol.x'' or change fjaco.m function directly)\n');
error('fjaco.m: numerical two-sided difference method yields errors in solution algorithm');
end
end %main function end

959
matlab/get_first_order_solution_params_deriv.m
View File

@ -1,959 +0,0 @@
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

16
matlab/get_identification_jacobians.m
View File

@ -71,7 +71,7 @@ function [J, G, D, dTAU, dLRE, dA, dOm, dYss, MOMENTS] = get_identification_jaco
% * get_minimal_state_representation
% * duplication
% * dyn_vech
% * get_first_order_solution_params_deriv (previously getH)
% * get_perturbation_params_deriv (previously getH)
% * get_all_parameters
% * fjaco
% * lyapunov_symm
@ -124,7 +124,13 @@ exo_nbr = M.exo_nbr;
endo_nbr = M.endo_nbr;
%% Construct dTAU, dLRE, dA, dOm, dYss
[dA, dB, dSig, dOm, dYss, dg1] = get_first_order_solution_params_deriv(A, B, estim_params, M, oo, options, kronflag, indpmodel, indpstderr, indpcorr, (1:endo_nbr)');
DERIVS = get_perturbation_params_derivs(M, options, estim_params, oo, indpmodel, indpstderr, indpcorr, 0);%d2flag=0
dA = zeros(M.endo_nbr, M.endo_nbr, size(DERIVS.dghx,3));
dA(:,M.nstatic+(1:M.nspred),:) = DERIVS.dghx;
dB = DERIVS.dghu;
dSig = DERIVS.dSigma_e;
dOm = DERIVS.dOm;
dYss = DERIVS.dYss;
% Collect terms for derivative of tau=[Yss; vec(A); vech(Om)]
dTAU = zeros(endo_nbr*endo_nbr+endo_nbr*(endo_nbr+1)/2, totparam_nbr);
@ -132,7 +138,7 @@ for j=1:totparam_nbr
dTAU(:,j) = [vec(dA(:,:,j)); dyn_vech(dOm(:,:,j))];
end
dTAU = [ [zeros(endo_nbr, stderrparam_nbr+corrparam_nbr) dYss]; dTAU ]; % add steady state
dLRE = [dYss; reshape(dg1,size(dg1,1)*size(dg1,2),size(dg1,3)) ]; %reshape dg1 and add steady state
dLRE = [dYss; reshape(DERIVS.dg1,size(DERIVS.dg1,1)*size(DERIVS.dg1,2),size(DERIVS.dg1,3)) ]; %reshape dg1 and add steady state
%State Space Matrices for VAROBS variables
C = A(indvobs,:);
@ -144,7 +150,7 @@ dD = dB(indvobs,:,:);
if ~no_identification_moments
if kronflag == -1
%numerical derivative of autocovariogram [outputflag=1]
J = fjaco('identification_numerical_objective', para0, 1, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvobs, useautocorr, nlags, grid_nbr);
J = fjaco(str2func('identification_numerical_objective'), para0, 1, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvobs, useautocorr, nlags, grid_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
@ -254,7 +260,7 @@ if ~no_identification_spectrum
IA = eye(size(A,1));
if kronflag == -1
%numerical derivative of spectral density [outputflag=2]
dOmega_tmp = fjaco('identification_numerical_objective', para0, 2, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvobs, useautocorr, nlags, grid_nbr);
dOmega_tmp = fjaco(str2func('identification_numerical_objective'), para0, 2, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvobs, useautocorr, nlags, grid_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

1533
matlab/get_perturbation_params_derivs.m Normal file
View File

File diff suppressed because it is too large Load Diff

109
matlab/get_perturbation_params_derivs_numerical_objective.m Normal file
View File

@ -0,0 +1,109 @@
function [out,info] = get_perturbation_params_derivs_numerical_objective(params, outputflag, estim_params, M, oo, options)
%function [out,info] = get_perturbation_params_derivs_numerical_objective(params, outputflag, estim_params, M, oo, options)
% -------------------------------------------------------------------------
% Objective function to compute numerically the Jacobians used for get_perturbation_params_derivs
% =========================================================================
% INPUTS
% params: [vector] parameter values at which to evaluate objective function
% stderr parameters come first, corr parameters second, model parameters third
% outputflag: [string] flag which objective to compute (see below)
% estim_params: [structure] storing the estimation information
% M: [structure] storing the model information
% oo: [structure] storing the solution results
% options: [structure] storing the options
% -------------------------------------------------------------------------
%
% OUTPUT (dependent on outputflag and order of approximation):
% - 'perturbation_solution': out = [vec(Sigma_e);vec(ghx);vec(ghu);vec(ghxx);vec(ghxu);vec(ghuu);vec(ghs2);vec(ghxxx);vec(ghxxu);vec(ghxuu);vec(ghuuu);vec(ghxss);vec(ghuss)]
% - 'dynamic_model': out = [Yss; vec(g1); vec(g2); vec(g3)]
% - 'Kalman_Transition': out = [Yss; vec(KalmanA); dyn_vech(KalmanB*Sigma_e*KalmanB')];
% all in DR-order
% -------------------------------------------------------------------------
% This function is called by
% * get_solution_params_deriv.m (previously getH.m)
% * get_identification_jacobians.m (previously getJJ.m)
% -------------------------------------------------------------------------
% This function calls
% * [M.fname,'.dynamic']
% * resol
% * dyn_vech
% =========================================================================
% Copyright (C) 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/>.
% =========================================================================
%% Update stderr, corr and model parameters and compute perturbation approximation and steady state with updated parameters
M = set_all_parameters(params,estim_params,M);
Sigma_e = M.Sigma_e;
[~,info,M,options,oo] = resol(0,M,options,oo);
if info(1) > 0
% there are errors in the solution algorithm
out = [];
return
else
ys = oo.dr.ys; %steady state of model variables in declaration order
ghx = oo.dr.ghx; ghu = oo.dr.ghu;
if options.order > 1
ghxx = oo.dr.ghxx; ghxu = oo.dr.ghxu; ghuu = oo.dr.ghuu; ghs2 = oo.dr.ghs2;
%ghxs = oo.dr.ghxs; ghus = oo.dr.ghus; %these are zero due to certainty equivalence and Gaussian shocks
end
if options.order > 2
ghxxx = oo.dr.ghxxx; ghxxu = oo.dr.ghxxu; ghxuu = oo.dr.ghxuu; ghxss = oo.dr.ghxss; ghuuu = oo.dr.ghuuu; ghuss = oo.dr.ghuss;
%ghxxs = oo.dr.ghxxs; ghxus = oo.dr.ghxus; ghuus = oo.dr.ghuus; ghsss = oo.dr.ghsss; %these are zero due to certainty equivalence and Gaussian shocks
end
end
Yss = ys(oo.dr.order_var); %steady state of model variables in DR order
%% out = [vec(Sigma_e);vec(ghx);vec(ghu);vec(ghxx);vec(ghxu);vec(ghuu);vec(ghs2);vec(ghxxx);vec(ghxxu);vec(ghxuu);vec(ghuuu);vec(ghxss);vec(ghuss)]
if strcmp(outputflag,'perturbation_solution')
out = [Sigma_e(:); ghx(:); ghu(:)];
if options.order > 1
out = [out; ghxx(:); ghxu(:); ghuu(:); ghs2(:);];
end
if options.order > 2
out = [out; ghxxx(:); ghxxu(:); ghxuu(:); ghuuu(:); ghxss(:); ghuss(:)];
end
end
%% out = [Yss; vec(g1); vec(g2); vec(g3)]; of all endogenous variables, in DR order
if strcmp(outputflag,'dynamic_model')
[I,~] = find(M.lead_lag_incidence'); %I is used to evaluate dynamic model files
if options.order == 1
[~, g1] = feval([M.fname,'.dynamic'], ys(I), oo.exo_steady_state', M.params, ys, 1);
out = [Yss; g1(:)];
elseif options.order == 2
[~, g1, g2] = feval([M.fname,'.dynamic'], ys(I), oo.exo_steady_state', M.params, ys, 1);
out = [Yss; g1(:); g2(:)];
elseif options.order == 3
[~, g1, g2, g3] = feval([M.fname,'.dynamic'], ys(I), oo.exo_steady_state', M.params, ys, 1);
g3 = unfold_g3(g3, length(ys(I))+M.exo_nbr);
out = [Yss; g1(:); g2(:); g3(:)];
end
end
%% out = [Yss; vec(KalmanA); dyn_vech(KalmanB*Sigma_e*KalmanB')]; in DR order, where A and B are Kalman transition matrices
if strcmp(outputflag,'Kalman_Transition')
if options.order == 1
KalmanA = zeros(M.endo_nbr,M.endo_nbr);
KalmanA(:,M.nstatic+(1:M.nspred)) = ghx;
Om = ghu*Sigma_e*transpose(ghu);
out = [Yss; KalmanA(:); dyn_vech(Om)];
else
error('''get_perturbation_params_derivs_numerical_objective.m'': Kalman_Transition works only at order=1');
end
end

23
matlab/identification_numerical_objective.m
View File

@ -92,17 +92,6 @@ end
ys = oo.dr.ys; %steady state of model variables in declaration order
y0 = ys(oo.dr.order_var); %steady state of model variables in DR order
%% out = [Yss; vec(A); vec(B); dyn_vech(Sig_e)]; of indvar variables only, in DR order
if outputflag == 0
out = [y0(indvar); vec(A(indvar,indvar)); vec(B(indvar,:)); dyn_vech(M.Sigma_e)];
end
%% out = [Yss; vec(A); dyn_vech(Om)]; of indvar variables only, in DR order
if outputflag == -2
Om = B*M.Sigma_e*transpose(B);
out = [y0(indvar); vec(A(indvar,indvar)); dyn_vech(Om(indvar,indvar))];
end
%% out = [vech(cov(Y_t,Y_t)); vec(cov(Y_t,Y_{t-1}); ...; vec(cov(Y_t,Y_{t-nlags})] of indvar variables, in DR order. This is Iskrev (2010)'s J matrix.
if outputflag == 1
% Denote Ezz0 = E_t(z_t * z_t'), then the following Lyapunov equation defines the autocovariagrom: Ezz0 -A*Ezz*A' = B*Sig_e*B'
@ -145,14 +134,4 @@ if outputflag == 2
kk = kk+1;
out(1 + (kk-1)*var_nbr^2 : kk*var_nbr^2) = g_omega(:);
end
end
%% out = [Yss; vec(g1)]; of all endogenous variables, in DR order
if outputflag == -1
[I,~] = find(M.lead_lag_incidence'); %I is used to evaluate dynamic model files
yy0 = oo.dr.ys(I); %steady state of dynamic model variables in DR order
ex0 = oo.exo_steady_state';
[~, g1] = feval([M.fname,'.dynamic'], yy0, ex0, M.params, ys, 1);
out = [y0; g1(:)];
end
end

46
matlab/unfold_g4.m Normal file
View File

@ -0,0 +1,46 @@
function g4_unfolded = unfold_g4(g4, ny)
% Given the 4th order derivatives stored in a sparse matrix and without
% symmetric elements (as returned by the static/dynamic files) and the number
% of (static or dynamic) variables in the jacobian, returns
% an unfolded version of the same matrix (i.e. with symmetric elements).
% Copyright (C) 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/>.
[i, j, v] = find(g4);
i_unfolded = [];
j_unfolded = [];
v_unfolded = [];
for k = 1:length(v)
l1 = rem(j(k)-1, ny);
j2 = floor((j(k)-1)/ny);
l2 = rem(j2, ny);
j3 = floor((j(k)-1)/ny^2);
l3 = rem(j3,ny);
l4 = floor(j3/ny);
p = unique(perms([l1 l2 l3 l4]), 'rows');
np = rows(p);
i_unfolded = [i_unfolded; repmat(i(k), np, 1)];
j_unfolded = [j_unfolded; 1 + p(:,1) + ny*(p(:,2) + ny*(p(:,3) + ny*p(:,4)))];
v_unfolded = [v_unfolded; repmat(v(k), np, 1)];
end
g4_unfolded = sparse(i_unfolded, j_unfolded, v_unfolded, size(g4, 1), size(g4, 2));

2
tests/.gitignore vendored
View File

@ -33,6 +33,7 @@ wsOct
!/AIM/fs2000_b1L1L_steadystate.m
!/AIM/fsdat.m
!/analytic_derivatives/fsdat_simul.m
!/analytic_derivatives/nBrockMirmanSYM.mat
!/block_bytecode/run_ls2003.m
!/bvar_a_la_sims/bvar_sample.m
!/conditional_forecasts/2/fsdat_simul.m
@ -67,6 +68,7 @@ wsOct
!/gsa/ls2003scr_results.mat
!/gsa/morris/nk_est_data.m
!/identification/as2007/as2007_steadystate.m
!/identification/as2007/G_QT.mat
!/identification/kim/kim2_steadystate.m
!/initval_file/ramst_initval_file_data.m
!/internals/tests.m

9
tests/Makefile.am
View File

@ -19,6 +19,8 @@ MODFILES = \
optimizers/fs2000_9.mod \
optimizers/fs2000_10.mod \
analytic_derivatives/fs2000_analytic_derivation.mod \
analytic_derivatives/BrockMirman_PertParamsDerivs.mod \
analytic_derivatives/burnside_3_order_PertParamsDerivs.mod \
measurement_errors/fs2000_corr_me_ml_mcmc/fs2000_corr_ME.mod \
TeX/fs2000_corr_ME.mod \
estimation/MH_recover/fs2000_recover_tarb.mod \
@ -383,6 +385,7 @@ XFAIL_MODFILES = ramst_xfail.mod \
estimation/fs2000_mixed_ML_xfail.mod \
estimation/fs2000_stochastic_singularity_xfail.mod \
identification/ident_unit_root/ident_unit_root_xfail.mod \
identification/LindeTrabandt/LindeTrabandt2019_xfail.mod \
steady_state/Linear_steady_state_xfail.mod \
optimal_policy/Ramsey/ramsey_histval_xfail.mod \
example1_extra_exo_xfail.mod
@ -696,6 +699,10 @@ identification: m/identification o/identification
m/identification: $(patsubst %.mod, %.m.trs, $(filter identification/%.mod, $(MODFILES)))
o/identification: $(patsubst %.mod, %.o.trs, $(filter identification/%.mod, $(MODFILES)))
analytic_derivatives: m/analytic_derivatives o/analytic_derivatives
m/analytic_derivatives: $(patsubst %.mod, %.m.trs, $(filter analytic_derivatives/%.mod, $(MODFILES)))
o/analytic_derivatives: $(patsubst %.mod, %.o.trs, $(filter analytic_derivatives/%.mod, $(MODFILES)))
fs2000: m/fs2000 o/fs2000
m/fs2000: $(patsubst %.mod, %.m.trs, $(filter fs2000/%.mod, $(MODFILES)))
o/fs2000: $(patsubst %.mod, %.o.trs, $(filter fs2000/%.mod, $(MODFILES)))
@ -813,6 +820,7 @@ EXTRA_DIST = \
data/test.xlsx \
gsa/morris/nk_est_data.m \
analytic_derivatives/fsdat_simul.m \
analytic_derivatives/nBrockMirmanSYM.mat \
fs2000/fsdat_simul.m \
ls2003/data_ca1.m \
measurement_errors/data_ca1.m \
@ -848,6 +856,7 @@ EXTRA_DIST = \
kalman/likelihood_from_dynare/fs2000_estimation_check.inc \
kalman/likelihood_from_dynare/fs2000ns_model.inc \
kalman/likelihood_from_dynare/fs2000ns_estimation_check.inc \
identification/as2007/G_QT.mat \
estimation/fsdat_simul.m \
ep/mean_preserving_spread.m \
decision_rules/example1_results_dyn_432.mat \

474
tests/analytic_derivatives/BrockMirman_PertParamsDerivs.mod Normal file
View File

@ -0,0 +1,474 @@
% "Augmented" Brock Mirman model
% True policy functions and their exact derivatives (ghx,ghu,ghxx,ghxu,ghuu,ghs2,ghxxx,ghxxu,ghxuu,ghuuu,ghxss,ghuss) are computed using Matlab's symbolic toolbox and saved to a mat file
% Created by @wmutschl (Willi Mutschler, willi@mutschler.eu)
% =========================================================================
% Copyright (C) 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/>.
% =========================================================================
@#define CREATE_SYMBOLIC = 0
@#define ORDER = 3
%define parameter values which are used for calibration and estimated_params block
@#define value_SE_A = 0.3
@#define value_SE_X = 0.1
@#define value_SE_Y = 0.9
@#define value_rho_AX = 0.4
@#define value_alph = 0.35
@#define value_betta = 0.99
@#define value_rhoA = 0.9
@#define value_sigA = 0.6
@#define value_sigX = 0.25
@#define value_Xss = 1.2
@#define value_dumA = 2
@#define value_dumK = 1
@#define value_dumepsA = 4
@#define value_dumepsX = 5
%define parameter values which are used for estimated_params block, note that all are different
@#define prior_value_SE_A = 0.8
@#define prior_value_SE_X = 0.6
@#define prior_value_SE_Y = 0.1
@#define prior_value_rho_AX = 0.1
@#define prior_value_alph = 0.15
@#define prior_value_betta = 0.92
@#define prior_value_rhoA = 0.3
@#define prior_value_sigA = 0.3
@#define prior_value_sigX = 0.15
@#define prior_value_Xss = 1.1
@#define prior_value_dumA = 2
@#define prior_value_dumK = 1
@#define prior_value_dumepsA = 5
@#define prior_value_dumepsX = 6
var Y X K C A W Z; %note that declaration order is different from DR order on purpose
varobs C K A;
varexo epsA epsX epsY;
parameters alph betta rhoA sigA sigX Xss dumA dumK dumepsA dumepsX;
%Calibration
alph = @{value_alph};
betta = @{value_betta};
rhoA = @{value_rhoA};
sigA = @{value_sigA};
sigX = @{value_sigX};
Xss = @{value_Xss};
dumA = @{value_dumA};
dumK = @{value_dumK};
dumepsA = @{value_dumepsA};
dumepsX = @{value_dumepsX};
model;
%this is original Brock-Mirman model equations with AR(1) shock
C^(-1) = alph*betta*C(+1)^(-1)*A(+1)*K^(alph-1);
K = A*K(-1)^alph - C;
log(A) = rhoA*log(A(-1)) + sigA*epsA;
Y = A*K(-1)^alph + epsY;
%augmented auxiliary equations to get nonzero ghs2, ghxss and ghuss, they have no economic meaning
log(X) = log(Xss) + sigX*epsX;
W = X(+1)*exp(dumA*A(-1))*exp(dumK*K(-1)); %this will get a nonzero ghs2 and ghxss
Z = X(+1)*exp(dumepsA*epsA)*exp(dumepsX*epsX); %this will get a nonzero ghs2 and ghuss
end;
shocks;
var epsA = (@{value_SE_A})^2;
var epsX = (@{value_SE_X})^2;
var epsA, epsX = @{value_rho_AX}*@{value_SE_A}*@{value_SE_X};
var epsY = (@{value_SE_Y})^2;
end;
steady_state_model;
X = Xss;
A = 1;
K = (alph*betta*A)^(1/(1-alph));
C = A*K^alph-K;
Y = A*K^alph;
W = Xss*exp(dumA*A)*exp(dumK*K);
Z = Xss*exp(dumepsA*0)*exp(dumepsX*0);
end;
estimated_params;%note that the parameter ordering is different than declaration order
rhoA, normal_pdf, @{prior_value_rhoA}, 0.1;
betta, normal_pdf, @{prior_value_betta}, 0.1;
alph, normal_pdf, @{prior_value_alph}, 0.1;
corr epsA, epsX, normal_pdf, @{prior_value_rho_AX}, 0.1;
sigA, normal_pdf, @{prior_value_sigA}, 0.1;
stderr epsA, normal_pdf, @{prior_value_SE_A}, 0.1;
stderr epsX, normal_pdf, @{prior_value_SE_X}, 0.1;
stderr epsY, normal_pdf, @{prior_value_SE_Y}, 0.1;
sigX, normal_pdf, @{prior_value_sigX}, 0.1;
Xss, normal_pdf, @{prior_value_Xss}, 0.1;
dumepsX, normal_pdf, @{prior_value_dumepsX}, 0.1;
dumK, normal_pdf, @{prior_value_dumK}, 0.1;
dumepsA, normal_pdf, @{prior_value_dumepsA}, 0.1;
dumA, normal_pdf, @{prior_value_dumA}, 0.1;
end;
%save calibration parameters as these get overwritten through stoch_simul and identification command
calib_params = M_.params;
calib_Sigma_e = M_.Sigma_e;
stoch_simul(order=@{ORDER},nograph,irf=0,periods=0);
identification(order=@{ORDER},nograph,no_identification_strength);
indpmodel = estim_params_.param_vals(:,1);
indpstderr = estim_params_.var_exo(:,1);
indpcorr = estim_params_.corrx(:,1:2);
[I,~] = find(M_.lead_lag_incidence');
%% Parameter derivatives of perturbation
@#if CREATE_SYMBOLIC == 1
syms Y_ Y0 Yp C_ C0 Cp K_ K0 Kp A_ A0 Ap X_ X0 Xp W_ W0 Wp Z_ Z0 Zp;
syms epsA0 epsX0 epsY0;
syms RHO_AX SE_A SE_X SE_Y ALPH BETTA RHOA SIGA SIGX XSS DUMA DUMK DUMEPSA DUMEPSX;
syms sig;
SYM.corr_params = [RHO_AX];
SYM.stderr_params_decl = [SE_A SE_X SE_Y];
SYM.modparams_decl = [ALPH BETTA RHOA SIGA SIGX XSS DUMA DUMK DUMEPSA DUMEPSX];
SYM.endo_decl_ = [Y_ X_ K_ C_ A_ W_ Z_];
SYM.endo_decl0 = [Y0 X0 K0 C0 A0 W0 Z0];
SYM.endo_declp = [Yp Xp Kp Cp Ap Wp Zp];
SYM.ex0 = [epsA0 epsX0 epsY0];
SYM.model_eqs = [C0^(-1)-ALPH*BETTA*Cp^(-1)*Ap*K0^(ALPH-1);
K0-A0*K_^ALPH+C0;
log(A0)-RHOA*log(A_)-SIGA*epsA0;
Y0 - A0*K_^ALPH - epsY0;
log(X0)-log(XSS)-SIGX*epsX0;
W0 - Xp*exp(DUMA*A_)*exp(DUMK*K_);
Z0 - Xp*exp(DUMEPSA*epsA0)*exp(DUMEPSX*epsX0);
];
SYM.Sigma_e = [SE_A^2 RHO_AX*SE_A*SE_X 0; RHO_AX*SE_A*SE_X SE_X^2 0; 0 0 SE_Y^2];
SYM.Correlation_matrix = [1 RHO_AX 0; RHO_AX 1 0; 0 0 1];
%steady states
SYM.epsAbar = sym(0);
SYM.epsXbar = sym(0);
SYM.epsYbar = sym(0);
SYM.Abar = sym(1);
SYM.Xbar = XSS;
SYM.Kbar = (ALPH*BETTA*SYM.Abar)^(1/(1-ALPH));
SYM.Cbar = SYM.Abar*SYM.Kbar^ALPH-SYM.Kbar;
SYM.Ybar = SYM.Abar*SYM.Kbar^ALPH;
SYM.Wbar = XSS*exp(DUMA*SYM.Abar)*exp(DUMK*SYM.Kbar);
SYM.Zbar = XSS*exp(DUMEPSA*0)*exp(DUMEPSX*0);
SYM.endo_bar_decl = [SYM.Ybar SYM.Xbar SYM.Kbar SYM.Cbar SYM.Abar SYM.Wbar SYM.Zbar];
SYM.ex0bar = [SYM.epsAbar SYM.epsXbar SYM.epsYbar];
%True policy functions (in declaration order): gs is used for derivatives wrt to perturbation parameter
SYM.g = [A_^RHOA*K_^ALPH*exp(SIGA*epsA0) + epsY0; %Y
XSS*exp(SIGX*epsX0); %X
ALPH*BETTA*A_^RHOA*K_^ALPH*exp(SIGA*epsA0); %K
(1-ALPH*BETTA)*A_^RHOA*K_^ALPH*exp(SIGA*epsA0); %C
A_^RHOA*exp(SIGA*epsA0); %A
XSS*exp(SIGX*0)*exp(DUMA*A_)*exp(DUMK*K_); %W
XSS*exp(SIGX*0)*exp(DUMEPSA*epsA0)*exp(DUMEPSX*epsX0); %Z
];
SYM.gs = [A_^RHOA*K_^ALPH*exp(SIGA*sig*epsA0) + sig*epsY0; %Y
XSS*exp(SIGX*sig*epsX0); %X
ALPH*BETTA*A_^RHOA*K_^ALPH*exp(SIGA*sig*epsA0); %K
(1-ALPH*BETTA)*A_^RHOA*K_^ALPH*exp(SIGA*sig*epsA0); %C
A_^RHOA*exp(SIGA*sig*epsA0); %A
XSS*(1+1/2*SIGX^2*sig^2*SE_X^2)*exp(DUMA*A_)*exp(DUMK*K_); %W
XSS*(1+1/2*SIGX^2*sig^2*SE_X^2)*exp(DUMEPSA*epsA0)*exp(DUMEPSX*epsX0); %Z
];
%put into in DR-order
SYM.g = SYM.g(oo_.dr.order_var);
SYM.gs = SYM.gs(oo_.dr.order_var);
SYM.endo_DR_ = SYM.endo_decl_(oo_.dr.order_var);
SYM.endo_DR0 = SYM.endo_decl0(oo_.dr.order_var);
SYM.endo_DRp = SYM.endo_declp(oo_.dr.order_var);
SYM.Yss = SYM.endo_bar_decl(oo_.dr.order_var);
SYM.yy0 = [SYM.endo_decl_(M_.lead_lag_incidence(1,:)~=0) SYM.endo_decl0(M_.lead_lag_incidence(2,:)~=0) SYM.endo_declp(M_.lead_lag_incidence(3,:)~=0)];
SYM.yy0ex0 = [SYM.yy0 SYM.ex0];
SYM.yy0ex0bar = [SYM.endo_bar_decl(I) SYM.ex0bar];
SYM.x = SYM.endo_DR_(M_.nstatic + (1:M_.nspred));
SYM.xbar = SYM.Yss(M_.nstatic + (1:M_.nspred));
SYM.stderr_params = SYM.stderr_params_decl(indpstderr);
SYM.modparams = SYM.modparams_decl(indpmodel);
SYM.params = [SYM.stderr_params SYM.corr_params SYM.modparams];
SYM.yy0ex0_nbr = length(SYM.yy0ex0);
%% Parameter derivatives
SYM.dYss = jacobian(SYM.Yss,SYM.modparams);
SYM.d2Yss = jacobian(SYM.dYss(:),SYM.modparams);
SYM.g1 = jacobian(SYM.model_eqs,SYM.yy0ex0);
SYM.g2 = reshape(jacobian(vec(SYM.g1),SYM.yy0ex0),M_.endo_nbr,SYM.yy0ex0_nbr^2);
for j = 1:M_.endo_nbr
SYM.g3(j,:) = vec(transpose(jacobian(vec(SYM.g2(j,:)),SYM.yy0ex0))); %dynare ordering
end
SYM.dg1 = jacobian(vec(subs(SYM.g1,SYM.yy0ex0,SYM.yy0ex0bar)),SYM.modparams);
SYM.dg2 = jacobian(vec(subs(SYM.g2,SYM.yy0ex0,SYM.yy0ex0bar)),SYM.modparams);
SYM.dg3 = jacobian(vec(subs(SYM.g3,SYM.yy0ex0,SYM.yy0ex0bar)),SYM.modparams);
SYM.dSigma_e = jacobian(SYM.Sigma_e(:),SYM.params);
SYM.dCorrelation_matrix = jacobian(SYM.Correlation_matrix(:),SYM.params);
SYM.ghx = jacobian(SYM.g,SYM.x);
SYM.KalmanA = sym(zeros(M_.endo_nbr,M_.endo_nbr));
SYM.KalmanA(:,M_.nstatic + (1:M_.nspred)) = SYM.ghx;
SYM.dghx = jacobian(vec(subs(SYM.ghx, [SYM.x SYM.ex0], [SYM.xbar SYM.ex0bar])),SYM.params);
SYM.dKalmanA = jacobian(vec(subs(SYM.KalmanA, [SYM.x SYM.ex0], [SYM.xbar SYM.ex0bar])),SYM.params);
SYM.d2KalmanA = jacobian(SYM.dKalmanA(:),SYM.params);
SYM.ghu = jacobian(SYM.g,SYM.ex0);
SYM.dghu = jacobian(vec(subs(SYM.ghu, [SYM.x SYM.ex0], [SYM.xbar SYM.ex0bar])),SYM.params);
SYM.Om = SYM.ghu*SYM.Sigma_e*transpose(SYM.ghu);
SYM.dOm = jacobian(vec(subs(SYM.Om,[SYM.x SYM.ex0], [SYM.xbar SYM.ex0bar])),SYM.params);
SYM.d2Om = jacobian(SYM.dOm(:),SYM.params);
SYM.ghxx = jacobian(SYM.ghx(:),SYM.x);
SYM.ghxx = reshape(SYM.ghxx,M_.endo_nbr,M_.nspred^2);
SYM.dghxx = jacobian(vec(subs(SYM.ghxx, [SYM.x SYM.ex0], [SYM.xbar SYM.ex0bar])),SYM.params);
SYM.ghxu = jacobian(SYM.ghu(:),SYM.x);
SYM.ghxu = reshape(SYM.ghxu,M_.endo_nbr,M_.nspred*M_.exo_nbr);
SYM.dghxu = jacobian(vec(subs(SYM.ghxu, [SYM.x SYM.ex0], [SYM.xbar SYM.ex0bar])),SYM.params);
SYM.ghuu = jacobian(SYM.ghu(:),SYM.ex0);
SYM.ghuu = reshape(SYM.ghuu,M_.endo_nbr,M_.exo_nbr*M_.exo_nbr);
SYM.dghuu = jacobian(vec(subs(SYM.ghuu, [SYM.x SYM.ex0], [SYM.xbar SYM.ex0bar])),SYM.params);
SYM.ghs2 = jacobian(jacobian(SYM.gs,sig),sig);
SYM.dghs2 = jacobian(vec(subs(SYM.ghs2, [SYM.x SYM.ex0], [SYM.xbar SYM.ex0bar])),SYM.params);
for j = 1:M_.endo_nbr
SYM.ghxxx(j,:) = vec(transpose(jacobian(vec(SYM.ghxx(j,:)),SYM.x))); %dynare ordering
SYM.ghxxu(j,:) = vec(transpose(jacobian(vec(SYM.ghxx(j,:)),SYM.ex0))); %dynare ordering
SYM.ghxuu(j,:) = vec(transpose(jacobian(vec(SYM.ghxu(j,:)),SYM.ex0))); %dynare ordering
SYM.ghuuu(j,:) = vec(transpose(jacobian(vec(SYM.ghuu(j,:)),SYM.ex0))); %dynare ordering
SYM.ghxss(j,:) = vec(transpose(jacobian(vec(SYM.ghs2(j,:)),SYM.x))); %dynare ordering
SYM.ghuss(j,:) = vec(transpose(jacobian(vec(SYM.ghs2(j,:)),SYM.ex0))); %dynare ordering
end
SYM.dghxxx = jacobian(vec(subs(SYM.ghxxx, [SYM.x SYM.ex0], [SYM.xbar SYM.ex0bar])),SYM.params);
SYM.dghxxu = jacobian(vec(subs(SYM.ghxxu, [SYM.x SYM.ex0], [SYM.xbar SYM.ex0bar])),SYM.params);
SYM.dghxuu = jacobian(vec(subs(SYM.ghxuu, [SYM.x SYM.ex0], [SYM.xbar SYM.ex0bar])),SYM.params);
SYM.dghuuu = jacobian(vec(subs(SYM.ghuuu, [SYM.x SYM.ex0], [SYM.xbar SYM.ex0bar])),SYM.params);
SYM.dghxss = jacobian(vec(subs(SYM.ghxss, [SYM.x SYM.ex0], [SYM.xbar SYM.ex0bar])),SYM.params);
SYM.dghuss = jacobian(vec(subs(SYM.ghuss, [SYM.x SYM.ex0], [SYM.xbar SYM.ex0bar])),SYM.params);
% Evaluate numerically
for jj = 1:2
if jj == 1
RHO_AX = @{prior_value_rho_AX};
SE_A = @{prior_value_SE_A};
SE_X = @{prior_value_SE_X};
SE_Y = @{prior_value_SE_Y};
ALPH = @{prior_value_alph};
BETTA = @{prior_value_betta};
RHOA = @{prior_value_rhoA};
SIGA = @{prior_value_sigA};
SIGX = @{prior_value_sigX};
XSS = @{prior_value_Xss};
DUMA = @{prior_value_dumA};
DUMK = @{prior_value_dumK};
DUMEPSA = @{prior_value_dumepsA};
DUMEPSX = @{prior_value_dumepsX};
elseif jj == 2
RHO_AX = @{value_rho_AX};
SE_A = @{value_SE_A};
SE_X = @{value_SE_X};
SE_Y = @{value_SE_Y};
ALPH = @{value_alph};
BETTA = @{value_betta};
RHOA = @{value_rhoA};
SIGA = @{value_sigA};
SIGX = @{value_sigX};
XSS = @{value_Xss};
DUMA = @{value_dumA};
DUMK = @{value_dumK};
DUMEPSA = @{value_dumepsA};
DUMEPSX = @{value_dumepsX};
end
sig = 1;
nSYM.Yss = transpose(double(subs(SYM.Yss)));
nSYM.dYss = reshape(double(subs(SYM.dYss)), M_.endo_nbr, length(SYM.modparams));
nSYM.d2Yss = double(reshape(subs(SYM.d2Yss), [M_.endo_nbr length(SYM.modparams) length(SYM.modparams)]));
nSYM.g1 = double(subs(subs(SYM.g1, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dg1 = reshape(double(subs(SYM.dg1)),M_.endo_nbr, SYM.yy0ex0_nbr, length(SYM.modparams));
nSYM.g2 = sparse(double(subs(subs(SYM.g2, SYM.yy0ex0, SYM.yy0ex0bar))));
nSYM.dg2 = reshape(sparse(double(subs(SYM.dg2))), M_.endo_nbr, SYM.yy0ex0_nbr^2*length(SYM.modparams));
nSYM.g3 = sparse(double(subs(subs(SYM.g3, SYM.yy0ex0, SYM.yy0ex0bar))));
nSYM.dg3 = reshape(sparse(double(subs(SYM.dg3))), M_.endo_nbr, SYM.yy0ex0_nbr^3*length(SYM.modparams));
nSYM.Sigma_e = double(subs(SYM.Sigma_e));
nSYM.dSigma_e = double(reshape(subs(SYM.dSigma_e), M_.exo_nbr, M_.exo_nbr,length(SYM.params)));
nSYM.Correlation_matrix = double(subs(SYM.Correlation_matrix));
nSYM.dCorrelation_matrix = double(reshape(subs(SYM.dCorrelation_matrix), M_.exo_nbr, M_.exo_nbr,length(SYM.params)));
nSYM.ghx = double(subs(subs(SYM.ghx, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dghx = reshape(double(subs(SYM.dghx)), M_.endo_nbr, M_.nspred, length(SYM.params));
nSYM.ghu = double(subs(subs(SYM.ghu, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dghu = reshape(double(subs(SYM.dghu)), M_.endo_nbr, M_.exo_nbr, length(SYM.params));
nSYM.Om = double(subs(subs(SYM.Om, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dOm = reshape(double(subs(SYM.dOm)), M_.endo_nbr, M_.endo_nbr, length(SYM.params));
nSYM.d2Om = reshape(sparse(double(subs(SYM.d2Om))), M_.endo_nbr^2, length(SYM.params)^2);
SYM.indp2 = reshape(1:length(SYM.params)^2,length(SYM.params),length(SYM.params));
SYM.indx2 = reshape(1:M_.endo_nbr^2, M_.endo_nbr, M_.endo_nbr);
nSYM.d2Om = nSYM.d2Om(dyn_vech(SYM.indx2),dyn_vech(SYM.indp2));
nSYM.KalmanA = double(subs(subs(SYM.KalmanA, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dKalmanA = reshape(double(subs(SYM.dKalmanA)), M_.endo_nbr, M_.endo_nbr, length(SYM.params));
nSYM.d2KalmanA = reshape(double(subs(SYM.d2KalmanA)), M_.endo_nbr^2, length(SYM.params)^2);
nSYM.d2KalmanA = nSYM.d2KalmanA(:,dyn_vech(SYM.indp2));
nSYM.ghxx = double(subs(subs(SYM.ghxx, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dghxx = reshape(double(subs(SYM.dghxx)), M_.endo_nbr, M_.nspred^2, length(SYM.params));
nSYM.ghxu = double(subs(subs(SYM.ghxu, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dghxu = reshape(double(subs(SYM.dghxu)), M_.endo_nbr, M_.nspred*M_.exo_nbr, length(SYM.params));
nSYM.ghuu = double(subs(subs(SYM.ghuu, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dghuu = reshape(double(subs(SYM.dghuu)), M_.endo_nbr, M_.exo_nbr^2, length(SYM.params));
nSYM.ghs2 = double(subs(subs(SYM.ghs2, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dghs2 = reshape(double(subs(SYM.dghs2)), M_.endo_nbr, length(SYM.params));
nSYM.ghxxx = double(subs(subs(SYM.ghxxx, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dghxxx = reshape(double(subs(SYM.dghxxx)), M_.endo_nbr, M_.nspred^3, length(SYM.params));
nSYM.ghxxu = double(subs(subs(SYM.ghxxu, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dghxxu = reshape(double(subs(SYM.dghxxu)), M_.endo_nbr, M_.nspred^2*M_.exo_nbr, length(SYM.params));
nSYM.ghxuu = double(subs(subs(SYM.ghxuu, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dghxuu = reshape(double(subs(SYM.dghxuu)), M_.endo_nbr, M_.nspred*M_.exo_nbr^2, length(SYM.params));
nSYM.ghuuu = double(subs(subs(SYM.ghuuu, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dghuuu = reshape(double(subs(SYM.dghuuu)), M_.endo_nbr, M_.exo_nbr^3, length(SYM.params));
nSYM.ghxss = double(subs(subs(SYM.ghxss, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dghxss = reshape(double(subs(SYM.dghxss)), M_.endo_nbr, M_.nspred, length(SYM.params));
nSYM.ghuss = double(subs(subs(SYM.ghuss, SYM.yy0ex0, SYM.yy0ex0bar)));
nSYM.dghuss = reshape(double(subs(SYM.dghuss)), M_.endo_nbr, M_.exo_nbr, length(SYM.params));
if jj == 1
nSYMprior = nSYM;
save('nBrockMirmanSYM.mat','nSYMprior')
elseif jj==2
nSYMcalib = nSYM;
save('nBrockMirmanSYM.mat','nSYMcalib','-append')
end
end
@#endif
clc;
tol_vars.Yss = 1e-14;
tol_vars.Sigma_e = 1e-15;
tol_vars.Correlation_matrix = 1e-16;
tol_vars.g1 = 1e-13;
tol_vars.ghx = 1e-13;
tol_vars.ghu = 1e-13;
@#if ORDER > 1
tol_vars.g2 = 1e-12;
tol_vars.ghxx = 1e-12;
tol_vars.ghxu = 1e-12;
tol_vars.ghuu = 1e-12;
tol_vars.ghs2 = 1e-12;
@#endif
@#if ORDER > 2
tol_vars.g3 = 1e-11;
tol_vars.ghxxx = 1e-11;
tol_vars.ghxxu = 1e-11;
tol_vars.ghxuu = 1e-11;
tol_vars.ghuuu = 1e-11;
tol_vars.ghxss = 1e-11;
tol_vars.ghuss = 1e-11;
@#endif
tol_dvars.dYss = [1e-9 1e-9 1e-14 1e-14];
tol_dvars.dSigma_e = [1e-9 1e-15 1e-14 1e-14];
tol_dvars.dCorrelation_matrix = [1e-9 1e-15 1e-14 1e-14];
tol_dvars.dg1 = [1e-6 1e-6 1e-13 1e-13];
tol_dvars.dghx = [1e-8 1e-8 1e-13 1e-13];
tol_dvars.dghu = [1e-8 1e-8 1e-13 1e-13];
@#if ORDER > 1
tol_dvars.dg2 = [1e-5 1e-5 1e-11 1e-11];
tol_dvars.dghxx = [1e-6 1e-6 1e-12 1e-12];
tol_dvars.dghxu = [1e-6 1e-6 1e-12 1e-12];
tol_dvars.dghuu = [1e-6 1e-6 1e-12 1e-12];
tol_dvars.dghs2 = [1e-6 1e-6 1e-12 1e-12];
@#endif
@#if ORDER > 2
tol_dvars.dg3 = [1e-3 1e-3 1e-9 1e-9];
tol_dvars.dghxxx = [1e-5 1e-5 1e-11 1e-11];
tol_dvars.dghxxu = [1e-5 1e-5 1e-11 1e-11];
tol_dvars.dghxuu = [1e-5 1e-5 1e-11 1e-11];
tol_dvars.dghuuu = [1e-5 1e-5 1e-11 1e-11];
tol_dvars.dghxss = [1e-5 1e-5 1e-11 1e-11];
tol_dvars.dghuss = [1e-5 1e-5 1e-11 1e-11];
@#endif
tol_dvars.d2KalmanA = [1e-3 1e-3 1e-13 1e-13];
tol_dvars.d2Om = [1e-3 1e-3 1e-13 1e-13];
tol_dvars.d2Yss = [1e-3 1e-3 1e-13 1e-13];
options_.dynatol.x = eps.^(1/3); %set numerical differentiation step in fjaco.m
for jj = 1:2
lst_vars = {'Yss', 'Sigma_e', 'Correlation_matrix','g1','ghx','ghu'};
@#if ORDER > 1
lst_vars = [lst_vars, 'g2','ghxx','ghxu','ghuu','ghs2'];
@#endif
@#if ORDER > 2
lst_vars = [lst_vars, 'g3','ghxxx','ghxxu','ghxuu','ghuuu','ghxss','ghuss'];
@#endif
lst_dvars = {'dYss','dSigma_e','dg1','dghx','dghu'};
@#if ORDER > 1
lst_dvars = [lst_dvars, 'dg2','dghxx','dghxu','dghuu','dghs2'];
@#endif
@#if ORDER > 2
lst_dvars = [lst_dvars, 'dg3','dghxxx','dghxxu','dghxuu','dghuuu','dghxss','dghuss'];
@#endif
load('nBrockMirmanSYM.mat');
if jj==1
strparamset = 'PRIOR';
nSYM = nSYMprior;
xparam_prior = set_prior(estim_params_,M_,options_);
M_ = set_all_parameters(xparam_prior,estim_params_,M_);
elseif jj==2
strparamset = 'CALIBRATION';
nSYM = nSYMcalib;
xparam1_calib = [];
for j = 1:length(indpstderr)
xparam1_calib = [xparam1_calib; sqrt(calib_Sigma_e(j,j))];
end
for j = 1:size(indpcorr,1)
xparam1_calib = [xparam1_calib; calib_Sigma_e(indpcorr(j,1),indpcorr(j,2))/( sqrt(calib_Sigma_e(indpcorr(j,1),indpcorr(j,1))) * sqrt(calib_Sigma_e(indpcorr(j,2),indpcorr(j,2))) )];
end
xparam1_calib = [xparam1_calib; calib_params(indpmodel)];
M_ = set_all_parameters(xparam1_calib,estim_params_,M_);
end
[~,info,M_,options_,oo_] = resol(0,M_, options_, oo_);
%For convenience we save the objects to compare into oo_.dr.
oo_.dr.Yss = oo_.dr.ys(oo_.dr.order_var);
oo_.dr.Sigma_e = M_.Sigma_e;
oo_.dr.Correlation_matrix = M_.Correlation_matrix;
ex0 = oo_.exo_steady_state';
[~, oo_.dr.g1, oo_.dr.g2, oo_.dr.g3] = feval([M_.fname,'.dynamic'], oo_.dr.ys(I), oo_.exo_steady_state', M_.params, oo_.dr.ys, 1);
oo_.dr.g3 = unfold_g3(oo_.dr.g3, length(oo_.dr.ys(I))+length(oo_.exo_steady_state')); %add symmetric elements to g3
fprintf('***** %s: SOME COMMON OBJECTS *****\n', strparamset)
for id_var = 1:size(lst_vars,2)
dx = norm( nSYM.(sprintf('%s',lst_vars{id_var})) - oo_.dr.(sprintf('%s',lst_vars{id_var})), Inf);
fprintf('Max absolute deviation for %s: %e\n', lst_vars{id_var}, dx);
if dx > tol_vars.(sprintf('%s',lst_vars{id_var}))
error('Something wrong in steady state computation, solution algorithm or preprocessor')
end
end
for d2flag = [0 1];
if d2flag
lst_dvars = [lst_dvars {'d2KalmanA', 'd2Om', 'd2Yss'}];
end
KRONFLAG = [-1 -2 0 1];
for id_kronflag = 1:length(KRONFLAG)
fprintf('***** %s: d2flag=%d and kronflag=%d *****\n',strparamset, d2flag,KRONFLAG(id_kronflag))
options_.analytic_derivation_mode = KRONFLAG(id_kronflag);
DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, oo_, indpmodel, indpstderr, indpcorr, d2flag);
for id_var = 1:size(lst_dvars,2)
dx = norm( vec(nSYM.(sprintf('%s',lst_dvars{id_var}))) - vec(DERIVS.(sprintf('%s',lst_dvars{id_var}))), Inf);
fprintf('Max absolute deviation for %s: %e\n', lst_dvars{id_var}, dx);
if dx > tol_dvars.(sprintf('%s',lst_dvars{id_var}))(id_kronflag)
error('Something wrong in get_perturbation_params_derivs.m')
end
end
end
end
end

230
tests/analytic_derivatives/burnside_3_order_PertParamsDerivs.mod Normal file
View File

@ -0,0 +1,230 @@
% =========================================================================
% Copyright (C) 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/>.
% =========================================================================
/*
Check the policy functions obtained by perturbation at a high approximation
order, using the Burnside (1998, JEDC) model (for which the analytical form of
the policy function is known).
As shown by Burnside, the policy function for yₜ is:
yₜ = βⁱ exp[aᵢ+bᵢ(xₜ−xₛₛ)]
where:
θ² ⎛ 2ρ 1−ρ²ⁱ⎞
— aᵢ = iθxₛₛ + σ² ─────── ⎢i − ────(1−ρⁱ) + ρ² ─────⎥
2(1−ρ)² ⎝ 1−ρ 1−ρ² ⎠
θρ
— bᵢ = ───(1−ρⁱ)
1−ρ
— xₛₛ is the steady state of x
— σ is the standard deviation of e.
With some algebra, it can be shown that the derivative of yₜ at the deterministic
steady state is equal to:
∂ᵐ⁺ⁿ⁺²ᵖ yₜ ∞ (2p)!
──────────────── = ∑ βⁱ bᵢᵐ⁺ⁿ ρᵐ ───── cᵢᵖ exp(iθxₛₛ)
∂ᵐxₜ₋₁ ∂ⁿeₜ ∂²ᵖs ⁱ⁼¹ p!
where:
— s is the stochastic scale factor
θ² ⎛ 2ρ 1−ρ²ⁱ⎞
— cᵢ = ─────── ⎢i − ────(1−ρⁱ) + ρ² ─────⎥
2(1−ρ)² ⎝ 1−ρ 1−ρ² ⎠
Note that derivatives with respect to an odd order for s (i.e. ∂²ᵖ⁺¹s) are always
equal to zero.
The policy function as returned in the oo_.dr.g_* matrices has the following properties:
— its elements are pre-multiplied by the Taylor coefficients;
— derivatives w.r.t. the stochastic scale factor have already been summed up;
— symmetric elements are folded (and they are not pre-multiplied by the number of repetitions).
As a consequence, the element gₘₙ corresponding to the m-th derivative w.r.t.
to xₜ₋₁ and the n-th derivative w.r.t. to eₜ is given by:
1 ∞ cᵢᵖ
gₘₙ = ────── ∑ ∑ βⁱ bᵢᵐ⁺ⁿ ρᵐ ──── exp(iθxₛₛ)
(m+n)! 0≤2p≤k-m-n ⁱ⁼¹ p!
where k is the order of approximation.
*/
@#define ORDER = 3
var y x;
varobs y;
varexo e;
parameters beta theta rho xbar;
xbar = 0.0179;
rho = -0.139;
theta = -1.5;
theta = -10;
beta = 0.95;
model;
y = beta*exp(theta*x(+1))*(1+y(+1));
x = (1-rho)*xbar + rho*x(-1)+e;
end;
shocks;
var e; stderr 0.0348;
end;
steady_state_model;
x = xbar;
y = beta*exp(theta*xbar)/(1-beta*exp(theta*xbar));
end;
estimated_params;
stderr e, normal_pdf, 0.0348,0.01;
beta, normal_pdf, 0.95, 0.01;
theta, normal_pdf, -10, 0.01;
rho, normal_pdf, -0.139, 0.01;
xbar, normal_pdf, 0.0179, 0.01;
end;
steady;check;model_diagnostics;
stoch_simul(order=@{ORDER},k_order_solver,irf=0,drop=0,periods=0,nograph);
identification(order=@{ORDER},nograph,no_identification_strength);
%make sure everything is computed at prior mean
xparam_prior = set_prior(estim_params_,M_,options_);
M_ = set_all_parameters(xparam_prior,estim_params_,M_);
[~,info,M_,options_,oo_] = resol(0,M_, options_, oo_);
indpmodel = estim_params_.param_vals(:,1);
indpstderr = estim_params_.var_exo(:,1);
indpcorr = estim_params_.corrx(:,1:2);
totparam_nbr = length(indpmodel) + length(indpstderr) + size(indpcorr,1);
%% Verify that the policy function coefficients are correct
i = 1:800;
SE_e=sqrt(M_.Sigma_e);
aux1 = rho*(1-rho.^i)/(1-rho);
aux2 = (1-rho.^(2*i))/(1-rho^2);
aux3 = 1/((1-rho)^2);
aux4 = (i-2*aux1+rho^2*aux2);
aux5=aux3*aux4;
b = theta*aux1;
c = 1/2*theta^2*SE_e^2*aux5;
%derivatives wrt to rho only
daux1_drho = zeros(1,length(i));
daux2_drho = zeros(1,length(i));
daux3_drho = 2/((1-rho)^3);
daux4_drho = zeros(1,length(i));
daux5_drho = zeros(1,length(i));
for ii = 1:length(i)
if ii == 1
daux1_drho(ii) = 1;
daux2_drho(ii) = 0;
else
daux1_drho(ii) = rho/(rho^2 - 2*rho + 1) - 1/(rho - 1) + ((ii+1)*rho^ii)/(rho - 1) - rho^(ii+1)/(rho^2 - 2*rho + 1);
daux2_drho(ii) = (2*rho)/(rho^4 - 2*rho^2 + 1) + (2*ii*rho^(2*ii-1))/(rho^2 - 1) - (2*rho^(2*ii+1))/(rho^4 - 2*rho^2 + 1);
end
daux4_drho(ii) = -2*daux1_drho(ii) + 2*rho*aux2(ii) + rho^2*daux2_drho(ii);
daux5_drho(ii) = daux3_drho*aux4(ii) + aux3*daux4_drho(ii);
end
%derivatives of b and c wrt to all parameters
db = zeros(size(b,1),size(b,2),M_.exo_nbr+M_.param_nbr);
db(:,:,3) = aux1;%wrt theta
db(:,:,4) = theta*daux1_drho;%wrt rho
dc = zeros(size(c,1),size(c,2),M_.exo_nbr+M_.param_nbr);
dc(:,:,1) = theta^2*SE_e*aux3*aux4;%wrt SE_e
dc(:,:,3) = theta*SE_e^2*aux3*aux4;%wrt theta
dc(:,:,4) = 1/2*theta^2*SE_e^2*daux5_drho; %wrt rho
d2flag=0;
g_0 = 1/2*oo_.dr.ghs2; if ~isequal(g_0,oo_.dr.g_0); error('something wrong'); end
g_1 = [oo_.dr.ghx oo_.dr.ghu] +3/6*[oo_.dr.ghxss oo_.dr.ghuss]; if ~isequal(g_1,oo_.dr.g_1); error('something wrong'); end
g_2 = 1/2*[oo_.dr.ghxx oo_.dr.ghxu oo_.dr.ghuu]; if ~isequal(g_2,oo_.dr.g_2); error('something wrong'); end
g_3 = 1/6*[oo_.dr.ghxxx oo_.dr.ghxxu oo_.dr.ghxuu oo_.dr.ghuuu]; if ~isequal(g_3,oo_.dr.g_3); error('something wrong'); end
tols = [1e-4 1e-4 1e-12 1e-12];
KRONFLAGS = [-1 -2 0 1];
for k = 1:length(KRONFLAGS)
fprintf('KRONFLAG=%d\n',KRONFLAGS(k));
options_.analytic_derivation_mode = KRONFLAGS(k);
DERIVS = get_perturbation_params_derivs(M_, options_, estim_params_, oo_, indpmodel, indpstderr, indpcorr, d2flag);
oo_.dr.dg_0 = permute(1/2*DERIVS.dghs2,[1 3 2]);
oo_.dr.dg_1 = cat(2,DERIVS.dghx,DERIVS.dghu) + 3/6*cat(2,DERIVS.dghxss,DERIVS.dghuss);
oo_.dr.dg_2 = 1/2*cat(2,DERIVS.dghxx,DERIVS.dghxu,DERIVS.dghuu);
oo_.dr.dg_3 = 1/6*[DERIVS.dghxxx DERIVS.dghxxu DERIVS.dghxuu DERIVS.dghuuu];
for ord = 0:@{ORDER}
g = oo_.dr.(['g_' num2str(ord)])(2,:); % Retrieve computed policy function for variable y
dg = oo_.dr.(['dg_' num2str(ord)])(2,:,:);
for m = 0:ord % m is the derivation order with respect to x(-1)
v = 0;
dv = zeros(1,M_.exo_nbr + M_.param_nbr);
for p = 0:floor((@{ORDER}-ord)/2) % 2p is the derivation order with respect to s
if ord+2*p > 0 % Skip the deterministic steady state constant
v = v + sum(beta.^i.*exp(theta*xbar*i).*b.^ord.*rho^m.*c.^p)/factorial(ord)/factorial(p);
%derivatives
dv(:,1) = dv(:,1) + sum( beta.^i.*exp(theta*xbar*i).*ord.*b.^(ord-1).*db(:,:,1).*rho^m.*c.^p...
+beta.^i.*exp(theta*xbar*i).*b.^ord.*rho^m.*p.*c.^(p-1).*dc(:,:,1)...
)/factorial(ord)/factorial(p);%wrt SE_E
dv(:,2) = dv(:,2) + sum( i.*beta.^(i-1).*exp(theta*xbar*i).*b.^ord.*rho^m.*c.^p...
+beta.^i.*exp(theta*xbar*i).*ord.*b.^(ord-1).*db(:,:,2).*rho^m.*c.^p...
+beta.^i.*exp(theta*xbar*i).*b.^ord.*rho^m.*p.*c.^(p-1).*dc(:,:,2)...
)/factorial(ord)/factorial(p);%wrt beta
dv(:,3) = dv(:,3) + sum( beta.^i.*exp(theta*xbar*i).*xbar.*i.*b.^ord.*rho^m.*c.^p...
+beta.^i.*exp(theta*xbar*i).*ord.*b.^(ord-1).*db(:,:,3).*rho^m.*c.^p...
+beta.^i.*exp(theta*xbar*i).*b.^ord.*rho^m.*p.*c.^(p-1).*dc(:,:,3)...
)/factorial(ord)/factorial(p);%wrt theta
dv(:,4) = dv(:,4) + sum( beta.^i.*exp(theta*xbar*i).*b.^ord.*m.*rho^(m-1).*c.^p...
+beta.^i.*exp(theta*xbar*i).*ord.*b.^(ord-1).*db(:,:,4).*rho^m.*c.^p...
+beta.^i.*exp(theta*xbar*i).*b.^ord.*rho^m.*p.*c.^(p-1).*dc(:,:,4)...
)/factorial(ord)/factorial(p);%wrt rho
dv(:,5) = dv(:,5) + sum( beta.^i.*exp(theta*xbar*i).*theta.*i.*b.^ord.*rho^m.*c.^p...
+beta.^i.*exp(theta*xbar*i).*ord.*b.^(ord-1).*db(:,:,5).*rho^m.*c.^p...
+beta.^i.*exp(theta*xbar*i).*b.^ord.*rho^m.*p.*c.^(p-1).*dc(:,:,5)...
)/factorial(ord)/factorial(p);%wrt xbar
end
end
if abs(v-g(ord+1-m)) > 1e-14
error(['Error in matrix oo_.dr.g_' num2str(ord)])
end
chk_dg = squeeze(dg(:,ord+1-m,:))';
if isempty(indpstderr)
chk_dv = dv(:,M_.exo_nbr+indpmodel);
elseif isempty(indpmodel)
chk_dv = dv(:,1:M_.exo_nbr);
else
chk_dv = dv;
end
fprintf('Max absolute deviation for dg_%d(2,%d,:): %e\n',ord,ord+1-m,norm( chk_dv - chk_dg, Inf));
if norm( chk_dv - chk_dg, Inf) > tols(k)
error(['Error in matrix dg_' num2str(ord)])
chk_dv
chk_dg
end
end
end
fprintf('\n');
end

66
tests/analytic_derivatives/ls2003.mod
View File

@ -1,66 +0,0 @@
var y y_s R pie dq pie_s de A y_obs pie_obs R_obs;
varexo e_R e_q e_ys e_pies e_A;
parameters psi1 psi2 psi3 rho_R tau alpha rr k rho_q rho_A rho_ys rho_pies;
psi1 = 1.54;
psi2 = 0.25;
psi3 = 0.25;
rho_R = 0.5;
alpha = 0.3;
rr = 2.51;
k = 0.5;
tau = 0.5;
rho_q = 0.4;
rho_A = 0.2;
rho_ys = 0.9;
rho_pies = 0.7;
model(linear);
y = y(+1) - (tau +alpha*(2-alpha)*(1-tau))*(R-pie(+1))-alpha*(tau +alpha*(2-alpha)*(1-tau))*dq(+1) + alpha*(2-alpha)*((1-tau)/tau)*(y_s-y_s(+1))-A(+1);
pie = exp(-rr/400)*pie(+1)+alpha*exp(-rr/400)*dq(+1)-alpha*dq+(k/(tau+alpha*(2-alpha)*(1-tau)))*y+k*alpha*(2-alpha)*(1-tau)/(tau*(tau+alpha*(2-alpha)*(1-tau)))*y_s;
pie = de+(1-alpha)*dq+pie_s;
R = rho_R*R(-1)+(1-rho_R)*(psi1*pie+psi2*(y+alpha*(2-alpha)*((1-tau)/tau)*y_s)+psi3*de)+e_R;
dq = rho_q*dq(-1)+e_q;
y_s = rho_ys*y_s(-1)+e_ys;
pie_s = rho_pies*pie_s(-1)+e_pies;
A = rho_A*A(-1)+e_A;
y_obs = y-y(-1)+A;
pie_obs = 4*pie;
R_obs = 4*R;
end;
shocks;
var e_R = 1.25^2;
var e_q = 2.5^2;
var e_A = 1.89;
var e_ys = 1.89;
var e_pies = 1.89;
end;
varobs y_obs R_obs pie_obs dq de;
estimated_params;
psi1 , gamma_pdf,1.5,0.5;
psi2 , gamma_pdf,0.25,0.125;
psi3 , gamma_pdf,0.25,0.125;
rho_R ,beta_pdf,0.5,0.2;
alpha ,beta_pdf,0.3,0.1;
rr ,gamma_pdf,2.5,1;
k , gamma_pdf,0.5,0.25;
tau ,gamma_pdf,0.5,0.2;
rho_q ,beta_pdf,0.4,0.2;
rho_A ,beta_pdf,0.5,0.2;
rho_ys ,beta_pdf,0.8,0.1;
rho_pies,beta_pdf,0.7,0.15;
stderr e_R,inv_gamma_pdf,1.2533,0.6551;
stderr e_q,inv_gamma_pdf,2.5066,1.3103;
stderr e_A,inv_gamma_pdf,1.2533,0.6551;
stderr e_ys,inv_gamma_pdf,1.2533,0.6551;
stderr e_pies,inv_gamma_pdf,1.88,0.9827;
end;
estimation(datafile='../ls2003/data_ca1.m',first_obs=8,nobs=79,mh_nblocks=2,
prefilter=1,mh_jscale=0.5,mh_replic=0, mode_compute=1, analytic_derivation);

BIN
tests/analytic_derivatives/nBrockMirmanSYM.mat Normal file
View File

Binary file not shown.

363
tests/identification/LindeTrabandt/LindeTrabandt2019_xfail.mod Normal file
View File

@ -0,0 +1,363 @@
% Original model by: J. Linde, M. Trabandt (2019: Should We Use Linearized Models to Calculate Fiscal Multipliers?
% Journal of Applied Econometrics, 2018, 33: 937-965. http://dx.doi.org/10.1002/jae.2641
% This version has some additional dynamics for capital and investment
% Created by @wmutschl (Willi Mutschler, willi@mutschler.eu)
% =========================================================================
% Declare endogenous variables
% =========================================================================
var
% Staggered-price economy
c ${c}$ (long_name='consumption')
lam ${\lambda}$ (long_name='lagrange multiplier budget')
n ${n}$ (long_name='labor')
w ${w}$ (long_name='real wage')
r ${r}$ (long_name='interest rate')
pie ${\pi}$ (long_name='inflation')
y ${y}$ (long_name='output')
pstar ${p^\ast}$ (long_name='price dispersion')
vartheta ${\vartheta}$ (long_name='aggregate price index')
mc ${mc}$ (long_name='marginal costs')
s ${s}$ (long_name='auxiliary variable for nonlinear pricing 1')
f ${f}$ (long_name='auxiliary variable for nonlinear pricing 2')
a ${a}$ (long_name='auxiliary variable for nonlinear pricing 3')
ptilde ${\tilde{p}}$ (long_name='reoptimized price')
delta1 ${\Delta_1}$ (long_name='price dispersion 1')
delta2 ${\Delta_2}$ (long_name='price dispersion 2')
delta3 ${\Delta_3}$ (long_name='price dispersion 3')
b ${b}$ (long_name='bonds')
tau ${\tau}$ (long_name='lump-sum tax')
g
nu %consumption preference shock
% Flex-price economy
cpot ${c^{pot}}$ (long_name='flex-price consumption')
npot ${n^{pot}}$ (long_name='flex-price labor')
wpot ${w^{pot}}$ (long_name='flex-price real wage')
rrpot ${r^{pot}}$ (long_name='flex-price interest rate')
ypot ${y^{pot}}$ (long_name='flex-price output')
bpot ${b^{pot}}$ (long_name='flex-price bonds')
taupot ${\tau^{pot}}$ (long_name='flex-price lump-sum tax')
% Added variables for capital and investment
k ${k}$ (long_name='capital')
rk ${r^{K}}$ (long_name='rental rate on capital')
iv ${i}$ (long_name='investment')
q ${q}$ (long_name='Tobins Q')
kpot ${k^{pot}}$ (long_name='flex-price capital')
rkpot ${r^{K,pot}}$ (long_name='flex-price rental rate on capital')
ivpot ${i^{pot}}$ (long_name='flex-price investment')
;
% =========================================================================
% Declare observable variables
% =========================================================================
varobs c y;
% =========================================================================
% Declare exogenous variables
% =========================================================================
varexo
epsg
epsnu
;
% =========================================================================
% Declare parameters
% =========================================================================
parameters
ALPHA ${\alpha}$ (long_name='capital share')
BETA ${\beta}$ (long_name='discount rate')
PIEBAR ${\bar{\pi}}$ (long_name='target inflation rate')
IOTA ${\iota}$ (long_name='indexation')
CHI ${\chi}$ (long_name='inverse frisch elasticity')
THETAP ${\theta_p}$ (long_name='net markup')
XIP ${\xi_p}$ (long_name='Calvo probability')
PSI ${\psi}$ (long_name='Kimbal curvature')
GAMMAPIE ${\gamma_\pi}$ (long_name='Taylor rule parameter inflation')
GAMMAX ${\gamma_x}$ (long_name='Taylor rule parameter output')
DELTA ${\delta}$ (long_name='depreciation rate')
QBAR ${\bar{Q}}$ (long_name='steady state Q')
NUBAR ${\bar{\nu}}$ (long_name='steady state consumption preference shock')
BBAR %government debt to annualized output ratio
GBAR_o_YBAR % government consumption ratio
TAUBAR
TAUNBAR
PHI
RHONU
RHOG
SIGNU
SIGG
;
% =========================================================================
% Calibration
% =========================================================================
BETA = 0.995;
PIEBAR = 1 + 0.005;
ALPHA = 0.3;
NUBAR = 0;%0.01;
BBAR = 0;%0.6;
GBAR_o_YBAR = 0;%0.2;
TAUBAR = 0;
PHI = 0.0101;
IOTA = 0.5;
CHI = 1/0.4;
XIP = 2/3;
PSI = -12.2;
THETAP = 0.1;
GAMMAPIE = 1.5;
GAMMAX = 0.125;
QBAR=1;
DELTA = 0;%0.025;
RHONU = 0.95;
RHOG = 0.95;
SIGG = 0.6/100;
SIGNU = 0.6/100;
TAUNBAR = 0;
% ========================================================================
% Model equations
% =========================================================================
model;
% -------------------------------------------------------------------------
% Auxiliary expressions
% -------------------------------------------------------------------------
#OMEGAP = (1+THETAP)*(1+PSI)/(1+PSI+THETAP*PSI); %gross markup
#gammap = pie^IOTA*PIEBAR^(1-IOTA); %indexation rule
%#gammap = PIEBAR;
% -------------------------------------------------------------------------
% Added equations for capital and investment
% -------------------------------------------------------------------------
[name='foc household wrt i']
lam = lam*q;
[name='foc household wrt k']
lam*q = BETA*lam(+1)*(rk(+1)+ (1-DELTA)*q(+1) );
[name='capital accumulation']
k = (1-DELTA)*k(-1) + iv;
[name='capital demand']
rk=ALPHA*mc*y/k(-1);
[name='flex price foc household wrt k']
(cpot-steady_state(cpot)*nu)^(-1) = BETA*(cpot(+1)-steady_state(cpot)*nu(+1))^(-1)*(rkpot(+1) + (1-DELTA) );
[name='flex price capital accumulation']
kpot = (1-DELTA)*kpot(-1) + ivpot;
[name='flex price capital demand']
rkpot=ALPHA/(1+THETAP)*(npot/kpot(-1))^(1-ALPHA);
% -------------------------------------------------------------------------
% Actual model equations
% -------------------------------------------------------------------------
[name='foc household wrt c (marginal utility)']
(c-steady_state(c)*nu)^(-1) = lam;
[name='foc household wrt l (leisure vs. labor)']
n^CHI = (1-TAUNBAR)*lam*w;
[name='foc household wrt b (Euler equation)']
lam = BETA*r/pie(+1)*lam(+1);
[name='aggregate resource constraint']
y = c + g + iv;
[name='production function']
y=pstar^(-1)*k(-1)^ALPHA*n^(1-ALPHA);
[name='Nonlinear pricing 1']
s = OMEGAP*lam*y*vartheta^(OMEGAP/(OMEGAP-1))*mc + BETA*XIP*(gammap/pie(+1))^(OMEGAP/(1-OMEGAP))*s(+1);
[name='Nonlinear pricing 2']
f = lam*y*vartheta^(OMEGAP/(OMEGAP-1)) + BETA*XIP*(gammap/pie(+1))^(1/(1-OMEGAP))*f(+1);
[name='Nonlinear pricing 3']
a = PSI*(OMEGAP-1)*lam*y + BETA*XIP*(gammap/pie(+1))*a(+1);
[name='Nonlinear pricing 4']
s = f*ptilde -a*ptilde^(1+OMEGAP/(OMEGAP-1));
[name='zero profit condition']
vartheta = 1 + PSI - PSI*delta2;
[name='aggregate price index']
vartheta = delta3;
[name='overall price dispersion']
pstar = 1/(1+PSI)*vartheta^(OMEGAP/(OMEGAP-1))*delta1^(OMEGAP/(1-OMEGAP)) + PSI/(1+PSI);
[name='price dispersion 1']
delta1^(OMEGAP/(1-OMEGAP)) = (1-XIP)*ptilde^(OMEGAP/(1-OMEGAP)) + XIP*(gammap/pie*delta1(-1))^(OMEGAP/(1-OMEGAP));
[name='price dispersion 2']
delta2 = (1-XIP)*ptilde+XIP*gammap/pie*delta2(-1);
[name='price dispersion 3']
delta3^(1/(1-OMEGAP)) = (1-XIP)*ptilde^(1/(1-OMEGAP)) + XIP*(gammap/pie*delta3(-1))^(1/(1-OMEGAP));
[name='marginal costs']
(1-ALPHA)*mc = w*k(-1)^(-ALPHA)*n^ALPHA;
[name='Taylor Rule']
r = steady_state(r)*(pie/PIEBAR)^GAMMAPIE*( (y/steady_state(y))/(ypot/steady_state(ypot)) )^GAMMAX;
[name='government budget']
b = r(-1)/pie*b(-1) + g/steady_state(y) -TAUNBAR*w*n/steady_state(y) - tau;
[name='fiscal rule']
tau = TAUBAR + PHI*(b(-1)-BBAR);
[name='flex price Euler equation']
(cpot-steady_state(cpot)*nu)^(-1) = BETA*rrpot*(cpot(+1)-steady_state(cpot)*nu(+1))^(-1);
[name='flex price foc household wrt l (leisure vs. labor)']
npot^CHI = (1-TAUNBAR)*(cpot-steady_state(cpot)*nu)^(-1)*wpot;
[name='flex price wage']
(1-ALPHA)/(1+THETAP)*kpot(-1)^ALPHA = wpot*npot^ALPHA;
[name='flex price aggregate resource constraint']
ypot = cpot + g + ivpot;
[name='flex price production function']
ypot=kpot(-1)^ALPHA*npot^(1-ALPHA);
[name='flex price government budget']
bpot = rrpot(-1)*bpot(-1) + g/steady_state(y) -TAUNBAR*wpot*npot/steady_state(y) - taupot;
[name='fiscal rule ']
taupot = TAUBAR + PHI*(bpot(-1)-BBAR);
g/steady_state(y) - GBAR_o_YBAR = RHOG*(g(-1)/steady_state(y) - GBAR_o_YBAR) + epsg;
nu - NUBAR = RHONU*(nu(-1) - NUBAR) + epsnu;
end;
% =========================================================================
% Steady state using a steady_state_model block
% =========================================================================
steady_state_model;
OMEGP = (1+THETAP)*(1+PSI)/(1+PSI+THETAP*PSI);
q = QBAR;
pie = PIEBAR;
GAMMAP = PIEBAR^IOTA*PIEBAR^(1-IOTA);
b = BBAR;
tau = TAUBAR;
taun = TAUNBAR;
nu = NUBAR;
r = pie/BETA;
RK = (1/BETA+DELTA-1)*q; % foc k
aux1 = ( (1-XIP)/(1-XIP*(GAMMAP/pie)^(1/(1-OMEGP)) ) )^(1-OMEGP);
aux2 = (1-XIP)/(1-XIP*(GAMMAP/pie));
ptilde = (1+PSI)/ ( aux1 + PSI*aux2 );
delta2 = aux2*ptilde;
delta3 = aux1*ptilde;
vartheta = delta3;
delta1 = ( (1-XIP)/(1-XIP*(GAMMAP/pie)^(OMEGP/(1-OMEGP) )) )^((1-OMEGP)/OMEGP)*ptilde;
pstar = 1/(1+PSI)*vartheta^(OMEGP/(OMEGP-1))*delta1^(OMEGP/(1-OMEGP)) + PSI/(1+PSI);
f_o_a = ( vartheta^(OMEGP/(OMEGP-1)) / (1-BETA*XIP*(GAMMAP/pie)^(1/(1-OMEGP))) ) / ( PSI*(OMEGP-1) /(1-BETA*XIP*GAMMAP/pie) );
s_o_a = ( OMEGP*vartheta^(OMEGP/(OMEGP-1)) / (1-BETA*XIP*(GAMMAP/pie)^(OMEGP/(1-OMEGP))) ) / ( PSI*(OMEGP-1) /(1-BETA*XIP*GAMMAP/pie) );
mc = (f_o_a*ptilde-ptilde^(1+OMEGP/(OMEGP-1)))*s_o_a^(-1);
k_o_n = (RK/(mc*ALPHA))^(1/(ALPHA-1));
w = (1-ALPHA)*mc*k_o_n^ALPHA; % labor demand
y_o_n = pstar^(-1)*k_o_n^ALPHA; % production function
iv_o_n = DELTA*k_o_n;
iv_o_y = iv_o_n / y_o_n;
g_o_y = GBAR_o_YBAR;
c_o_n = y_o_n - iv_o_n - g_o_y; % market clearing
% The labor level, closed-form solution for labor
n = (c_o_n^(-1)*w)^(1/(CHI+1)) ;
% Value of remaining variables
c = c_o_n*n;
rk = RK;
iv = iv_o_n *n;
k = k_o_n*n;
y = y_o_n*n;
lam = n^CHI/w;
g = g_o_y*y;
a = PSI*(OMEGP-1)*y*lam/(1-BETA*XIP*GAMMAP/pie);
f = f_o_a*a;
s = f*ptilde - a*ptilde^(1+OMEGP/(OMEGP-1));
%flex price economy
bpot = BBAR;
taupot = TAUBAR;
taun = TAUNBAR;
rrpot = 1/BETA;
RKpot = (1/BETA+DELTA-1)*QBAR; % foc k
rkpot = RKpot;
kpot_o_npot = (RKpot/(ALPHA/(1+THETAP)))^(1/(ALPHA-1));
wpot = (1-ALPHA)/(1+THETAP)*kpot_o_npot^ALPHA; % labor demand
ypot_o_npot = kpot_o_npot^ALPHA; % production function
ivpot_o_npot = DELTA*kpot_o_npot;
ivpot_o_ypot = ivpot_o_npot / ypot_o_npot;
cpot_o_npot = ypot_o_npot - ivpot_o_npot - g_o_y; % market clearing
% The labor level, closed-form solution for labor
npot = (cpot_o_npot^(-1)*wpot)^(1/(CHI+1)) ;
% Value of remaining variables
cpot = cpot_o_npot*npot;
ivpot = ivpot_o_npot *npot;
kpot = kpot_o_npot*npot;
ypot = ypot_o_npot*npot;
end;
% =========================================================================
% Declare settings for shocks
% =========================================================================
shocks;
var epsg = 1;
var epsnu = 1;
end;
estimated_params;
ALPHA, 0.3, normal_pdf, 0.3 , 0.1;
BETA, 0.995, normal_pdf, 0.995, 0.1;
PIEBAR, 1.005, normal_pdf, 1.005, 0.1;
IOTA, 0.5, normal_pdf, 0.5, 0.1;
CHI, 1/0.4, normal_pdf, 1/0.4, 0.1;
THETAP, 0.1, normal_pdf, 0.1, 0.1;
XIP, 2/3, normal_pdf, 2/3, 0.1;
PSI, -12.2, normal_pdf, -12.2, 0.1;
GAMMAPIE, 1.5, normal_pdf, 1.5, 0.1;
GAMMAX, 0.125, normal_pdf, 0.125, 0.1;
DELTA, 0, normal_pdf, 0, 0.1;
QBAR, 1, normal_pdf, 1, 0.1;
NUBAR, 0, normal_pdf, 0, 0.1;
BBAR, 0, normal_pdf, 0, 0.1;
GBAR_o_YBAR, 0, normal_pdf, 0, 0.1; %commenting this solves the analytic_derivation_mode=-2|-1 problem
TAUBAR, 0, normal_pdf, 0, 0.1;
TAUNBAR, 0, normal_pdf, 0, 0.1;
PHI, 0.0101, normal_pdf, 0.0101, 0.1;
RHONU, 0.95, normal_pdf, 0.95, 0.1;
RHOG, 0.95, normal_pdf, 0.95, 0.1;
SIGNU, 0.6/100, normal_pdf, 0.6/100, 0.1;
SIGG, 0.6/100, normal_pdf, 0.6/100, 0.1;
end;
% =========================================================================
% Computations
% =========================================================================
model_diagnostics;
steady;
resid;
check;
identification(order=1,no_identification_strength,analytic_derivation_mode= 0,ar=5); %works
%identification(no_identification_strength,analytic_derivation_mode= 1,ar=5); %works, this takes the longest due to Kronecker products
options_.dynatol.x = 1e-9;
identification(order=1,no_identification_strength,analytic_derivation_mode=-1,ar=5); %works, but tolerance is way too small
identification(order=1,no_identification_strength,analytic_derivation_mode=-2,ar=5); %works, but tolerance is way to small
options_.dynatol.x = 1e-5; %this is the default value
identification(order=1,no_identification_strength,analytic_derivation_mode=-1,ar=5); %works only if GBAR_o_YBAR is uncommented
identification(order=1,no_identification_strength,analytic_derivation_mode=-2,ar=5); %works only if GBAR_o_YBAR is uncommented

BIN
tests/identification/as2007/G_QT.mat Normal file
View File

Binary file not shown.

41
tests/identification/as2007/as2007_QT.mod
View File

@ -74,24 +74,23 @@ identification(parameter_set=calibration,
checks_via_subsets=1,
max_dim_subsets_groups=4);
% load('G_QT'); %note that this is computed using replication files of Qu and Tkachenko (2012)
% temp = load([M_.dname filesep 'identification' filesep M_.fname '_identif']);
% G_dynare = temp.ide_spectrum_point.G;
%
% % Compare signs
% if ~isequal(sign(G_dynare),sign(G_QT))
% error('signs of normalized G are note equal');
% end
%
% % Compare normalized versions
% tilda_G_dynare = temp.ide_spectrum_point.tilda_G;
% ind_G_QT = (find(max(abs(G_QT'),[],1) > temp.store_options_ident.tol_deriv));
% tilda_G_QT = zeros(size(G_QT));
% delta_G_QT = sqrt(diag(G_QT(ind_G_QT,ind_G_QT)));
% tilda_G_QT(ind_G_QT,ind_G_QT) = G_QT(ind_G_QT,ind_G_QT)./((delta_G_QT)*(delta_G_QT'));
% if ~isequal(rank(tilda_G_QT,temp.store_options_ident.tol_rank),rank(tilda_G_dynare,temp.store_options_ident.tol_rank))
% error('ranks are not the same for normalized version')
% end
%
% max(max(abs(abs(tilda_G_dynare)-abs(tilda_G_QT))))
% norm(tilda_G_dynare - tilda_G_QT)
load('G_QT'); %note that this is computed using replication files of Qu and Tkachenko (2012)
temp = load([M_.dname filesep 'identification' filesep M_.fname '_identif']);
G_dynare = temp.ide_spectrum_point.G;
% Compare signs
if ~isequal(sign(G_dynare),sign(G_QT))
error('signs of normalized G are note equal');
end
% Compare normalized versions
tilda_G_dynare = temp.ide_spectrum_point.tilda_G;
ind_G_QT = (find(max(abs(G_QT'),[],1) > temp.store_options_ident.tol_deriv));
tilda_G_QT = zeros(size(G_QT));
delta_G_QT = sqrt(diag(G_QT(ind_G_QT,ind_G_QT)));
tilda_G_QT(ind_G_QT,ind_G_QT) = G_QT(ind_G_QT,ind_G_QT)./((delta_G_QT)*(delta_G_QT'));
if ~isequal(rank(tilda_G_QT,temp.store_options_ident.tol_rank),rank(tilda_G_dynare,temp.store_options_ident.tol_rank))
error('ranks are not the same for normalized version')
end
max(max(abs(abs(tilda_G_dynare)-abs(tilda_G_QT))))
norm(tilda_G_dynare - tilda_G_QT)