diff --git a/matlab/dynare_sensitivity.m b/matlab/dynare_sensitivity.m
index 1d9eeef7f..d85ca1f22 100644
--- a/matlab/dynare_sensitivity.m
+++ b/matlab/dynare_sensitivity.m
@@ -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
diff --git a/matlab/endogenous_prior_restrictions.m b/matlab/endogenous_prior_restrictions.m
index e28a726ae..d5b62104b 100644
--- a/matlab/endogenous_prior_restrictions.m
+++ b/matlab/endogenous_prior_restrictions.m
@@ -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 .
-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)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)endo_prior_restrictions.irf{j,4}(1)) && (RR(iendo,iexo)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).
+
+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(probanbr_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(proba0 && 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)ksstat);
- indstab=find(probaksstat);
- indindet=find(probaksstat);
+ indstab=find(probaksstat);
- indunst=find(probaksstat);
+ indindet=find(probaksstat);
+ indunst=find(probaksstat);
+ indwrong=find(proba10,
+ 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(proba10,
+ 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
diff --git a/tests/gsa/ls2003a.mod b/tests/gsa/ls2003a.mod
new file mode 100644
index 000000000..fcbd9b75c
--- /dev/null
+++ b/tests/gsa/ls2003a.mod
@@ -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;
+*/