Commit ed4d3734 authored by Marco Ratto's avatar Marco Ratto
Browse files

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;
parent 379972d7
......@@ -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
%------------------------------------------------------------------------------
......@@ -490,9 +494,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;
......
......@@ -191,7 +191,7 @@ end
bayestopt_.penalty = 1e8;
% Get informations about the variables of the model.
dr = set_state_space(oo_.dr,M_,options_);
dr = set_state_space(oo_.dr,M_);
oo_.dr = dr;
nstatic = dr.nstatic; % Number of static variables.
npred = dr.npred; % Number of predetermined variables.
......@@ -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
......
......@@ -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);
......
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
......
......@@ -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),
......
......@@ -315,6 +315,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;
......
......@@ -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!')
......
......@@ -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);
......
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