stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
dynare/tests/analytic_derivatives/BrockMirman_PertParamsDeriv...

474 lines
23 KiB
Modula-2
Raw Permalink Blame History

% "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 © 2019-2020 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 <https://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, estim_params_]= 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
[oo_.dr,info,M_.params] = resol(0,M_, options_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state);
%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 = identification.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 = identification.get_perturbation_params_derivs(M_, options_, estim_params_, oo_.dr, oo_.steady_state, oo_.exo_steady_state, oo_.exo_det_steady_state, 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
Reference in New Issue View Git Blame Copy Permalink