From e37c2532ee116d163064daa0e845e712b7ac5167 Mon Sep 17 00:00:00 2001
From: Marco Ratto <marco.ratto@jrc.ec.europa.eu>
Date: Fri, 6 Jul 2012 12:19:25 +0200
Subject: [PATCH] * Fix problem with models where steadystate files change
parameter values.
1) allow to compute derivatives starting from NUMERICAL derivatives of jacobian and steady state: this has a minor cost in accuracy and allow apply without errors identification and estimation with numerical derivatives;
2) added trap in dynare_estimation_init: if steadystate changes param values, automaticly shifts to numerical derivs of jacoban and steady state + analytic derivatives of all the rest;
3) bug fixes for 2nd order derivatives w.r.t. model parameters;
**** (Manually cherry picked from ed4d37341c707ec80e8c3291a451e5e1388021f2)
---
matlab/dsge_likelihood.m | 8 ++-
matlab/dynare_estimation_init.m | 25 ++++++++--
matlab/dynare_identification.m | 2 +
matlab/getH.m | 83 ++++++++++++++++++++------------
matlab/getJJ.m | 13 ++++-
matlab/global_initialization.m | 1 +
matlab/identification_analysis.m | 8 ++-
matlab/thet2tau.m | 7 +++
8 files changed, 107 insertions(+), 40 deletions(-)
diff --git a/matlab/dsge_likelihood.m b/matlab/dsge_likelihood.m
index a646f82af4..da63a35cae 100644
--- a/matlab/dsge_likelihood.m
+++ b/matlab/dsge_likelihood.m
@@ -181,6 +181,10 @@ if nargout==1,
analytic_derivation=0;
end
+if analytic_derivation,
+ kron_flag=DynareOptions.analytic_derivation_mode;
+end
+
%------------------------------------------------------------------------------
% 1. Get the structural parameters & define penalties
%------------------------------------------------------------------------------
@@ -498,9 +502,9 @@ if analytic_derivation,
end
if full_Hess,
- [dum, DT, DOm, DYss, dum2, D2T, D2Om, D2Yss] = getH(A, B, Model,DynareResults,0,indparam,indexo);
+ [dum, DT, DOm, DYss, dum2, D2T, D2Om, D2Yss] = getH(A, B, Model,DynareResults,DynareOptions,kron_flag,indparam,indexo);
else
- [dum, DT, DOm, DYss] = getH(A, B, Model,DynareResults,0,indparam,indexo);
+ [dum, DT, DOm, DYss] = getH(A, B, Model,DynareResults,DynareOptions,kron_flag,indparam,indexo);
end
else
DT = derivatives_info.DT;
diff --git a/matlab/dynare_estimation_init.m b/matlab/dynare_estimation_init.m
index a66ed7821f..e1d42f0c13 100644
--- a/matlab/dynare_estimation_init.m
+++ b/matlab/dynare_estimation_init.m
@@ -285,10 +285,29 @@ else
end;
if options_.analytic_derivation,
- if ~(exist('sylvester3','file')==2),
+ if ~(exist('sylvester3','file')==2),
dynareroot = strrep(which('dynare'),'dynare.m','');
addpath([dynareroot 'gensylv'])
end
+ if estim_params_.np,
+ % check if steady state changes param values
+ M=M_;
+ M.params(estim_params_.param_vals(:,1)) = M.params(estim_params_.param_vals(:,1))*1.01;
+ if options_.diffuse_filter
+ steadystate_check_flag = 0;
+ else
+ steadystate_check_flag = 1;
+ end
+ [tmp1, params] = evaluate_steady_state(oo_.steady_state,M,options_,oo_,steadystate_check_flag);
+ change_flag=any(find(params-M.params));
+ if change_flag,
+ disp('The steadystate file changed the values for the following parameters: '),
+ disp(M.param_names(find(params-M.params),:))
+ disp('The derivatives of jacobian and steady-state will be computed numerically'),
+ disp('(re-set options_.analytic_derivation_mode= -2)'),
+ options_.analytic_derivation_mode= -2;
+ end
+ end
end
% Test if the data file is declared.
@@ -335,10 +354,10 @@ nvx = estim_params_.nvx;
ncx = estim_params_.ncx;
nvn = estim_params_.nvn;
ncn = estim_params_.ncn;
-if estim_params_.np
+if estim_params_.np,
M.params(estim_params_.param_vals(:,1)) = xparam1(nvx+ncx+nvn+ncn+1:end);
end;
-oo_.steady_state = evaluate_steady_state(oo_.steady_state,M,options_,oo_,steadystate_check_flag);
+[oo_.steady_state, params] = evaluate_steady_state(oo_.steady_state,M,options_,oo_,steadystate_check_flag);
if all(abs(oo_.steady_state(bayestopt_.mfys))<1e-9)
options_.noconstant = 1;
else
diff --git a/matlab/dynare_identification.m b/matlab/dynare_identification.m
index fb8174a855..46d75d1b8c 100644
--- a/matlab/dynare_identification.m
+++ b/matlab/dynare_identification.m
@@ -60,6 +60,7 @@ options_ident = set_default_option(options_ident,'replic',100);
options_ident = set_default_option(options_ident,'advanced',0);
options_ident = set_default_option(options_ident,'normalize_jacobians',1);
options_ident = set_default_option(options_ident,'lik_init',1);
+options_ident = set_default_option(options_ident,'analytic_derivation',1);
if options_ident.gsa_sample_file,
GSAFolder = checkpath('gsa',M_.dname);
if options_ident.gsa_sample_file==1,
@@ -117,6 +118,7 @@ options_ = set_default_option(options_,'datafile',[]);
options_.mode_compute = 0;
options_.plot_priors = 0;
[dataset_,xparam1, M_, options_, oo_, estim_params_,bayestopt_]=dynare_estimation_init(M_.endo_names,fname_,1, M_, options_, oo_, estim_params_, bayestopt_);
+options_ident.analytic_derivation_mode = options_.analytic_derivation_mode;
if isempty(dataset_),
dataset_.info.ntobs = periods;
dataset_.info.nvobs = rows(options_.varobs);
diff --git a/matlab/getH.m b/matlab/getH.m
index 34097cb58b..825d251382 100644
--- a/matlab/getH.m
+++ b/matlab/getH.m
@@ -1,4 +1,4 @@
-function [H, dA, dOm, Hss, gp, d2A, d2Om, H2ss] = getH(A, B, M_,oo_,kronflag,indx,indexo)
+function [H, dA, dOm, Hss, gp, d2A, d2Om, H2ss] = getH(A, B, M_,oo_,options_,kronflag,indx,indexo)
% computes derivative of reduced form linear model w.r.t. deep params
%
@@ -20,47 +20,73 @@ function [H, dA, dOm, Hss, gp, d2A, d2Om, H2ss] = getH(A, B, M_,oo_,kronflag,ind
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
if nargin<3 || isempty(kronflag), kronflag = 0; end
-if nargin<4 || isempty(indx), indx = [1:M_.param_nbr];, end,
-if nargin<5 || isempty(indexo), indexo = [];, end,
-
+if nargin<4 || isempty(indx), indx = [1:M_.param_nbr]; end,
+if nargin<5 || isempty(indexo), indexo = []; end,
[I,J]=find(M_.lead_lag_incidence');
yy0=oo_.dr.ys(I);
+param_nbr = length(indx);
+if nargout>5,
+ param_nbr_2 = param_nbr*(param_nbr+1)/2;
+end
+
+m = size(A,1);
+n = size(B,2);
+
+if kronflag==-2,
+ if nargout>5,
+ [residual, g1, g2 ] = feval([M_.fname,'_dynamic'],yy0, oo_.exo_steady_state', ...
+ M_.params, oo_.dr.ys, 1);
+ g22 = hessian('thet2tau',[M_.params(indx)],options_.gstep,M_, oo_, indx,[],-1);
+ H2ss=g22(1:M_.endo_nbr,:);
+ H2ss = reshape(H2ss,[M_.endo_nbr param_nbr param_nbr]);
+ g22=g22(M_.endo_nbr+1:end,:);
+ g22 = reshape(g22,[size(g1) param_nbr param_nbr]);
+ else
+ [residual, g1 ] = feval([M_.fname,'_dynamic'],yy0, oo_.exo_steady_state', ...
+ M_.params, oo_.dr.ys, 1);
+ end
+ gp = fjaco('thet2tau',[M_.params(indx)],M_, oo_, indx,[],-1);
+ Hss=gp(1:M_.endo_nbr,:);
+ gp=gp(M_.endo_nbr+1:end,:);
+ gp = reshape(gp,[size(g1) param_nbr]);
+else
+
% yy0=[];
% for j=1:size(M_.lead_lag_incidence,1);
% yy0 = [ yy0; oo_.dr.ys(find(M_.lead_lag_incidence(j,:)))];
% end
dyssdtheta=zeros(length(oo_.dr.ys),M_.param_nbr);
+d2yssdtheta=zeros(length(oo_.dr.ys),M_.param_nbr,M_.param_nbr);
if nargout>5,
[residual, gg1, gg2] = feval([M_.fname,'_static'],oo_.dr.ys, oo_.exo_steady_state', M_.params);
[df, gp, d2f] = feval([M_.fname,'_params_derivs'],yy0, oo_.exo_steady_state', ...
- M_.params, oo_.dr.ys, 1, dyssdtheta);
+ M_.params, oo_.dr.ys, 1, dyssdtheta, d2yssdtheta);
d2f = get_all_resid_2nd_derivs(d2f,length(oo_.dr.ys),M_.param_nbr);
d2f = d2f(:,indx,indx);
if isempty(find(gg2)),
for j=1:length(indx),
- d2yssdtheta(:,:,j) = -gg1\d2f(:,:,j);
+ d2yssdtheta(:,indx,j) = -gg1\d2f(:,:,j);
end
else
gam = d2f*0;
for j=1:length(indx),
- d2yssdtheta(:,:,j) = -gg1\(d2f(:,:,j)+gam(:,:,j));
+ d2yssdtheta(:,indx,j) = -gg1\(d2f(:,:,j)+gam(:,:,j));
end
end
else
[residual, gg1] = feval([M_.fname,'_static'],oo_.dr.ys, oo_.exo_steady_state', M_.params);
df = feval([M_.fname,'_params_derivs'],yy0, oo_.exo_steady_state', ...
- M_.params, oo_.dr.ys, 1, dyssdtheta);
+ M_.params, oo_.dr.ys, 1, dyssdtheta, d2yssdtheta);
end
dyssdtheta = -gg1\df;
if any(any(isnan(dyssdtheta))),
[U,T] = schur(gg1);
- global options_
qz_criterium=options_.qz_criterium;
e1 = abs(ordeig(T)) < qz_criterium-1;
k = sum(e1); % Number non stationary variables.
- n = length(e1)-k; % Number of stationary variables.
+% n = length(e1)-k; % Number of stationary variables.
[U,T] = ordschur(U,T,e1);
T = T(k+1:end,k+1:end);
dyssdtheta = -U(:,k+1:end)*(T\U(:,k+1:end)')*df;
@@ -71,22 +97,24 @@ if any(any(isnan(dyssdtheta))),
end
end
if nargout>5,
-[df, gp, d2f, gpp] = feval([M_.fname,'_params_derivs'],yy0, oo_.exo_steady_state', ...
- M_.params, oo_.dr.ys, 1, dyssdtheta);
-H2ss = d2yssdtheta(oo_.dr.order_var,:,:);
-[residual, g1, g2, g3] = feval([M_.fname,'_dynamic'],yy0, oo_.exo_steady_state', ...
- M_.params, oo_.dr.ys, 1);
+ [df, gp, d2f, gpp] = feval([M_.fname,'_params_derivs'],yy0, oo_.exo_steady_state', ...
+ M_.params, oo_.dr.ys, 1, dyssdtheta, d2yssdtheta);
+ H2ss = d2yssdtheta(oo_.dr.order_var,:,:);
+ [residual, g1, g2, g3] = feval([M_.fname,'_dynamic'],yy0, oo_.exo_steady_state', ...
+ M_.params, oo_.dr.ys, 1);
+ nelem=size(g1,2);
+ g22 = get_all_2nd_derivs(gpp,m,nelem,M_.param_nbr);
+ g22 = g22(:,:,indx,indx);
else
-[df, gp] = feval([M_.fname,'_params_derivs'],yy0, oo_.exo_steady_state', ...
- M_.params, oo_.dr.ys, 1, dyssdtheta);
-[residual, g1, g2 ] = feval([M_.fname,'_dynamic'],yy0, oo_.exo_steady_state', ...
- M_.params, oo_.dr.ys, 1);
+ [df, gp] = feval([M_.fname,'_params_derivs'],yy0, oo_.exo_steady_state', ...
+ M_.params, oo_.dr.ys, 1, dyssdtheta,d2yssdtheta);
+ [residual, g1, g2 ] = feval([M_.fname,'_dynamic'],yy0, oo_.exo_steady_state', ...
+ M_.params, oo_.dr.ys, 1);
end
Hss = dyssdtheta(oo_.dr.order_var,indx);
dyssdtheta = dyssdtheta(I,:);
[nr, nc]=size(g2);
-[m, nelem]=size(g1);
nc = sqrt(nc);
ns = max(max(M_.lead_lag_incidence));
gp2 = gp*0;
@@ -105,14 +133,8 @@ end
gp = gp+gp2;
gp = gp(:,:,indx);
-param_nbr = length(indx);
-if nargout>5,
- param_nbr_2 = param_nbr*(param_nbr+1)/2;
end
-m = size(A,1);
-n = size(B,2);
-
klen = M_.maximum_endo_lag + M_.maximum_endo_lead + 1;
@@ -120,9 +142,6 @@ k11 = M_.lead_lag_incidence(find([1:klen] ~= M_.maximum_endo_lag+1),:);
a = g1(:,nonzeros(k11'));
da = gp(:,nonzeros(k11'),:);
if nargout > 5,
- nelem = size(g1,2);
- g22 = get_all_2nd_derivs(gpp,m,nelem,M_.param_nbr);
- g22 = g22(:,:,indx,indx);
d2a = g22(:,nonzeros(k11'),:,:);
end
kstate = oo_.dr.kstate;
@@ -250,8 +269,8 @@ end
% B0=B;
% B = Bx; clear Bx B1;
-m = size(A,1);
-n = size(B,2);
+% m = size(A,1);
+% n = size(B,2);
% Dg1 = zeros(m,m,param_nbr);
% Dg1(:, nf, :) = -gp(:,M_.lead_lag_incidence(3,nf),:);
@@ -499,7 +518,7 @@ is=find(gpp(:,3)==i);
is=is(find(gpp(is,4)==j));
if ~isempty(is),
- g22(gpp(is,1),gpp(is,2))=gpp(is,5);
+ g22(sub2ind([m,n],gpp(is,1),gpp(is,2)))=gpp(is,5)';
end
return
diff --git a/matlab/getJJ.m b/matlab/getJJ.m
index 859edf0dab..b717858a59 100644
--- a/matlab/getJJ.m
+++ b/matlab/getJJ.m
@@ -33,10 +33,21 @@ if kronflag == -1,
params0 = M_.params;
H = fjaco(fun,[sqrt(diag(M_.Sigma_e(indexo,indexo))); M_.params(indx)],M_, oo_, indx,indexo,0,mf,nlags,useautocorr);
M_.params = params0;
+ params0 = M_.params;
+ gp = fjaco(fun,[sqrt(diag(M_.Sigma_e(indexo,indexo))); M_.params(indx)],M_, oo_, indx,indexo,-1);
+ M_.params = params0;
+ offset = length(indexo);
+ gp = gp(1:M_.endo_nbr,offset+1:end);
+ dYss = H(1:M_.endo_nbr,offset+1:end);
+ dA = reshape(H(M_.orig_endo_nbr+[1:numel(A)],:),[size(A),size(H,2)]);
+ dOm = dA*0;
+ for j=1:size(H,2),
+ dOm(:,:,j) = dyn_unvech(H(M_.orig_endo_nbr+numel(A)+1:end,j));
+ end
assignin('base','M_', M_);
assignin('base','oo_', oo_);
else
- [H, dA, dOm, dYss, gp] = getH(A, B, M_,oo_,kronflag,indx,indexo);
+ [H, dA, dOm, dYss, gp] = getH(A, B, M_,oo_,options_,kronflag,indx,indexo);
gp = reshape(gp,size(gp,1)*size(gp,2),size(gp,3));
gp = [dYss; gp];
% if isempty(H),
diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index aa2f677039..4ab1d3bdc6 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -320,6 +320,7 @@ options_.MaxNumberOfBytes = 1e6;
options_.MaximumNumberOfMegaBytes = 111;
options_.PosteriorSampleSize = 1000;
options_.analytic_derivation = 0;
+options_.analytic_derivation_mode = 0;
options_.bayesian_irf = 0;
options_.bayesian_th_moments = 0;
options_.diffuse_filter = 0;
diff --git a/matlab/identification_analysis.m b/matlab/identification_analysis.m
index 12862b7df2..2c0b701f68 100644
--- a/matlab/identification_analysis.m
+++ b/matlab/identification_analysis.m
@@ -57,6 +57,8 @@ replic = options_ident.replic;
periods = options_ident.periods;
max_dim_cova_group = options_ident.max_dim_cova_group;
normalize_jacobians = options_ident.normalize_jacobians;
+kron_flag = options_ident.analytic_derivation_mode;
+
[I,J]=find(M_.lead_lag_incidence');
ide_hess = struct();
@@ -75,7 +77,7 @@ if info(1)==0,
oo_.dr.ys, 1);
vg1 = [oo_.dr.ys(oo_.dr.order_var); vec(g1)];
- [JJ, H, gam, gp, dA, dOm, dYss] = getJJ(A, B, M_,oo0,options_,0,indx,indexo,bayestopt_.mf2,nlags,useautocorr);
+ [JJ, H, gam, gp, dA, dOm, dYss] = getJJ(A, B, M_,oo0,options_,kron_flag,indx,indexo,bayestopt_.mf2,nlags,useautocorr);
derivatives_info.DT=dA;
derivatives_info.DOm=dOm;
derivatives_info.DYss=dYss;
@@ -85,7 +87,7 @@ if info(1)==0,
disp('The number of moments with non-zero derivative is smaller than the number of parameters')
disp(['Try increasing ar = ', int2str(nlags+1)])
nlags=nlags+1;
- [JJ, H, gam, gp, dA, dOm, dYss] = getJJ(A, B, M_,oo0,options_,0,indx,indexo,bayestopt_.mf2,nlags,useautocorr);
+ [JJ, H, gam, gp, dA, dOm, dYss] = getJJ(A, B, M_,oo0,options_,kron_flag,indx,indexo,bayestopt_.mf2,nlags,useautocorr);
derivatives_info.DT=dA;
derivatives_info.DOm=dOm;
derivatives_info.DYss=dYss;
@@ -133,6 +135,7 @@ if info(1)==0,
if options_.kalman_algo > 2,
options_.kalman_algo = 1;
end
+ analytic_derivation = options_.analytic_derivation;
options_.analytic_derivation = -2;
info = stoch_simul(options_.varobs);
data_info.data=oo_.endo_simul(options_.varobs_id,100+1:end);
@@ -140,6 +143,7 @@ if info(1)==0,
derivatives_info.no_DLIK=1;
[fval,DLIK,AHess,cost_flag,ys,trend_coeff,info,M_,options_,bayestopt_,oo_] = dsge_likelihood(params',data_info,options_,M_,estim_params_,bayestopt_,oo_,derivatives_info);
% fval = DsgeLikelihood(xparam1,data_info,options_,M_,estim_params_,bayestopt_,oo_);
+ options_.analytic_derivation = analytic_derivation;
AHess=-AHess;
if min(eig(AHess))<0,
error('Analytic Hessian is not positive semi-definite!')
diff --git a/matlab/thet2tau.m b/matlab/thet2tau.m
index a59b374b8d..79c593fad8 100644
--- a/matlab/thet2tau.m
+++ b/matlab/thet2tau.m
@@ -39,6 +39,13 @@ end
[A,B,tele,tubbies,M_,options_,oo_] = dynare_resolve(M_,options_,oo_);
if flagmoments==0,
tau = [oo_.dr.ys(oo_.dr.order_var); A(:); dyn_vech(B*M_.Sigma_e*B')];
+elseif flagmoments==-1
+ [I,J]=find(M_.lead_lag_incidence');
+ yy0=oo_.dr.ys(I);
+ [residual, g1] = feval([M_.fname,'_dynamic'],yy0, oo_.exo_steady_state', ...
+ M_.params, oo_.dr.ys, 1);
+ tau=[oo_.dr.ys(oo_.dr.order_var); g1(:)];
+
else
GAM = lyapunov_symm(A,B*M_.Sigma_e*B',options_.qz_criterium,options_.lyapunov_complex_threshold);
k = find(abs(GAM) < 1e-12);
--
GitLab