diff --git a/matlab/dynare_identification.m b/matlab/dynare_identification.m index 7ce0a2b8e..b98b549a1 100644 --- a/matlab/dynare_identification.m +++ b/matlab/dynare_identification.m @@ -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_; diff --git a/matlab/identification_analysis.m b/matlab/identification_analysis.m index e77b9cfa1..855c16049 100644 --- a/matlab/identification_analysis.m +++ b/matlab/identification_analysis.m @@ -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)1.e-8)); end if length(indJJ)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'); diff --git a/matlab/simulated_moment_uncertainty.m b/matlab/simulated_moment_uncertainty.m index 09441561d..1fc16307c 100644 --- a/matlab/simulated_moment_uncertainty.m +++ b/matlab/simulated_moment_uncertainty.m @@ -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 . -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!')