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

Merge pull request #1191 from JohannesPfeifer/identification 

Identification fixes
time-shift
Stéphane Adjemian 2016-05-19 21:34:30 +02:00
parent 4ea52066b9 89fd14f026
commit 5f8d204c02
3 changed files with 87 additions and 24 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
matlab/dynare_identification.m
View File

@ -38,6 +38,8 @@ function [pdraws, TAU, GAM, LRE, gp, H, JJ] = dynare_identification(options_iden
global M_ options_ oo_ bayestopt_ estim_params_
options0_ = options_;
if isoctave
warning('off'),
else
@ -280,7 +282,7 @@ if iload <=0,
parameters_TeX = 'Current parameter values';
disp('Testing current parameter values')
end
[idehess_point, idemoments_point, idemodel_point, idelre_point, derivatives_info_point, info] = ...
[idehess_point, idemoments_point, idemodel_point, idelre_point, derivatives_info_point, info, options_ident] = ...
identification_analysis(params,indx,indexo,options_ident,dataset_, dataset_info, prior_exist, name_tex,1);
if info(1)~=0,
skipline()
@ -323,7 +325,7 @@ if iload <=0,
while kk<50 && info(1),
kk=kk+1;
params = prior_draw();
[idehess_point, idemoments_point, idemodel_point, idelre_point, derivatives_info_point, info] = ...
[idehess_point, idemoments_point, idemodel_point, idelre_point, derivatives_info_point, info, options_ident] = ...
identification_analysis(params,indx,indexo,options_ident,dataset_,dataset_info, prior_exist, name_tex,1);
end
end
@ -376,7 +378,7 @@ if iload <=0,
else
params = prior_draw();
end
[dum1, ideJ, ideH, ideGP, dum2 , info] = ...
[dum1, ideJ, ideH, ideGP, dum2 , info, options_MC] = ...
identification_analysis(params,indx,indexo,options_MC,dataset_, dataset_info, prior_exist, name_tex,0);
if iteration==0 && info(1)==0,
MAX_tau = min(SampleSize,ceil(MaxNumberOfBytes/(size(ideH.siH,1)*nparam)/8));
@ -547,7 +549,7 @@ if SampleSize > 1,
fprintf('\n')
disp(['Testing ',tittxt, '. Press ENTER']), pause(5),
if ~iload,
[idehess_max, idemoments_max, idemodel_max, idelre_max, derivatives_info_max] = ...
[idehess_max, idemoments_max, idemodel_max, idelre_max, derivatives_info_max, info_max, options_ident] = ...
identification_analysis(pdraws(jmax,:),indx,indexo,options_ident,dataset_,dataset_info, prior_exist, name_tex,1);
save([IdentifDirectoryName '/' M_.fname '_identif.mat'], 'idehess_max', 'idemoments_max','idemodel_max', 'idelre_max', 'jmax', '-append');
end
@ -562,7 +564,7 @@ if SampleSize > 1,
fprintf('\n')
disp(['Testing ',tittxt, '. Press ENTER']), pause(5),
if ~iload,
[idehess_min, idemoments_min, idemodel_min, idelre_min, derivatives_info_min] = ...
[idehess_min, idemoments_min, idemodel_min, idelre_min, derivatives_info_min, info_min, options_ident] = ...
identification_analysis(pdraws(jmin,:),indx,indexo,options_ident,dataset_, dataset_info, prior_exist, name_tex,1);
save([IdentifDirectoryName '/' M_.fname '_identif.mat'], 'idehess_min', 'idemoments_min','idemodel_min', 'idelre_min', 'jmin', '-append');
end
@ -577,7 +579,7 @@ if SampleSize > 1,
fprintf('\n')
disp(['Testing ',tittxt, '. Press ENTER']), pause(5),
if ~iload,
[idehess_(j), idemoments_(j), idemodel_(j), idelre_(j), derivatives_info_(j)] = ...
[idehess_(j), idemoments_(j), idemodel_(j), idelre_(j), derivatives_info_(j), info_resolve, options_ident] = ...
identification_analysis(pdraws(jcrit(j),:),indx,indexo,options_ident,dataset_, dataset_info, prior_exist, name_tex,1);
end
disp_identification(pdraws(jcrit(j),:), idemodel_(j), idemoments_(j), name,1);
@ -603,3 +605,5 @@ end
skipline()
disp(['==== Identification analysis completed ====' ]),
skipline(2)
options_ = options0_;

12
matlab/identification_analysis.m
View File

@ -1,4 +1,4 @@
function [ide_hess, ide_moments, ide_model, ide_lre, derivatives_info, info] = identification_analysis(params,indx,indexo,options_ident,dataset_,dataset_info, prior_exist, name_tex, init)
function [ide_hess, ide_moments, ide_model, ide_lre, derivatives_info, info, options_ident] = identification_analysis(params,indx,indexo,options_ident,dataset_,dataset_info, prior_exist, name_tex, init)
% function [ide_hess, ide_moments, ide_model, ide_lre, derivatives_info, info] = identification_analysis(params,indx,indexo,options_ident,data_info, prior_exist, name_tex, init)
% given the parameter vector params, wraps all identification analyses
%
@ -88,6 +88,8 @@ if info(1)==0,
indJJ = (find(max(abs(JJ'),[],1)>1.e-8));
if isempty(indJJ) && any(any(isnan(JJ)))
error('There are NaN in the JJ matrix. Please check whether your model has units roots and you forgot to set diffuse_filter=1.' )
elseif any(any(isnan(gam)))
error('There are NaN''s in the theoretical moments: make sure that for non-stationary models stationary transformations of non-stationary observables are used for checking identification. [TIP: use first differences].')
end
while length(indJJ)<nparam && nlags<10,
disp('The number of moments with non-zero derivative is smaller than the number of parameters')
@ -97,14 +99,14 @@ if info(1)==0,
derivatives_info.DT=dA;
derivatives_info.DOm=dOm;
derivatives_info.DYss=dYss;
evalin('caller',['options_ident.ar=',int2str(nlags),';']);
options_.ar=nlags;
indJJ = (find(max(abs(JJ'),[],1)>1.e-8));
end
if length(indJJ)<nparam,
disp('The number of moments with non-zero derivative is smaller than the number of parameters')
disp('up to 10 lags: check your model')
disp('Either further increase ar or reduce the list of estimated parameters')
error('IDETooManyParams',''),
error('identification_analysis: there are not enough moments and too many parameters'),
end
indH = (find(max(abs(H'),[],1)>1.e-8));
indLRE = (find(max(abs(gp'),[],1)>1.e-8));
@ -154,7 +156,7 @@ if info(1)==0,
options_.analytic_derivation = analytic_derivation;
AHess=-AHess;
if min(eig(AHess))<0,
error('Analytic Hessian is not positive semi-definite!')
error('identification_analysis: Analytic Hessian is not positive semi-definite!')
end
% chol(AHess);
ide_hess.AHess= AHess;
@ -183,7 +185,7 @@ if info(1)==0,
catch,
% replic = max([replic, nparam*(nparam+1)/2*10]);
replic = max([replic, length(indJJ)*3]);
cmm = simulated_moment_uncertainty(indJJ, periods, replic);
cmm = simulated_moment_uncertainty(indJJ, periods, replic,options_,M_,oo_);
% [V,D,W]=eig(cmm);
sd=sqrt(diag(cmm));
cc=cmm./(sd*sd');

83
matlab/simulated_moment_uncertainty.m
View File

@ -1,7 +1,17 @@
function [cmm, mm] = simulated_moment_uncertainty(indx, periods, replic)
% function [cmm, mm] = simulated_moment_uncertainty(indx, periods, replic)
%
function [cmm, mm] = simulated_moment_uncertainty(indx, periods, replic,options_,M_,oo_)
% function [cmm, mm] = simulated_moment_uncertainty(indx, periods, replic,options_,M_,oo_)
% Compute the uncertainty around simulated moments
% Inputs
% - indx [n_moments by 1] index vector of moments
% - periods [scalar] number of simulation periods
% - replic [scalar] number of simulation replications
% - options_ Dynare options structure
% - M_ Dynare Model structure
% - oo_ Dynare results structure
% Outputs:
% - cmm: [n_moments by n_moments] covariance matrix of simulated moments
% - mm: [n_moments by replic] matrix of moments
%
% Copyright (C) 2009-2016 Dynare Team
%
% This file is part of Dynare.
@ -19,20 +29,64 @@ function [cmm, mm] = simulated_moment_uncertainty(indx, periods, replic)
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
global options_ oo_
mm=zeros(length(indx),replic);
disp('Evaluating simulated moment uncertainty ... please wait')
disp(['Doing ',int2str(replic),' replicas of length ',int2str(periods),' periods.'])
noprint0 = options_.noprint;
h = dyn_waitbar(0,'Simulated moment uncertainty ...');
%Do check whether simulation is possible
if options_.periods == 0
error('simulated_moment_uncertainty: Periods must be bigger than 0')
end
if options_.periods <= options_.drop
error('simulated_moment_uncertainty: The horizon of simulation is shorter than the number of observations to be dropped. Either increase options_.periods or decrease options_.drop.')
end
%locally set options
options_.TeX=0;
options_.noprint = 1;
options_.order = 1;
options_.periods = periods;
if isempty(options_.qz_criterium)
options_.qz_criterium = 1+1e-6;
end
if M_.exo_nbr > 0
oo_.exo_simul= ones(max(options_.periods,1) + M_.maximum_lag + M_.maximum_lead,1) * oo_.exo_steady_state';
end
oo_.dr=set_state_space(oo_.dr,M_,options_);
if options_.logged_steady_state %if steady state was previously logged, undo this
oo_.dr.ys=exp(oo_.dr.ys);
oo_.steady_state=exp(oo_.steady_state);
options_.logged_steady_state=0;
logged_steady_state_indicator=1;
evalin('base','options_.logged_steady_state=0;')
end
[dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
oo_.dr=dr;
if info(1)
fprintf('\nsimulated_moment_uncertainty: model could not be solved')
print_info(info,0,options_);
end
%set starting point of simulations
if isempty(M_.endo_histval)
y0 = oo_.dr.ys;
else
y0 = M_.endo_histval;
end
for j=1:replic;
options_.irf = 0;
options_.noprint = 1;
options_.order = 1;
options_.periods = periods;
info = stoch_simul(char(options_.varobs));
[ys, oo_] = simult(y0,oo_.dr,M_,options_,oo_);%do simulation
oo_=disp_moments(ys,char(options_.varobs),M_,options_,oo_); %get moments
dum=[oo_.mean; dyn_vech(oo_.var)];
sd = sqrt(diag(oo_.var));
for i=1:options_.ar;
@ -42,6 +96,9 @@ for j=1:replic;
dyn_waitbar(j/replic,h,['Simulated moment uncertainty. Replic ',int2str(j),'/',int2str(replic)])
end;
dyn_waitbar_close(h);
options_.noprint = noprint0;
if logged_steady_state_indicator
evalin('base','options_.logged_steady_state=1;') %reset base workspace option to conform to base oo_
end
cmm = cov(mm');
disp('Simulated moment uncertainty ... done!')