Provisions for moment calibration.

Sensitivity analysis for irf and moment calibration, with new function /gsa/map_calibration.m
Added new tests/gsa/ls2003a.mod for testsuite of irf/moment calibration

To be done:
- extend SA of calibration to MC and posterior samples (currently uses prior sample only)
- documentaion
- partial progress to close #267
time-shift
Marco Ratto 2014-09-04 16:39:59 +02:00
parent 07d321c159
commit 62b4cfe631
6 changed files with 500 additions and 133 deletions

View File

@ -36,6 +36,10 @@ x0=[];
options_gsa = set_default_option(options_gsa,'datafile',[]);
options_gsa = set_default_option(options_gsa,'rmse',0);
options_gsa = set_default_option(options_gsa,'useautocorr',0);
options_gsa = set_default_option(options_gsa,'moment_calibration',options_.endogenous_prior_restrictions.moment);
options_gsa = set_default_option(options_gsa,'irf_calibration',options_.endogenous_prior_restrictions.irf);
if isfield(options_gsa,'nograph'),
options_.nograph=options_gsa.nograph;
end
@ -100,7 +104,6 @@ if options_gsa.identification,
options_gsa = set_default_option(options_gsa,'trans_ident',0);
options_gsa = set_default_option(options_gsa,'load_ident_files',0);
options_gsa = set_default_option(options_gsa,'ar',1);
options_gsa = set_default_option(options_gsa,'useautocorr',0);
options_.ar = options_gsa.ar;
if options_gsa.morris==0,
disp('The option morris = 0 is no longer supported (Type I errors)')
@ -239,6 +242,10 @@ end
% redform_map(namendo, namlagendo, namexo, icomp, pprior, ilog, threshold)
options_.opt_gsa = options_gsa;
if ~isempty(options_gsa.moment_calibration) || ~isempty(options_gsa.irf_calibration),
map_calibration(OutputDirectoryName, M_, options_, oo_, estim_params_);
end
if options_gsa.identification,
map_ident_(OutputDirectoryName,options_gsa);
end

View File

@ -1,8 +1,8 @@
function info = endogenous_prior_restrictions(T,R,Model,DynareOptions,DynareResults);
% Check for prior (sign) restrictions on irf's
function [info, info_irf, info_moment, data_irf, data_moment] = endogenous_prior_restrictions(T,R,Model,DynareOptions,DynareResults);
% Check for prior (sign) restrictions on irf's and theoretical moments
%
% INPUTS
% T [double] n*n state space matrix
% T [double] n*n state space matrix
% R [double] n*k matrix of shocks
% Model [structure]
% DynareOptions [structure]
@ -11,6 +11,9 @@ function info = endogenous_prior_restrictions(T,R,Model,DynareOptions,DynareResu
% OUTPUTS
% info [double] check if prior restrictions are matched by the
% model and related info
% info_irf [double] array of test checks for all individual irf restrictions
% info_moment [double] array of test checks for all individual moment restrictions
%
% Copyright (C) 2013 Dynare Team
%
@ -29,45 +32,136 @@ function info = endogenous_prior_restrictions(T,R,Model,DynareOptions,DynareResu
% You should have received a copy of the GNU General Public License
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
endo_prior_restrictions= DynareOptions.endogenous_prior_restrictions;
info=0;
if isempty(endo_prior_restrictions.irf),
% no restriction to be checked
return
end
if DynareOptions.order>1,
error('The algorithm for prior (sign) restrictions on irf''s is currently restricted to first-order decision rules')
return
end
infos=[0 0];
varlist=Model.endo_names(DynareResults.dr.order_var,:);
varlist=varlist(DynareResults.dr.restrict_var_list,:);
NT=1;
for j=1:size(endo_prior_restrictions.irf,1),
NT=max(NT,endo_prior_restrictions.irf{j,3});
end
for t=1:NT,
RR = T^(t-1)*R;
for j=1:size(endo_prior_restrictions.irf,1),
if endo_prior_restrictions.irf{j,3}~=t,
continue,
end
iendo=strmatch(endo_prior_restrictions.irf{j,1},varlist,'exact');
iexo=strmatch(endo_prior_restrictions.irf{j,2},Model.exo_names,'exact');
if (RR(iendo,iexo)>endo_prior_restrictions.irf{j,4}(1)) && (RR(iendo,iexo)<endo_prior_restrictions.irf{j,4}(2)),
infos(j,:)=[0, 0];
else
if RR(iendo,iexo)<endo_prior_restrictions.irf{j,4}(1),
delt = (RR(iendo,iexo)-endo_prior_restrictions.irf{j,4}(1))^2;
else
delt = (RR(iendo,iexo)-endo_prior_restrictions.irf{j,4}(2))^2;
end
infos(j,:)=[49, delt];
info=[0 0];
info_irf=[];
info_moment=[];
data_irf=[];
data_moment=[];
endo_prior_restrictions.irf= DynareOptions.endogenous_prior_restrictions.irf;
endo_prior_restrictions.moment= DynareOptions.endogenous_prior_restrictions.moment;
if ~isempty(endo_prior_restrictions.irf),
data_irf=cell(size(endo_prior_restrictions.irf,1),1);
if DynareOptions.order>1,
error('The algorithm for prior (sign) restrictions on irf''s is currently restricted to first-order decision rules')
return
end
varlist=Model.endo_names(DynareResults.dr.order_var,:);
if isempty(T),
[T,R,SteadyState,infox,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults);
else % check if T and R are given in the restricted form!!!
if size(T,1)<size(varlist,1),
varlist=varlist(DynareResults.dr.restrict_var_list,:);
end
end
NT=1;
for j=1:size(endo_prior_restrictions.irf,1),
NT=max(NT,max(endo_prior_restrictions.irf{j,3}));
end
info_irf=ones(size(endo_prior_restrictions.irf,1),2);
for t=1:NT,
if ~DynareOptions.relative_irf,
RR = T^(t-1)*R*diag(sqrt(diag(Model.Sigma_e)));
else
RR = T^(t-1)*R*100;
end
for j=1:size(endo_prior_restrictions.irf,1),
if endo_prior_restrictions.irf{j,3}~=t,
continue,
end
iendo=strmatch(endo_prior_restrictions.irf{j,1},varlist,'exact');
iexo=strmatch(endo_prior_restrictions.irf{j,2},Model.exo_names,'exact');
data_irf{j}=[data_irf{j}; [t RR(iendo,iexo)]];
if (RR(iendo,iexo)>endo_prior_restrictions.irf{j,4}(1)) && (RR(iendo,iexo)<endo_prior_restrictions.irf{j,4}(2)),
info_irf(j,:)=info_irf(j,:).*[0, 0];
else
if RR(iendo,iexo)<endo_prior_restrictions.irf{j,4}(1),
delt = (RR(iendo,iexo)-endo_prior_restrictions.irf{j,4}(1))^2;
else
delt = (RR(iendo,iexo)-endo_prior_restrictions.irf{j,4}(2))^2;
end
info_irf(j,:)=info_irf(j,:).*[49, delt];
end
end
end
if any(info_irf),
info=[49,sum(info_irf(:,2))];
end
end
if any(infos),
info=[49,sum(infos(:,2))];
if ~isempty(endo_prior_restrictions.moment),
if DynareOptions.order>1,
error('The algorithm for prior (sign) restrictions on moments is currently restricted to first-order decision rules')
return
end
data_moment=cell(size(endo_prior_restrictions.moment,1),1);
var_list_=endo_prior_restrictions.moment{1,1};
for j=1:size(endo_prior_restrictions.moment,1),
tmp=endo_prior_restrictions.moment{j,1};
if ~ismember(tmp,cellstr(var_list_)),
var_list_ = char(var_list_, tmp);
end
tmp=endo_prior_restrictions.moment{j,2};
if ~ismember(tmp,cellstr(var_list_)),
var_list_ = char(var_list_, tmp);
end
end
NTmax=0;
NTmin=0;
for j=1:size(endo_prior_restrictions.moment,1),
NTmax=max(NTmax,max(endo_prior_restrictions.moment{j,3}));
NTmin=min(NTmin,min(endo_prior_restrictions.moment{j,3}));
end
info_moment=ones(size(endo_prior_restrictions.moment,1),2);
nvar = size(var_list_,1);
ivar=zeros(nvar,1);
for i=1:nvar
i_tmp = strmatch(var_list_(i,:),Model.endo_names,'exact');
if isempty(i_tmp)
error (['One of the variable specified does not exist']) ;
else
ivar(i) = i_tmp;
end
end
DynareOptions.ar = max(abs(NTmin),NTmax);
[gamma_y,stationary_vars] = th_autocovariances(DynareResults.dr, ivar, Model, DynareOptions,1);
for t=NTmin:NTmax,
RR = gamma_y{abs(t)+1};
if t==0,
RR = RR./(sqrt(diag(RR))*sqrt(diag(RR))')-eye(nvar)+diag(diag(gamma_y{t+1})); % becomes correlation
end
for j=1:size(endo_prior_restrictions.moment,1),
if endo_prior_restrictions.moment{j,3}~=t,
continue,
end
iendo1 = strmatch(endo_prior_restrictions.moment{j,1},var_list_,'exact');
iendo2 = strmatch(endo_prior_restrictions.moment{j,2},var_list_,'exact');
if t>0
tmp0 = iendo1;
iendo1=iendo2;
iendo2=tmp0;
end
data_moment{j}=[data_moment{j}; [t RR(iendo1,iendo2)]];
if (RR(iendo1,iendo2)>endo_prior_restrictions.moment{j,4}(1)) && (RR(iendo1,iendo2)<endo_prior_restrictions.moment{j,4}(2)),
info_moment(j,:)=info_moment(j,:).*[0, 0];
else
if RR(iendo1,iendo2)<endo_prior_restrictions.moment{j,4}(1),
delt = (RR(iendo1,iendo2)-endo_prior_restrictions.moment{j,4}(1))^2;
else
delt = (RR(iendo1,iendo2)-endo_prior_restrictions.moment{j,4}(2))^2;
end
info_moment(j,:)=info_moment(j,:).*[49, delt];
end
end
end
if any(info_moment),
info=[49, info(2) + sum(info_moment(:,2))];
end
end
return

View File

@ -609,6 +609,7 @@ options_.risky_steadystate = 0;
% endogenous prior
options_.endogenous_prior = 0;
options_.endogenous_prior_restrictions.irf={};
options_.endogenous_prior_restrictions.moment={};
% OSR Optimal Simple Rules
options_.osr.tolf=1e-7;

View File

@ -0,0 +1,168 @@
function map_calibration(OutputDirectoryName, Model, DynareOptions, DynareResults, EstimatedParameters)
% Copyright (C) 2014 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 <http://www.gnu.org/licenses/>.
fname_ = Model.fname;
pnames = Model.param_names(EstimatedParameters.param_vals(:,1),:);
pvalue_ks = DynareOptions.opt_gsa.pvalue_ks;
indx_irf = [];
indx_moment = [];
skipline()
disp('Sensitivity analysis for calibration criteria')
filetoload=[OutputDirectoryName '/' fname_ '_prior'];
load(filetoload,'lpmat','lpmat0','istable','iunstable','iindeterm','iwrong' ,'infox')
if ~isempty(lpmat0),
lpmatx=lpmat0(istable,:);
else
lpmatx=[];
end
[Nsam, np] = size(lpmat);
nshock = size(lpmat0,2);
npT = np+nshock;
nbr_irf_restrictions = size(DynareOptions.endogenous_prior_restrictions.irf,1);
mat_irf=cell(nbr_irf_restrictions,1);
for ij=1:nbr_irf_restrictions,
mat_irf{ij}=NaN(Nsam,length(DynareOptions.endogenous_prior_restrictions.irf{ij,3}));
end
nbr_moment_restrictions = size(DynareOptions.endogenous_prior_restrictions.moment,1);
mat_moment=cell(nbr_moment_restrictions,1);
for ij=1:nbr_moment_restrictions,
mat_moment{ij}=NaN(Nsam,length(DynareOptions.endogenous_prior_restrictions.moment{ij,3}));
end
irestrictions = [1:Nsam];
for j=1:Nsam,
Model = set_all_parameters([lpmat0(j,:) lpmat(j,:)]',EstimatedParameters,Model);
[Tt,Rr,SteadyState,info] = dynare_resolve(Model,DynareOptions,DynareResults,'restrict');
if info(1)==0,
[info, info_irf, info_moment, data_irf, data_moment]=endogenous_prior_restrictions(Tt,Rr,Model,DynareOptions,DynareResults);
if ~isempty(info_irf)
for ij=1:nbr_irf_restrictions,
mat_irf{ij}(j,:)=data_irf{ij}(:,2)';
end
indx_irf(j,:)=info_irf(:,1);
end
if ~isempty(info_moment)
for ij=1:nbr_moment_restrictions,
mat_moment{ij}(j,:)=data_moment{ij}(:,2)';
end
indx_moment(j,:)=info_moment(:,1);
end
else
irestrictions(j)=0;
end
end
irestrictions=irestrictions(find(irestrictions));
xmat=[lpmat0(irestrictions,:) lpmat(irestrictions,:)];
skipline()
if ~isempty(indx_irf),
indx_irf = indx_irf(irestrictions,:);
h1=dyn_figure(DynareOptions,'name','Evaluation of irf restrictions');
nrow=ceil(sqrt(nbr_irf_restrictions));
ncol=nrow;
if nrow*(nrow-1)>nbr_irf_restrictions,
ncol=nrow-1;
end
for ij=1:nbr_irf_restrictions,
figure(h1),
mat_irf{ij}=mat_irf{ij}(irestrictions,:);
subplot(nrow,ncol,ij),
for ik=1:size(mat_irf{ij},2),
cumplot(mat_irf{ij}(:,ik)),
hold all,
end
% hist(mat_irf{ij}),
a=axis;
x1val=max(DynareOptions.endogenous_prior_restrictions.irf{ij,4}(1),a(1));
x2val=min(DynareOptions.endogenous_prior_restrictions.irf{ij,4}(2),a(2));
hp = patch([x1val x2val x2val x1val],a([3 3 4 4]),'k');
set(hp,'FaceAlpha',[0.5])
hold off,
leg = num2str(DynareOptions.endogenous_prior_restrictions.irf{ij,3}(1));
if size(mat_irf{ij},2)>1,
leg = [leg,':' ,num2str(DynareOptions.endogenous_prior_restrictions.irf{ij,3}(end))];
end
title([DynareOptions.endogenous_prior_restrictions.irf{ij,1},' vs ',DynareOptions.endogenous_prior_restrictions.irf{ij,2}, '(', leg,')'],'interpreter','none'),
indx1 = find(indx_irf(:,ij)==0);
indx2 = find(indx_irf(:,ij)~=0);
atitle=[DynareOptions.endogenous_prior_restrictions.irf{ij,1},' vs ',DynareOptions.endogenous_prior_restrictions.irf{ij,2}, '(', leg,')'];
fprintf(['%4.1f%% of the prior support matches IRF ',atitle,' inside [%4.1f, %4.1f]\n'],length(indx1)/length(irestrictions)*100,DynareOptions.endogenous_prior_restrictions.irf{ij,4})
aname=['irf_calib_',int2str(ij)];
atitle=['IRF Calib: Parameter(s) driving ',DynareOptions.endogenous_prior_restrictions.irf{ij,1},' vs ',DynareOptions.endogenous_prior_restrictions.irf{ij,2}, '(', leg,')'];
[proba, dproba] = stab_map_1(xmat, indx1, indx2, aname, 0);
indplot=find(proba<pvalue_ks);
if ~isempty(indplot)
stab_map_1(xmat, indx1, indx2, aname, 1, indplot, OutputDirectoryName,[],atitle);
end
end
dyn_saveas(h1,[OutputDirectoryName,filesep,fname_,'_irf_restrictions'],DynareOptions);
skipline()
end
if ~isempty(indx_moment)
indx_moment = indx_moment(irestrictions,:);
h2=dyn_figure(DynareOptions,'name','Evaluation of moment restrictions');
nrow=ceil(sqrt(nbr_moment_restrictions));
ncol=nrow;
if nrow*(nrow-1)>nbr_moment_restrictions,
ncol=nrow-1;
end
for ij=1:nbr_moment_restrictions,
figure(h2),
mat_moment{ij}=mat_moment{ij}(irestrictions,:);
subplot(nrow,ncol,ij),
for ik=1:size(mat_moment{ij},2),
cumplot(mat_moment{ij}(:,ik)),
hold all,
end
% hist(mat_moment{ij}),
a=axis;
x1val=max(DynareOptions.endogenous_prior_restrictions.moment{ij,4}(1),a(1));
x2val=min(DynareOptions.endogenous_prior_restrictions.moment{ij,4}(2),a(2));
hp = patch([x1val x2val x2val x1val],a([3 3 4 4]),'k');
set(hp,'FaceAlpha',[0.5])
hold off,
leg = num2str(DynareOptions.endogenous_prior_restrictions.moment{ij,3}(1));
if size(mat_moment{ij},2)>1,
leg = [leg,':' ,num2str(DynareOptions.endogenous_prior_restrictions.moment{ij,3}(end))];
end
title([DynareOptions.endogenous_prior_restrictions.moment{ij,1},' vs ',DynareOptions.endogenous_prior_restrictions.moment{ij,2},'(',leg,')'],'interpreter','none'),
indx1 = find(indx_moment(:,ij)==0);
indx2 = find(indx_moment(:,ij)~=0);
atitle=[DynareOptions.endogenous_prior_restrictions.moment{ij,1},' vs ',DynareOptions.endogenous_prior_restrictions.moment{ij,2}, '(', leg,')'];
fprintf(['%4.1f%% of the prior support matches MOMENT ',atitle,' inside [%4.1f, %4.1f]\n'],length(indx1)/length(irestrictions)*100,DynareOptions.endogenous_prior_restrictions.moment{ij,4})
aname=['moment_calib_',int2str(ij)];
atitle=['MOMENT Calib: Parameter(s) driving ',DynareOptions.endogenous_prior_restrictions.moment{ij,1},' vs ',DynareOptions.endogenous_prior_restrictions.moment{ij,2}, '(', leg,')'];
[proba, dproba] = stab_map_1(xmat, indx1, indx2, aname, 0);
indplot=find(proba<pvalue_ks);
if ~isempty(indplot)
stab_map_1(xmat, indx1, indx2, aname, 1, indplot, OutputDirectoryName,[],atitle);
end
end
dyn_saveas(h2,[OutputDirectoryName,filesep,fname_,'_moment_restrictions'],DynareOptions);
skipline()
end
return

View File

@ -253,7 +253,7 @@ if fload==0,
iindeterm=zeros(1,Nsam);
iwrong=zeros(1,Nsam);
for j=1:Nsam,
M_.params(estim_params_.param_vals(:,1)) = lpmat(j,:)';
M_ = set_all_parameters([lpmat0(j,:) lpmat(j,:)]',estim_params_,M_);
%try stoch_simul([]);
try
[Tt,Rr,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_,'restrict');
@ -496,7 +496,7 @@ delete([OutputDirectoryName,filesep,fname_,'_',auname,'_corr_*.*']);
delete([OutputDirectoryName,filesep,fname_,'_',aunstname,'_corr_*.*']);
delete([OutputDirectoryName,filesep,fname_,'_',aindname,'_corr_*.*']);
if length(iunstable)>0 && length(iunstable)<Nsam,
if length(iunstable)>0 ,
fprintf(['%4.1f%% of the prior support gives unique saddle-path solution.\n'],length(istable)/Nsam*100)
fprintf(['%4.1f%% of the prior support gives explosive dynamics.\n'],(length(iunstable)-length(iwrong)-length(iindeterm) )/Nsam*100)
if ~isempty(iindeterm),
@ -504,137 +504,137 @@ if length(iunstable)>0 && length(iunstable)<Nsam,
end
if ~isempty(iwrong),
skipline()
disp(['For ',num2str(length(iwrong)/Nsam*100,'%1.3f'),'\% of the prior support dynare could not find a solution.'])
disp(['For ',num2str(length(iwrong)/Nsam*100,'%4.1f'),'% of the prior support dynare could not find a solution.'])
skipline()
if any(infox==1),
disp([' For ',num2str(length(find(infox==1))/Nsam*100,'%1.3f'),'\% The model doesn''t determine the current variables uniquely.'])
disp([' For ',num2str(length(find(infox==1))/Nsam*100,'%4.1f'),'% The model doesn''t determine the current variables uniquely.'])
end
if any(infox==2),
disp([' For ',num2str(length(find(infox==2))/Nsam*100,'%1.3f'),'\% MJDGGES returned an error code.'])
disp([' For ',num2str(length(find(infox==2))/Nsam*100,'%4.1f'),'% MJDGGES returned an error code.'])
end
if any(infox==6),
disp([' For ',num2str(length(find(infox==6))/Nsam*100,'%1.3f'),'\% The jacobian evaluated at the deterministic steady state is complex.'])
disp([' For ',num2str(length(find(infox==6))/Nsam*100,'%4.1f'),'% The jacobian evaluated at the deterministic steady state is complex.'])
end
if any(infox==19),
disp([' For ',num2str(length(find(infox==19))/Nsam*100,'%1.3f'),'\% The steadystate routine thrown an exception (inconsistent deep parameters).'])
disp([' For ',num2str(length(find(infox==19))/Nsam*100,'%4.1f'),'% The steadystate routine thrown an exception (inconsistent deep parameters).'])
end
if any(infox==20),
disp([' For ',num2str(length(find(infox==20))/Nsam*100,'%1.3f'),'\% Cannot find the steady state.'])
disp([' For ',num2str(length(find(infox==20))/Nsam*100,'%4.1f'),'% Cannot find the steady state.'])
end
if any(infox==21),
disp([' For ',num2str(length(find(infox==21))/Nsam*100,'%1.3f'),'\% The steady state is complex.'])
disp([' For ',num2str(length(find(infox==21))/Nsam*100,'%4.1f'),'% The steady state is complex.'])
end
if any(infox==22),
disp([' For ',num2str(length(find(infox==22))/Nsam*100,'%1.3f'),'\% The steady has NaNs.'])
disp([' For ',num2str(length(find(infox==22))/Nsam*100,'%4.1f'),'% The steady has NaNs.'])
end
if any(infox==23),
disp([' For ',num2str(length(find(infox==23))/Nsam*100,'%1.3f'),'\% M_.params has been updated in the steadystate routine and has complex valued scalars.'])
disp([' For ',num2str(length(find(infox==23))/Nsam*100,'%4.1f'),'% M_.params has been updated in the steadystate routine and has complex valued scalars.'])
end
if any(infox==24),
disp([' For ',num2str(length(find(infox==24))/Nsam*100,'%1.3f'),'\% M_.params has been updated in the steadystate routine and has some NaNs.'])
disp([' For ',num2str(length(find(infox==24))/Nsam*100,'%4.1f'),'% M_.params has been updated in the steadystate routine and has some NaNs.'])
end
if any(infox==30),
disp([' For ',num2str(length(find(infox==30))/Nsam*100,'%1.3f'),'\% Ergodic variance can''t be computed.'])
disp([' For ',num2str(length(find(infox==30))/Nsam*100,'%4.1f'),'% Ergodic variance can''t be computed.'])
end
if any(infox==49),
disp([' For ',num2str(length(find(infox==49))/Nsam*100,'%1.3f'),'\% The model violates one (many) endogenous prior restriction(s).'])
disp([' For ',num2str(length(find(infox==49))/Nsam*100,'%4.1f'),'% The model violates one (many) endogenous prior restriction(s).'])
end
end
skipline()
% Blanchard Kahn
[proba, dproba] = stab_map_1(lpmat, istable, iunstable, aname,0);
% indstab=find(dproba>ksstat);
indstab=find(proba<pvalue_ks);
disp('Smirnov statistics in driving acceptable behaviour')
for j=1:length(indstab),
disp([M_.param_names(estim_params_.param_vals(indstab(j),1),:),' d-stat = ', num2str(dproba(indstab(j)),'%1.3f'),' p-value = ', num2str(proba(indstab(j)),'%1.3f')])
end
skipline()
if ~isempty(indstab)
stab_map_1(lpmat, istable, iunstable, aname, 1, indstab, OutputDirectoryName,[],atitle);
end
ixun=iunstable(find(~ismember(iunstable,[iindeterm,iwrong])));
if ~isempty(iindeterm),
[proba, dproba] = stab_map_1(lpmat, [1:Nsam], iindeterm, aindetname ,0);
% indindet=find(dproba>ksstat);
indindet=find(proba<pvalue_ks);
disp('Smirnov statistics in driving indeterminacy')
for j=1:length(indindet),
disp([M_.param_names(estim_params_.param_vals(indindet(j),1),:),' d-stat = ', num2str(dproba(indindet(j)),'%1.3f'),' p-value = ', num2str(proba(indindet(j)),'%1.3f')])
if length(iunstable)<Nsam
% Blanchard Kahn
[proba, dproba] = stab_map_1(lpmat, istable, iunstable, aname,0);
% indstab=find(dproba>ksstat);
indstab=find(proba<pvalue_ks);
disp('Smirnov statistics in driving acceptable behaviour')
for j=1:length(indstab),
disp([M_.param_names(estim_params_.param_vals(indstab(j),1),:),' d-stat = ', num2str(dproba(indstab(j)),'%1.3f'),' p-value = ', num2str(proba(indstab(j)),'%1.3f')])
end
skipline()
if ~isempty(indindet)
stab_map_1(lpmat, [1:Nsam], iindeterm, aindetname, 1, indindet, OutputDirectoryName,[],aindettitle);
if ~isempty(indstab)
stab_map_1(lpmat, istable, iunstable, aname, 1, indstab, OutputDirectoryName,[],atitle);
end
end
if ~isempty(ixun),
[proba, dproba] = stab_map_1(lpmat, [1:Nsam], ixun, aunstablename,0);
% indunst=find(dproba>ksstat);
indunst=find(proba<pvalue_ks);
disp('Smirnov statistics in driving instability')
for j=1:length(indunst),
disp([M_.param_names(estim_params_.param_vals(indunst(j),1),:),' d-stat = ', num2str(dproba(indunst(j)),'%1.3f'),' p-value = ', num2str(proba(indunst(j)),'%1.3f')])
ixun=iunstable(find(~ismember(iunstable,[iindeterm,iwrong])));
if ~isempty(iindeterm),
[proba, dproba] = stab_map_1(lpmat, [1:Nsam], iindeterm, aindetname ,0);
% indindet=find(dproba>ksstat);
indindet=find(proba<pvalue_ks);
disp('Smirnov statistics in driving indeterminacy')
for j=1:length(indindet),
disp([M_.param_names(estim_params_.param_vals(indindet(j),1),:),' d-stat = ', num2str(dproba(indindet(j)),'%1.3f'),' p-value = ', num2str(proba(indindet(j)),'%1.3f')])
end
skipline()
if ~isempty(indindet)
stab_map_1(lpmat, [1:Nsam], iindeterm, aindetname, 1, indindet, OutputDirectoryName,[],aindettitle);
end
end
if ~isempty(ixun),
[proba, dproba] = stab_map_1(lpmat, [1:Nsam], ixun, aunstablename,0);
% indunst=find(dproba>ksstat);
indunst=find(proba<pvalue_ks);
disp('Smirnov statistics in driving instability')
for j=1:length(indunst),
disp([M_.param_names(estim_params_.param_vals(indunst(j),1),:),' d-stat = ', num2str(dproba(indunst(j)),'%1.3f'),' p-value = ', num2str(proba(indunst(j)),'%1.3f')])
end
skipline()
if ~isempty(indunst)
stab_map_1(lpmat, [1:Nsam], ixun, aunstablename, 1, indunst, OutputDirectoryName,[],aunstabletitle);
end
end
if ~isempty(iwrong),
[proba, dproba] = stab_map_1(lpmat, [1:Nsam], iwrong, awronguniname,0);
% indwrong=find(dproba>ksstat);
indwrong=find(proba<pvalue_ks);
disp('Smirnov statistics in driving no solution')
for j=1:length(indwrong),
disp([M_.param_names(estim_params_.param_vals(indwrong(j),1),:),' d-stat = ', num2str(dproba(indwrong(j)),'%1.3f'),' p-value = ', num2str(proba(indwrong(j)),'%1.3f')])
end
skipline()
if ~isempty(indwrong)
stab_map_1(lpmat, [1:Nsam], iwrong, awronguniname, 1, indwrong, OutputDirectoryName,[],awrongunititle);
end
end
skipline()
if ~isempty(indunst)
stab_map_1(lpmat, [1:Nsam], ixun, aunstablename, 1, indunst, OutputDirectoryName,[],aunstabletitle);
disp('Starting bivariate analysis:')
c0=corrcoef(lpmat(istable,:));
c00=tril(c0,-1);
if length(istable)>10,
stab_map_2(lpmat(istable,:),alpha2, pvalue_corr, asname, OutputDirectoryName,xparam1,astitle);
end
end
if ~isempty(iwrong),
[proba, dproba] = stab_map_1(lpmat, [1:Nsam], iwrong, awronguniname,0);
% indwrong=find(dproba>ksstat);
indwrong=find(proba<pvalue_ks);
disp('Smirnov statistics in driving no solution')
for j=1:length(indwrong),
disp([M_.param_names(estim_params_.param_vals(indwrong(j),1),:),' d-stat = ', num2str(dproba(indwrong(j)),'%1.3f'),' p-value = ', num2str(proba(indwrong(j)),'%1.3f')])
if length(iunstable)>10,
stab_map_2(lpmat(iunstable,:),alpha2, pvalue_corr, auname, OutputDirectoryName,xparam1,autitle);
end
skipline()
if ~isempty(indwrong)
stab_map_1(lpmat, [1:Nsam], iwrong, awronguniname, 1, indwrong, OutputDirectoryName,[],awrongunititle);
if length(iindeterm)>10,
stab_map_2(lpmat(iindeterm,:),alpha2, pvalue_corr, aindname, OutputDirectoryName,xparam1,aindtitle);
end
end
skipline()
disp('Starting bivariate analysis:')
c0=corrcoef(lpmat(istable,:));
c00=tril(c0,-1);
if length(istable)>10,
stab_map_2(lpmat(istable,:),alpha2, pvalue_corr, asname, OutputDirectoryName,xparam1,astitle);
end
if length(iunstable)>10,
stab_map_2(lpmat(iunstable,:),alpha2, pvalue_corr, auname, OutputDirectoryName,xparam1,autitle);
end
if length(iindeterm)>10,
stab_map_2(lpmat(iindeterm,:),alpha2, pvalue_corr, aindname, OutputDirectoryName,xparam1,aindtitle);
end
if length(ixun)>10,
stab_map_2(lpmat(ixun,:),alpha2, pvalue_corr, aunstname, OutputDirectoryName,xparam1,aunsttitle);
end
if length(iwrong)>10,
stab_map_2(lpmat(iwrong,:),alpha2, pvalue_corr, awrongname, OutputDirectoryName,xparam1,awrongtitle);
end
x0=0.5.*(bayestopt_.ub(1:nshock)-bayestopt_.lb(1:nshock))+bayestopt_.lb(1:nshock);
x0 = [x0; lpmat(istable(1),:)'];
if istable(end)~=Nsam
M_.params(estim_params_.param_vals(:,1)) = lpmat(istable(1),:)';
[oo_.dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
% stoch_simul([]);
end
else
if length(iunstable)==0,
disp('All parameter values in the specified ranges give unique saddle-path solution!')
if length(ixun)>10,
stab_map_2(lpmat(ixun,:),alpha2, pvalue_corr, aunstname, OutputDirectoryName,xparam1,aunsttitle);
end
if length(iwrong)>10,
stab_map_2(lpmat(iwrong,:),alpha2, pvalue_corr, awrongname, OutputDirectoryName,xparam1,awrongtitle);
end
x0=0.5.*(bayestopt_.ub(1:nshock)-bayestopt_.lb(1:nshock))+bayestopt_.lb(1:nshock);
x0 = [x0; lpmat(istable(1),:)'];
if istable(end)~=Nsam
M_.params(estim_params_.param_vals(:,1)) = lpmat(istable(1),:)';
[oo_.dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
% stoch_simul([]);
end
else
disp('All parameter values in the specified ranges are not acceptable!')
x0=[];
end
else
disp('All parameter values in the specified ranges give unique saddle-path solution!')
x0=0.5.*(bayestopt_.ub(1:nshock)-bayestopt_.lb(1:nshock))+bayestopt_.lb(1:nshock);
x0 = [x0; lpmat(istable(1),:)'];
end

97
tests/gsa/ls2003a.mod Normal file
View File

@ -0,0 +1,97 @@
var y y_s R pie dq pie_s de A y_obs pie_obs R_obs;
varexo e_R e_q e_ys e_pies e_A;
parameters psi1 psi2 psi3 rho_R tau alpha rr k rho_q rho_A rho_ys rho_pies;
psi1 = 1.54;
psi2 = 0.25;
psi3 = 0.25;
rho_R = 0.5;
alpha = 0.3;
rr = 2.51;
k = 0.5;
tau = 0.5;
rho_q = 0.4;
rho_A = 0.2;
rho_ys = 0.9;
rho_pies = 0.7;
model(linear);
y = y(+1) - (tau +alpha*(2-alpha)*(1-tau))*(R-pie(+1))-alpha*(tau +alpha*(2-alpha)*(1-tau))*dq(+1) + alpha*(2-alpha)*((1-tau)/tau)*(y_s-y_s(+1))-A(+1);
pie = exp(-rr/400)*pie(+1)+alpha*exp(-rr/400)*dq(+1)-alpha*dq+(k/(tau+alpha*(2-alpha)*(1-tau)))*y+k*alpha*(2-alpha)*(1-tau)/(tau*(tau+alpha*(2-alpha)*(1-tau)))*y_s;
pie = de+(1-alpha)*dq+pie_s;
R = rho_R*R(-1)+(1-rho_R)*(psi1*pie+psi2*(y+alpha*(2-alpha)*((1-tau)/tau)*y_s)+psi3*de)+e_R;
dq = rho_q*dq(-1)+e_q;
y_s = rho_ys*y_s(-1)+e_ys;
pie_s = rho_pies*pie_s(-1)+e_pies;
A = rho_A*A(-1)+e_A;
y_obs = y-y(-1)+A;
pie_obs = 4*pie;
R_obs = 4*R;
end;
shocks;
var e_R = 1.25^2;
var e_q = 2.5^2;
var e_A = 1.89;
var e_ys = 1.89;
var e_pies = 1.89;
end;
varobs y_obs R_obs pie_obs dq de;
estimated_params;
psi1 , gamma_pdf,1.5,0.5;
psi2 , gamma_pdf,0.25,0.125;
psi3 , gamma_pdf,0.25,0.125;
rho_R ,beta_pdf,0.5,0.2;
alpha ,beta_pdf,0.3,0.1;
rr ,gamma_pdf,2.5,1;
k , gamma_pdf,0.5,0.25;
tau ,gamma_pdf,0.5,0.2;
rho_q ,beta_pdf,0.4,0.2;
rho_A ,beta_pdf,0.5,0.2;
rho_ys ,beta_pdf,0.8,0.1;
rho_pies,beta_pdf,0.7,0.15;
stderr e_R,inv_gamma_pdf,1.2533,0.6551;
stderr e_q,inv_gamma_pdf,2.5066,1.3103;
stderr e_A,inv_gamma_pdf,1.2533,0.6551;
stderr e_ys,inv_gamma_pdf,1.2533,0.6551;
stderr e_pies,inv_gamma_pdf,1.88,0.9827;
end;
// endogenous prior restrictions
options_.relative_irf=1;
irf_calibration;
y(1:4), e_ys, [ -50 50]; //[first year response]
//y(1:4), e_ys, [-inf -50]; //[first year response]
@#for ilag in 21:40
R_obs(@{ilag}), e_ys, [0 6]; //[response after 4th year to 10th year]
@#endfor
end;
/*
irf_calibration;
y(1:4), e_ys, [-inf -0.4]; //[first year response]
@#for ilag in 21:40
R_obs(@{ilag}), e_ys, [0 0.25]; //[response after 4th year to 10th year]
@#endfor
end;
*/
moment_calibration;
//y_obs,y_obs, [0.8 1.1]; //[unconditional variance]
y_obs,y_obs(-(1:4)), +; //[first year acf]
@#for ilag in 0:1
y_obs,R_obs(-@{ilag}), -; //[ccf]
@#endfor
end;
dynare_sensitivity(prior_range=0);
/*
estimation(datafile='data_ca1.m',first_obs=8,nobs=79,mh_nblocks=2, mode_file = ls2003a_mode,
prefilter=1,mh_jscale=0.5,mh_replic=5000, mode_compute=0, mh_drop=0.6, bayesian_irf);
//stoch_simul(irf=40, order=1, relative_irf) y_obs R_obs pie_obs dq de;
stoch_simul(irf=40, order=1) y_obs R_obs pie_obs dq de;
*/