Verified Commit 46c4dea5 authored by Willi Mutschler's avatar Willi Mutschler
Browse files

📄 Updated code comments

parent a62e69cf
This diff is collapsed.
......@@ -10,7 +10,7 @@ function fjac = fjaco(f,x,varargin)
% OUTPUT
% fjac : finite difference Jacobian
%
% Copyright (C) 2010-2017,2019 Dynare Team
% Copyright (C) 2010-2017,2019-2020 Dynare Team
%
% This file is part of Dynare.
%
......
function [E_y, dE_y, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params, M, oo, options, options_ident, indpmodel, indpstderr, indpcorr, indvobs)
function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params, M, oo, options, options_ident, indpmodel, indpstderr, indpcorr, indvobs)
% function [MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params, M, oo, options, options_ident, indpmodel, indpstderr, indpcorr, indvobs)
% previously getJJ.m
% % Sets up the Jacobians needed for identification analysis
% previously getJJ.m in Dynare 4.5
% Sets up the Jacobians needed for identification analysis
% =========================================================================
% INPUTS
% estim_params: [structure] storing the estimation information
......@@ -27,46 +27,50 @@ function [E_y, dE_y, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOM
%
% MEAN [endo_nbr by 1], in DR order. Expectation of all model variables
% * order==1: corresponds to steady state
% * order==2: computed from pruned state space system (as in e.g. Andreasen, Fernandez-Villaverde, Rubio-Ramirez, 2018)
% * order==2|3: corresponds to mean computed from pruned state space system (as in Andreasen, Fernandez-Villaverde, Rubio-Ramirez, 2018)
% dMEAN [endo_nbr by totparam_nbr], in DR Order, Jacobian (wrt all params) of MEAN
%
% REDUCEDFORM [rowredform_nbr by 1] in DR order. Steady state and reduced-form model solution matrices for all model variables
% * order==1: [Yss' vec(dghx)' vech(dghu*Sigma_e*dghu')']',
% where rowredform_nbr = endo_nbr+endo_nbr*nspred+endo_nbr*(endo_nbr+1)/2
% * order==2: [Yss' vec(dghx)' vech(dghu*Sigma_e*dghu')' vec(dghxx)' vec(dghxu)' vec(dghuu)' vec(dghs2)']',
% where rowredform_nbr = endo_nbr+endo_nbr*nspred+endo_nbr*(endo_nbr+1)/2+endo_nbr*nspred^2*endo_nbr*nspred*exo_nr+endo_nbr*exo_nbr^2+endo_nbr
% dREDUCEDFORM: [rowrf_nbr by totparam_nbr] in DR order, Jacobian (wrt all params) of REDUCEDFORM
% * order==1: [Yss' vec(ghx)' vech(ghu*Sigma_e*ghu')']',
% where rowredform_nbr = endo_nbr*(1+nspred+(endo_nbr+1)/2)
% * order==2: [Yss' vec(ghx)' vech(ghu*Sigma_e*ghu')' vec(ghxx)' vec(ghxu)' vec(ghuu)' vec(ghs2)']',
% where rowredform_nbr = endo_nbr*(1+nspred+(endo_nbr+1)/2+nspred^2+nspred*exo_nr+exo_nbr^2+1)
% * order==3: [Yss' vec(ghx)' vech(ghu*Sigma_e*ghu')' vec(ghxx)' vec(ghxu)' vec(ghuu)' vec(ghs2)' vec(ghxxx)' vec(ghxxu)' vec(ghxuu)' vec(ghuuu)' vec(ghxss)' vec(ghuss)']',
% where rowredform_nbr = endo_nbr*(1+nspred+(endo_nbr+1)/2+nspred^2+nspred*exo_nr+exo_nbr^2+1+nspred^3+nspred^2*exo_nbr+nspred*exo_nbr^2+exo_nbr^3+nspred+exo_nbr)
% dREDUCEDFORM: [rowredform_nbr by totparam_nbr] in DR order, Jacobian (wrt all params) of REDUCEDFORM
% * order==1: corresponds to Iskrev (2010)'s J_2 matrix
% * order==2: corresponds to Mutschler (2015)'s J matrix
%
% DYNAMIC [rowdyn_nbr by 1] in declaration order. Steady state and dynamic model derivatives for all model variables
% * order==1: [ys' vec(g1)']', rowdyn_nbr=endo_nbr+length(g1)
% * order==2: [ys' vec(g1)' vec(g2)']', rowdyn_nbr=endo_nbr+length(g1)+length(g2)
% * order==3: [ys' vec(g1)' vec(g2)' vec(g3)']', rowdyn_nbr=endo_nbr+length(g1)+length(g2)+length(g3)
% dDYNAMIC [rowdyn_nbr by modparam_nbr] in declaration order. Jacobian (wrt model parameters) of DYNAMIC
% * order==1: corresponds to Ratto and Iskrev (2011)'s J_\Gamma matrix (or LRE)
% * order==2: additionally includes derivatives of dynamic hessian
%
% MOMENTS: [obs_nbr+obs_nbr*(obs_nbr+1)/2+nlags*obs_nbr^2 by 1] in DR order. First two theoretical moments for VAROBS variables, i.e.
% [E[yobs]' vech(E[yobs*yobs'])' vec(E[yobs*yobs(-1)'])' ... vec(E[yobs*yobs(-nlag)'])']
% [E[varobs]' vech(E[varobs*varobs'])' vec(E[varobs*varobs(-1)'])' ... vec(E[varobs*varobs(-nlag)'])']
% dMOMENTS: [obs_nbr+obs_nbr*(obs_nbr+1)/2+nlags*obs_nbr^2 by totparam_nbr] in DR order. Jacobian (wrt all params) of MOMENTS
% * order==1: corresponds to Iskrev (2010)'s J matrix
% * order==2: corresponds to Mutschler (2015)'s \bar{M}_2 matrix, i.e. theoretical moments from the pruned state space system
%
% dSPECTRUM: [totparam_nbr by totparam_nbr] in DR order. Gram matrix of Jacobian (wrt all params) of spectral density for VAROBS variables, where
% spectral density at frequency w: f(w) = (2*pi)^(-1)*H(exp(-i*w))*E[xi*xi']*ctranspose(H(exp(-i*w)) with H being the Transfer function
% spectral density at frequency w: f(w) = (2*pi)^(-1)*H(exp(-i*w))*E[Inov*Inov']*ctranspose(H(exp(-i*w)) with H being the Transfer function
% dSPECTRUM = int_{-\pi}^\pi transpose(df(w)/dp')*(df(w)/dp') dw
% * order==1: corresponds to Qu and Tkachenko (2012)'s G matrix, where xi and H are computed from linear state space system
% * order==2: corresponds to Mutschler (2015)'s G_2 matrix, where xi and H are computed from pruned state space system
% * order==1: corresponds to Qu and Tkachenko (2012)'s G matrix, where Inov and H are computed from linear state space system
% * order==2: corresponds to Mutschler (2015)'s G_2 matrix, where Inov and H are computed from second-order pruned state space system
% * order==3: Inov and H are computed from third-order pruned state space system
%
% dMINIMAL: [obs_nbr+minx^2+minx*exo_nbr+obs_nbr*minx+obs_nbr*exo_nbr+exo_nbr*(exo_nbr+1)/2 by totparam_nbr+minx^2+exo_nbr^2]
% dMINIMAL: [obs_nbr+minx_nbr^2+minx_nbr*exo_nbr+obs_nbr*minx_nbr+obs_nbr*exo_nbr+exo_nbr*(exo_nbr+1)/2 by totparam_nbr+minx_nbr^2+exo_nbr^2]
% Jacobian (wrt all params, and similarity_transformation_matrices (T and U)) of observational equivalent minimal ABCD system,
% corresponds to Komunjer and Ng (2011)'s Deltabar matrix, where
% MINIMAL = [vec(E[yobs]' vec(minA)' vec(minB)' vec(minC)' vec(minD)' vech(Sigma_e)']'
% * order==1: E[yobs] is equal to steady state
% * order==2: E[yobs] is computed from the pruned state space system (second-order accurate), as noted in section 5 of Komunjer and Ng (2011)
% MINIMAL = [vec(E[varobs]' vec(minA)' vec(minB)' vec(minC)' vec(minD)' vech(Sigma_e)']'
% minA, minB, minC and minD is the minimal state space system computed in get_minimal_state_representation
% * order==1: E[varobs] is equal to steady state
% * order==2|3: E[varobs] is computed from the pruned state space system (second|third-order accurate), as noted in section 5 of Komunjer and Ng (2011)
%
% derivatives_info [structure] for use in dsge_likelihood to compute Hessian analytically.
% Contains dA, dB, and d(B*Sigma_e*B'), where A and B are from Kalman filter. Only used at order==1.
% derivatives_info [structure] for use in dsge_likelihood to compute Hessian analytically. Only used at order==1.
% Contains dA, dB, and d(B*Sigma_e*B'), where A and B are Kalman filter transition matrice.
%
% -------------------------------------------------------------------------
% This function is called by
......@@ -220,7 +224,7 @@ options.options_ident.indpmodel = indpmodel;
options.options_ident.indpstderr = indpstderr;
options.options_ident.indpcorr = indpcorr;
oo.dr = pruned_state_space_system(M, options, oo.dr);
E_y = oo.dr.pruned.E_y; dE_y = oo.dr.pruned.dE_y;
MEAN = oo.dr.pruned.E_y; dMEAN = oo.dr.pruned.dE_y;
A = oo.dr.pruned.A; dA = oo.dr.pruned.dA;
B = oo.dr.pruned.B; dB = oo.dr.pruned.dB;
C = oo.dr.pruned.C; dC = oo.dr.pruned.dC;
......@@ -247,17 +251,17 @@ end
% yhat = C*zhat(-1) + D*xi, where yhat = y - E(y)
if ~no_identification_moments
MOMENTS = identification_numerical_objective(para0, 1, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvobs, useautocorr, nlags, grid_nbr); %[outputflag=1]
MOMENTS = [E_y; MOMENTS];
MOMENTS = [MEAN; MOMENTS];
if kronflag == -1
%numerical derivative of autocovariogram
dMOMENTS = fjaco(str2func('identification_numerical_objective'), para0, 1, estim_params, M, oo, options, indpmodel, indpstderr, indpcorr, indvobs, useautocorr, nlags, grid_nbr); %[outputflag=1]
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
dMOMENTS = [dE_y; dMOMENTS]; %add Jacobian of steady state of VAROBS variables
dMOMENTS = [dMEAN; dMOMENTS]; %add Jacobian of steady state of VAROBS variables
else
dMOMENTS = zeros(obs_nbr + obs_nbr*(obs_nbr+1)/2 + nlags*obs_nbr^2 , totparam_nbr);
dMOMENTS(1:obs_nbr,:) = dE_y; %add Jacobian of first moments of VAROBS variables
dMOMENTS(1:obs_nbr,:) = dMEAN; %add Jacobian of first moments of VAROBS variables
% Denote Ezz0 = E[zhat*zhat'], then the following Lyapunov equation defines the autocovariagram: Ezz0 -A*Ezz0*A' = B*Sig_xi*B' = Om_z
[Ezz0,u] = lyapunov_symm(A, Om_z, options.lyapunov_fixed_point_tol, options.qz_criterium, options.lyapunov_complex_threshold, 1, options.debug);
stationary_vars = (1:length(indvobs))';
......@@ -426,7 +430,7 @@ if ~no_identification_spectrum
end
end
% Normalize Matrix and add steady state Jacobian, note that G is real and symmetric by construction
dSPECTRUM = 2*pi*dSPECTRUM./(2*length(freqs)-1) + dE_y'*dE_y;
dSPECTRUM = 2*pi*dSPECTRUM./(2*length(freqs)-1) + dMEAN'*dMEAN;
dSPECTRUM = real(dSPECTRUM);
else
dSPECTRUM = [];
......@@ -482,7 +486,7 @@ if ~no_identification_minimal
-2*Enu*kron(Sigma_e0,Inu)];
dMINIMAL = full([KomunjerNg_DL KomunjerNg_DT KomunjerNg_DU]);
%add Jacobian of steady state (here we also allow for higher-order perturbation, i.e. only the mean provides additional restrictions
dMINIMAL = [dE_y zeros(obs_nbr,minnx^2+exo_nbr^2); dMINIMAL];
dMINIMAL = [dMEAN zeros(obs_nbr,minnx^2+exo_nbr^2); dMINIMAL];
end
end
else
......
......@@ -13,22 +13,25 @@ function [out,info] = get_perturbation_params_derivs_numerical_objective(params,
% 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
% OUTPUT
% out (dependent on outputflag and order of approximation):
% - 'perturbation_solution': out = out1 = [vec(Sigma_e);vec(ghx);vec(ghu)]; (order==1)
% out = out2 = [out1;vec(ghxx);vec(ghxu);vec(ghuu);vec(ghs2)]; (order==2)
% out = out3 = [out1;out2;vec(ghxxx);vec(ghxxu);vec(ghxuu);vec(ghuuu);vec(ghxss);vec(ghuss)]; (order==3)
% - '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
% info [integer] output from resol
% -------------------------------------------------------------------------
% This function is called by
% * get_solution_params_deriv.m (previously getH.m)
% * get_identification_jacobians.m (previously getJJ.m)
% * get_perturbation_params_derivs.m (previously getH.m)
% -------------------------------------------------------------------------
% This function calls
% * [M.fname,'.dynamic']
% * resol
% * dyn_vech
% =========================================================================
% Copyright (C) 2019 Dynare Team
% Copyright (C) 2019-2020 Dynare Team
%
% This file is part of Dynare.
%
......@@ -48,8 +51,8 @@ function [out,info] = get_perturbation_params_derivs_numerical_objective(params,
%% 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);
Sigma_e = M.Sigma_e;
if info(1) > 0
% there are errors in the solution algorithm
......@@ -60,11 +63,9 @@ else
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
......
......@@ -28,7 +28,7 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide
% 1: prior exists. Identification is checked for stderr, corr and model parameters as declared in estimated_params block
% 0: prior does not exist. Identification is checked for all stderr and model parameters, but no corr parameters
% * init [integer]
% flag for initialization of persistent vars. This is needed if one may wants to make more calls to identification in the same dynare mod file
% flag for initialization of persistent vars. This is needed if one may want to make more calls to identification in the same mod file
% -------------------------------------------------------------------------
% OUTPUTS
% * ide_moments [structure]
......@@ -36,7 +36,7 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide
% * ide_spectrum [structure]
% identification results using spectrum (Qu and Tkachenko, 2012; Mutschler, 2015)
% * ide_minimal [structure]
% identification results using minimal system (Komunjer and Ng, 2011)
% identification results using theoretical mean and minimal system (Komunjer and Ng, 2011)
% * ide_hess [structure]
% identification results using asymptotic Hessian (Ratto and Iskrev, 2011)
% * ide_reducedform [structure]
......@@ -44,7 +44,7 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide
% * ide_dynamic [structure]
% identification results using steady state and dynamic model derivatives (Ratto and Iskrev, 2011)
% * derivatives_info [structure]
% info about derivatives, used in dsge_likelihood.m
% info about first-order perturbation derivatives, used in dsge_likelihood.m
% * info [integer]
% output from dynare_resolve
% * options_ident [structure]
......@@ -71,7 +71,7 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide
% * stoch_simul
% * vec
% =========================================================================
% Copyright (C) 2008-2019 Dynare Team
% Copyright (C) 2008-2020 Dynare Team
%
% This file is part of Dynare.
%
......@@ -91,16 +91,17 @@ function [ide_moments, ide_spectrum, ide_minimal, ide_hess, ide_reducedform, ide
global oo_ M_ options_ bayestopt_ estim_params_
persistent ind_dMOMENTS ind_dREDUCEDFORM ind_dDYNAMIC
% persistence is necessary, because in a MC loop the numerical threshold used may provide vectors of different length, leading to crashes in MC loops
% persistent indices are necessary, because in a MC loop the numerical threshold
% used may provide vectors of different length, leading to crashes in MC loops
%initialize output structures
ide_hess = struct(); %Jacobian of asymptotic/simulated information matrix
ide_reducedform = struct(); %Jacobian of steady state and reduced form solution
ide_dynamic = struct(); %Jacobian of steady state and dynamic model derivatives
ide_moments = struct(); %Jacobian of first two moments (Iskrev, 2010; Mutschler, 2015)
ide_spectrum = struct(); %Gram matrix of Jacobian of spectral density plus Gram matrix of Jacobian of steady state (Qu and Tkachenko, 2012; Mutschler, 2015)
ide_minimal = struct(); %Jacobian of minimal system (Komunjer and Ng, 2011)
derivatives_info = struct(); %storage for Jacobians used in dsge_likelihood.m for (analytical or simulated) asymptotic Gradient and Hession of likelihood
ide_hess = struct(); %Identification structure based on asymptotic/simulated information matrix
ide_reducedform = struct(); %Identification structure based on steady state and reduced form solution
ide_dynamic = struct(); %Identification structure based on steady state and dynamic model derivatives
ide_moments = struct(); %Identification structure based on first two moments (Iskrev, 2010; Mutschler, 2015)
ide_spectrum = struct(); %Identification structure based on Gram matrix of Jacobian of spectral density plus Gram matrix of Jacobian of steady state (Qu and Tkachenko, 2012; Mutschler, 2015)
ide_minimal = struct(); %Identification structure based on mean and minimal system (Komunjer and Ng, 2011)
derivatives_info = struct(); %storage for first-order perturbation Jacobians used in dsge_likelihood.m
totparam_nbr = length(params); %number of all parameters to be checked
modparam_nbr = length(indpmodel); %number of model parameters to be checked
......@@ -114,11 +115,6 @@ if ~isempty(estim_params_)
end
%get options (see dynare_identification.m for description of options)
no_identification_strength = options_ident.no_identification_strength;
no_identification_reducedform = options_ident.no_identification_reducedform;
no_identification_moments = options_ident.no_identification_moments;
no_identification_minimal = options_ident.no_identification_minimal;
no_identification_spectrum = options_ident.no_identification_spectrum;
order = options_ident.order;
nlags = options_ident.ar;
advanced = options_ident.advanced;
......@@ -130,11 +126,17 @@ checks_via_subsets = options_ident.checks_via_subsets;
tol_deriv = options_ident.tol_deriv;
tol_rank = options_ident.tol_rank;
tol_sv = options_ident.tol_sv;
no_identification_strength = options_ident.no_identification_strength;
no_identification_reducedform = options_ident.no_identification_reducedform;
no_identification_moments = options_ident.no_identification_moments;
no_identification_minimal = options_ident.no_identification_minimal;
no_identification_spectrum = options_ident.no_identification_spectrum;
%Compute linear approximation and fill dr structure
[oo_.dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
if info(1) == 0 %no errors in solution
% Compute parameter Jacobians for identification analysis
[MEAN, dMEAN, REDUCEDFORM, dREDUCEDFORM, DYNAMIC, dDYNAMIC, MOMENTS, dMOMENTS, dSPECTRUM, dMINIMAL, derivatives_info] = get_identification_jacobians(estim_params_, M_, oo_, options_, options_ident, indpmodel, indpstderr, indpcorr, indvobs);
if isempty(dMINIMAL)
% Komunjer and Ng is not computed if (1) minimality conditions are not fullfilled or (2) there are more shocks and measurement errors than observables, so we need to reset options
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment