diff --git a/doc/gsa/gsa.tex b/doc/gsa/gsa.tex
index 3e9bcdbdf30a957000b82836b69736b6662b492f..04eb8f042d09d1d3e7162764594d58e8a862ed2b 100644
--- a/doc/gsa/gsa.tex
+++ b/doc/gsa/gsa.tex
@@ -167,9 +167,9 @@ out-of-sample validation. \vspace{0.5cm}
& & of reduced form coefficients\\
& & [\verb"max" \verb"max"] = analyse filtered \\
& & entries within the range [\verb"max" \verb"max"] \\
- \verb"ksstat_redform"& 0.1& critical value for Smirnov statistics $d$ \\
+ \verb"ksstat_redform"& 0.001& p-value for Smirnov statistics $d$ \\
& & when reduced form entries are filtered\\
- \verb"alpha2_redform"& 0.3& critical value for correlation $\rho$ \\
+ \verb"alpha2_redform"& 0.001& p-value for correlation $\rho$ \\
& & when reduced form entries are filtered\\
\verb"namendo"& () & list of endogenous variables \\
& : & jolly character to indicate ALL endogenous \\
@@ -234,12 +234,11 @@ options of \verb"dynare_estimation". These are:
& & filter the best 10\% for each observed series\\
\verb"istart_rmse"& 1& start computing RMSE's from \verb"istart_rmse":\\
& & use 2 to avoid big initial error \\
- \verb"alpha_rmse"& 0.002& critical value for Smirnov statistics $d$:\\
- & & plot parameters with $d>$\verb"alpha_rmse"\\
- \verb"alpha2_rmse"& 1& critical value for correlation $\rho$\\
+ \verb"alpha_rmse"& 0.001& p-value for Smirnov statistics $d$:\\
+ & & plot parameters with $pvalue<$\verb"alpha_rmse"\\
+ \verb"alpha2_rmse"& 0.001& p-value for correlation $\rho$\\
& & plot couples of parameters with
- $|\rho|>$\verb"alpha2_rmse"\\
- \verb"glue"& 0& prepare for GLUE graphical interface\\\hline
+ $pvalue<$\verb"alpha2_rmse"\\
\end{tabular}
\subsection{Screening analysis}
diff --git a/matlab/dynare_sensitivity.m b/matlab/dynare_sensitivity.m
index 2b9542dd3dc1951bf77d1de2c12663acaa3f95f3..7aef12582fe3e1b88f37c24f98c9c652cb8c5223 100644
--- a/matlab/dynare_sensitivity.m
+++ b/matlab/dynare_sensitivity.m
@@ -57,6 +57,9 @@ if ~isempty(options_gsa.datafile) || isempty(bayestopt_),
if isfield(options_gsa,'loglinear'),
options_.loglinear=options_gsa.loglinear;
end
+ if isfield(options_gsa,'lik_init'),
+ options_.lik_init=options_gsa.lik_init;
+ end
options_.mode_compute = 0;
options_.filtered_vars = 1;
options_.plot_priors = 0;
@@ -121,8 +124,24 @@ options_gsa = set_default_option(options_gsa,'ksstat',0.1);
options_gsa = set_default_option(options_gsa,'pvalue_ks',0.001);
options_gsa = set_default_option(options_gsa,'pvalue_corr',0.001);
%options_gsa = set_default_option(options_gsa,'load_mh',0);
+% REDFORM mapping
+options_gsa = set_default_option(options_gsa,'logtrans_redform',0);
+options_gsa = set_default_option(options_gsa,'threshold_redform',[]);
+options_gsa = set_default_option(options_gsa,'ksstat_redform',0.001);
+options_gsa = set_default_option(options_gsa,'alpha2_redform',0.001);
+options_gsa = set_default_option(options_gsa,'namendo',[]);
+options_gsa = set_default_option(options_gsa,'namlagendo',[]);
+options_gsa = set_default_option(options_gsa,'namexo',[]);
+% RMSE mapping
+options_gsa = set_default_option(options_gsa,'rmse',0);
+options_gsa = set_default_option(options_gsa,'lik_only',0);
+options_gsa = set_default_option(options_gsa,'var_rmse',options_.varobs);
+options_gsa = set_default_option(options_gsa,'pfilt_rmse',0.1);
+options_gsa = set_default_option(options_gsa,'istart_rmse',options_.presample+1);
+options_gsa = set_default_option(options_gsa,'alpha_rmse',0.001);
+options_gsa = set_default_option(options_gsa,'alpha2_rmse',0.001);
-if options_gsa.redform,
+if options_gsa.redform && options_gsa.neighborhood_width==0 && isempty(options_gsa.threshold_redform),
options_gsa.pprior=1;
options_gsa.ppost=0;
end
@@ -145,6 +164,8 @@ if options_gsa.morris==1,
options_gsa.load_rmse=0;
options_gsa.alpha2_stab=1;
options_gsa.ksstat=1;
+ options_gsa.pvalue_ks=0;
+ options_gsa.pvalue_corr=0;
% if options_gsa.morris==3,
% options_gsa = set_default_option(options_gsa,'Nsam',256);
% OutputDirectoryName = CheckPath('gsa/identif',M_.dname);
@@ -155,7 +176,7 @@ else
OutputDirectoryName = CheckPath('gsa',M_.dname);
end
-options_.opt_gsa = options_gsa;
+% options_.opt_gsa = options_gsa;
if (options_gsa.load_stab || options_gsa.load_rmse || options_gsa.load_redform) ...
&& options_gsa.pprior,
@@ -187,22 +208,15 @@ if (options_gsa.load_stab || options_gsa.load_rmse || options_gsa.load_redform)
end
if options_gsa.stab && ~options_gsa.ppost,
- x0 = stab_map_(OutputDirectoryName);
+ x0 = stab_map_(OutputDirectoryName,options_gsa);
end
% reduced form
% redform_map(namendo, namlagendo, namexo, icomp, pprior, ilog, threshold)
-options_gsa = set_default_option(options_gsa,'logtrans_redform',0);
-options_gsa = set_default_option(options_gsa,'threshold_redform',[]);
-options_gsa = set_default_option(options_gsa,'ksstat_redform',0.1);
-options_gsa = set_default_option(options_gsa,'alpha2_redform',0.3);
-options_gsa = set_default_option(options_gsa,'namendo',[]);
-options_gsa = set_default_option(options_gsa,'namlagendo',[]);
-options_gsa = set_default_option(options_gsa,'namexo',[]);
options_.opt_gsa = options_gsa;
if options_gsa.identification,
- map_ident_(OutputDirectoryName);
+ map_ident_(OutputDirectoryName,options_gsa);
end
if options_gsa.redform && ~isempty(options_gsa.namendo) ...
@@ -216,30 +230,24 @@ if options_gsa.redform && ~isempty(options_gsa.namendo) ...
if strmatch(':',options_gsa.namlagendo,'exact'),
options_gsa.namlagendo=M_.endo_names;
end
- options_.opt_gsa = options_gsa;
+% options_.opt_gsa = options_gsa;
if options_gsa.morris==1,
- redform_screen(OutputDirectoryName);
+ redform_screen(OutputDirectoryName,options_gsa);
else
% check existence of the SS_ANOVA toolbox
- if ~(exist('gsa_sdp','file')==6 || exist('gsa_sdp','file')==2),
+ if isempty(options_gsa.threshold_redform) && ...
+ ~(exist('gsa_sdp','file')==6 || exist('gsa_sdp','file')==2),
disp('Download Mapping routines at:')
disp('http://eemc.jrc.ec.europa.eu/softwareDYNARE-Dowload.htm')
disp(' ' )
error('Mapping routines missing!')
end
- redform_map(OutputDirectoryName);
+ redform_map(OutputDirectoryName,options_gsa);
end
end
% RMSE mapping
% function [rmse_MC, ixx] = filt_mc_(vvarvecm, loadSA, pfilt, alpha, alpha2)
-options_gsa = set_default_option(options_gsa,'rmse',0);
-options_gsa = set_default_option(options_gsa,'lik_only',0);
-options_gsa = set_default_option(options_gsa,'var_rmse',options_.varobs);
-options_gsa = set_default_option(options_gsa,'pfilt_rmse',0.1);
-options_gsa = set_default_option(options_gsa,'istart_rmse',options_.presample+1);
-options_gsa = set_default_option(options_gsa,'alpha_rmse',0.002);
-options_gsa = set_default_option(options_gsa,'alpha2_rmse',1);
options_.opt_gsa = options_gsa;
if options_gsa.rmse,
if ~options_gsa.ppost
@@ -260,6 +268,12 @@ if options_gsa.rmse,
end
end
if isempty(a),
+ if options_gsa.lik_only,
+ options_.smoother=0;
+ options_.filter_step_ahead=[];
+ options_.forecast=0;
+ options_.filtered_vars=0;
+ end
% dynare_MC([],OutputDirectoryName,data,rawdata,data_info);
prior_posterior_statistics('gsa',dataset_);
if options_.bayesian_irf
@@ -284,8 +298,9 @@ if options_gsa.rmse,
end
clear a;
% filt_mc_(OutputDirectoryName,data_info);
- filt_mc_(OutputDirectoryName);
+ filt_mc_(OutputDirectoryName,options_gsa);
end
+options_.opt_gsa = options_gsa;
if options_gsa.glue,
diff --git a/matlab/gsa/filt_mc_.m b/matlab/gsa/filt_mc_.m
index 82ecebb3be08632bb3698cfaf3fa5e7d464cc8ac..6cc14ce896d79c007ec753fe5fb93daaf3275905 100644
--- a/matlab/gsa/filt_mc_.m
+++ b/matlab/gsa/filt_mc_.m
@@ -1,4 +1,4 @@
-function [rmse_MC, ixx] = filt_mc_(OutDir,data_info)
+function [rmse_MC, ixx] = filt_mc_(OutDir,options_gsa_)
% function [rmse_MC, ixx] = filt_mc_(OutDir)
% inputs (from opt_gsa structure)
% vvarvecm = options_gsa_.var_rmse;
@@ -35,13 +35,14 @@ function [rmse_MC, ixx] = filt_mc_(OutDir,data_info)
global bayestopt_ estim_params_ M_ options_ oo_
-options_gsa_=options_.opt_gsa;
+% options_gsa_=options_.opt_gsa;
vvarvecm = options_gsa_.var_rmse;
loadSA = options_gsa_.load_rmse;
pfilt = options_gsa_.pfilt_rmse;
alpha = options_gsa_.alpha_rmse;
-alpha2 = options_gsa_.alpha2_rmse;
-pvalue = options_gsa_.pvalue_corr;
+% alpha2 = options_gsa_.alpha2_rmse;
+alpha2 = 0;
+pvalue = options_gsa_.alpha2_rmse;
istart = options_gsa_.istart_rmse;
alphaPC = 0.5;
@@ -210,9 +211,9 @@ if ~loadSA,
save([OutDir,filesep,fnamtmp,'.mat'], 'x', 'logpo2', 'likelihood', 'rmse_MC', 'rmse_mode','rmse_pmean')
else
if options_.opt_gsa.lik_only
- save([OutDir,filesep,fnamtmp, '.mat'], 'likelihood', '-append')
+ save([OutDir,filesep,fnamtmp, '.mat'], 'x', 'logpo2','likelihood', '-append')
else
- save([OutDir,filesep,fnamtmp, '.mat'], 'likelihood', 'rmse_MC','-append')
+ save([OutDir,filesep,fnamtmp, '.mat'], 'x', 'logpo2','likelihood', 'rmse_MC','-append')
if exist('xparam1_mean','var')
save([OutDir,filesep,fnamtmp, '.mat'], 'rmse_pmean','-append')
end
@@ -222,7 +223,7 @@ if ~loadSA,
end
end
else
- if options_.opt_gsa.lik_only & options_.opt_gsa.ppost==0
+ if options_.opt_gsa.lik_only && options_.opt_gsa.ppost==0
load([OutDir,filesep,fnamtmp, '.mat'],'x','logpo2','likelihood');
else
load([OutDir,filesep,fnamtmp, '.mat'],'x','logpo2','likelihood','rmse_MC','rmse_mode','rmse_pmean');
@@ -238,21 +239,21 @@ if ~options_.opt_gsa.ppost
[dum, ipost]=sort(-logpo2);
[dum, ilik]=sort(-likelihood);
end
-if ~options_.opt_gsa.ppost & options_.opt_gsa.lik_only
+if ~options_.opt_gsa.ppost && options_.opt_gsa.lik_only
if options_.opt_gsa.pprior
anam='rmse_prior_post';
else
anam='rmse_mc_post';
end
stab_map_1(x, ipost(1:nfilt), ipost(nfilt+1:end), anam, 1,[],OutDir);
- stab_map_2(x(ipost(1:nfilt),:),alpha2,anam, OutDir);
+ stab_map_2(x(ipost(1:nfilt),:),alpha2,pvalue,anam, OutDir);
if options_.opt_gsa.pprior
anam='rmse_prior_lik';
else
anam='rmse_mc_lik';
end
stab_map_1(x, ilik(1:nfilt), ilik(nfilt+1:end), anam, 1,[],OutDir);
- stab_map_2(x(ilik(1:nfilt),:),alpha2,anam, OutDir);
+ stab_map_2(x(ilik(1:nfilt),:),alpha2,pvalue,anam, OutDir);
else
for i=1:size(vvarvecm,1),
[dum, ixx(:,i)]=sort(rmse_MC(:,i));
@@ -260,7 +261,7 @@ else
%nfilt0(i)=length(find(rmse_MC(:,i)<rmse_pmean(i)));
rmse_txt=rmse_pmean;
else
- if options_.opt_gsa.pprior | ~exist('rmse_pmean'),
+ if options_.opt_gsa.pprior || ~exist('rmse_pmean'),
if exist('rmse_mode'),
rmse_txt=rmse_mode;
else
@@ -298,7 +299,7 @@ else
h=cumplot(lnprior(ixx(nfilt0(i)+1:end,i)));
set(h,'color','green')
title(vvarvecm(i,:),'interpreter','none')
- if mod(i,9)==0 | i==size(vvarvecm,1)
+ if mod(i,9)==0 || i==size(vvarvecm,1)
if options_.opt_gsa.ppost
dyn_saveas(hh,[OutDir '/' fname_ '_rmse_post_lnprior',int2str(ifig)],options_);
else
@@ -329,7 +330,7 @@ else
if options_.opt_gsa.ppost==0,
set(gca,'xlim',[min( likelihood(ixx(1:nfilt0(i),i)) ) max( likelihood(ixx(1:nfilt0(i),i)) )])
end
- if mod(i,9)==0 | i==size(vvarvecm,1)
+ if mod(i,9)==0 || i==size(vvarvecm,1)
if options_.opt_gsa.ppost
dyn_saveas(hh,[OutDir '/' fname_ '_rmse_post_lnlik',int2str(ifig) ],options_);
else
@@ -360,7 +361,7 @@ else
if options_.opt_gsa.ppost==0,
set(gca,'xlim',[min( logpo2(ixx(1:nfilt0(i),i)) ) max( logpo2(ixx(1:nfilt0(i),i)) )])
end
- if mod(i,9)==0 | i==size(vvarvecm,1)
+ if mod(i,9)==0 || i==size(vvarvecm,1)
if options_.opt_gsa.ppost
dyn_saveas(hh,[OutDir '/' fname_ '_rmse_post_lnpost',int2str(ifig) ],options_);
else
@@ -402,7 +403,7 @@ else
rmse_MC=rmse_MC(:,ivar);
disp(' ')
- % if options_.opt_gsa.ppost==0 & options_.opt_gsa.pprior,
+ % if options_.opt_gsa.ppost==0 && options_.opt_gsa.pprior,
disp(['Sample filtered the ',num2str(pfilt*100),'% best RMSE''s for each observed series ...' ])
% else
% disp(['Sample filtered the best RMSE''s smaller than RMSE at the posterior mean ...' ])
@@ -414,7 +415,7 @@ else
disp(' ')
disp(' ')
disp('RMSE ranges after filtering:')
- if options_.opt_gsa.ppost==0 & options_.opt_gsa.pprior,
+ if options_.opt_gsa.ppost==0 && options_.opt_gsa.pprior,
disp([' best ',num2str(pfilt*100),'% filtered remaining 90%'])
disp([' min max min max posterior mode'])
else
@@ -478,13 +479,18 @@ else
h0=cumplot(x(:,nsnam(j)));
set(h0,'color',[0 0 0])
hold on,
- np=find(SP(nsnam(j),:));
+ npx=find(SP(nsnam(j),:)==0);
%a0=jet(nsp(nsnam(j)));
- a0=a00(np,:);
- for i=1:nsp(nsnam(j)), %size(vvarvecm,1),
+% a0=a00(np,:);
+ for i=1:size(vvarvecm,1),
%h0=cumplot(x(ixx(1:nfilt,np(i)),nsnam(j)+nshock));
- h0=cumplot(x(ixx(1:nfilt0(np(i)),np(i)),nsnam(j)));
- set(h0,'color',a0(i,:))
+% h0=cumplot(x(ixx(1:nfilt0(np(i)),np(i)),nsnam(j)));
+ if any(npx==i),
+ h0=cumplot(x(ixx(1:nfilt0(i),i),nsnam(j))*NaN);
+ else
+ h0=cumplot(x(ixx(1:nfilt0(i),i),nsnam(j)));
+ end
+ set(h0,'color',a00(i,:))
end
ydum=get(gca,'ylim');
%xdum=xparam1(nshock+nsnam(j));
@@ -498,9 +504,9 @@ else
end
%subplot(3,2,6)
if exist('OCTAVE_VERSION'),
- legend(char('base',vvarvecm(np,:)),'location','eastoutside');
+ legend(char('base',vvarvecm),'location','eastoutside');
else
- h0=legend(char('base',vvarvecm(np,:)),0);
+ h0=legend(char('base',vvarvecm),0);
set(h0,'fontsize',6,'position',[0.7 0.1 0.2 0.3],'interpreter','none');
end
%h0=legend({'base',vnam{np}}',0);
@@ -521,104 +527,6 @@ else
nsx(j)=length(find(SP(:,j)));
end
- number_of_grid_points = 2^9; % 2^9 = 512 !... Must be a power of two.
- bandwidth = 0; % Rule of thumb optimal bandwidth parameter.
- kernel_function = 'gaussian'; % Gaussian kernel for Fast Fourrier Transform approximaton.
- %kernel_function = 'uniform'; % Gaussian kernel for Fast Fourrier Transform approximaton.
-
- for ix=1:ceil(length(nsnam)/5),
- hh = dyn_figure(options_);
- for j=1+5*(ix-1):min(size(snam2,1),5*ix),
- subplot(2,3,j-5*(ix-1))
- optimal_bandwidth = mh_optimal_bandwidth(x(:,nsnam(j)),size(x,1),bandwidth,kernel_function);
- [x1,f1] = kernel_density_estimate(x(:,nsnam(j)),number_of_grid_points,...
- size(x,1),optimal_bandwidth,kernel_function);
- h0 = plot(x1, f1,'k');
- hold on,
- np=find(SP(nsnam(j),:));
- %a0=jet(nsp(nsnam(j)));
- a0=a00(np,:);
- for i=1:nsp(nsnam(j)), %size(vvarvecm,1),
- optimal_bandwidth = mh_optimal_bandwidth(x(ixx(1:nfilt0(np(i)),np(i)),nsnam(j)),nfilt,bandwidth,kernel_function);
- [x1,f1] = kernel_density_estimate(x(ixx(1:nfilt0(np(i)),np(i)),nsnam(j)),number_of_grid_points,...
- nfilt, optimal_bandwidth,kernel_function);
- h0 = plot(x1, f1);
- set(h0,'color',a0(i,:))
- end
- ydum=get(gca,'ylim');
- set(gca,'ylim',[0 ydum(2)]);
- if exist('xparam1')
- %xdum=xparam1(nshock+nsnam(j));
- xdum=xparam1(nsnam(j));
- h1=plot([xdum xdum],[0 ydum(2)]);
- set(h1,'color',[0.85 0.85 0.85],'linewidth',2)
- end
- xlabel('')
- title([pnam{nsnam(j)}],'interpreter','none')
- end
- if exist('OCTAVE_VERSION'),
- legend(char('base',vvarvecm(np,:)),'location','eastoutside');
- else
- h0=legend(char('base',vvarvecm(np,:)),0);
- set(h0,'fontsize',6,'position',[0.7 0.1 0.2 0.3],'interpreter','none')
- end
- %h0=legend({'base',vnam{np}}',0);
- %set(findobj(get(h0,'children'),'type','text'),'interpreter','none')
- if options_.opt_gsa.ppost
- dyn_saveas(hh,[ OutDir '/' fname_ '_rmse_post_dens_' int2str(ix) ],options_);
- else
- if options_.opt_gsa.pprior
- dyn_saveas(hh,[OutDir '/' fname_ '_rmse_prior_dens_' int2str(ix)],options_);
- else
- dyn_saveas(hh,[OutDir '/' fname_ '_rmse_mc_dens_' int2str(ix) ],options_);
- end
- end
- end
- close all
- % for j=1:size(SP,2),
- % nfig=0;
- % np=find(SP(:,j));
- % for i=1:nsx(j), %size(vvarvecm,1),
- % if mod(i,12)==1,
- % nfig=nfig+1;
- % %figure('name',['Sensitivity of fit of ',vnam{j}]),
- % figure('name',['Sensitivity of fit of ',deblank(vvarvecm(j,:)),' ',num2str(nfig)]),
- % end
- %
- % subplot(3,4,i-12*(nfig-1))
- % optimal_bandwidth = mh_optimal_bandwidth(x(ixx(1:nfilt,j),np(i)),nfilt,bandwidth,kernel_function);
- % [x1,f1] = kernel_density_estimate(x(ixx(1:nfilt,j),np(i)),number_of_grid_points,...
- % nfilt, optimal_bandwidth,kernel_function);
- % plot(x1, f1,':k','linewidth',2)
- % optimal_bandwidth = mh_optimal_bandwidth(x(ixx(nfilt+1:end,j),np(i)),nruns-nfilt,bandwidth,kernel_function);
- % [x1,f1] = kernel_density_estimate(x(ixx(nfilt+1:end,j),np(i)),number_of_grid_points,...
- % nruns-nfilt,optimal_bandwidth,kernel_function);
- % hold on, plot(x1, f1,'k','linewidth',2)
- % ydum=get(gca,'ylim');
- % %xdum=xparam1(nshock+np(i));
- % xdum=xparam1(np(i));
- % h1=plot([xdum xdum],ydum);
- % set(h1,'color',[0.85 0.85 0.85],'linewidth',2)
- % %xdum1=mean(x(ixx(1:nfilt,j),np(i)+nshock));
- % xdum1=mean(x(ixx(1:nfilt,j),np(i)));
- % h2=plot([xdum1 xdum1],ydum);
- % set(h2,'color',[0 1 0],'linewidth',2)
- % % h0=cumplot(x(nfilt+1:end,np(i)+nshock));
- % % set(h0,'color',[1 1 1])
- % % hold on,
- % % h0=cumplot(x(ixx(1:nfilt,j),np(i)+nshock));
- % % set(h0,'linestyle',':','color',[1 1 1])
- % %title([pnam{np(i)}])
- % title([pnam{np(i)},'. K-S prob ', num2str(PP(np(i),j))],'interpreter','none')
- % xlabel('')
- % if mod(i,12)==0 | i==nsx(j),
- % saveas(gcf,[fname_,'_rmse_',deblank(vvarvecm(j,:)),'_',int2str(nfig)])
- % close(gcf)
- % end
- % end
- % end
-
-
disp(' ')
disp(' ')
disp('Sensitivity table (significance and direction):')
diff --git a/matlab/gsa/map_ident_.m b/matlab/gsa/map_ident_.m
index e3e6b25d7ce787eaf7ab6657fc9a1262733a09d5..452d091d7e8f67c508b72e0306d7792aae6aa2bd 100644
--- a/matlab/gsa/map_ident_.m
+++ b/matlab/gsa/map_ident_.m
@@ -1,4 +1,4 @@
-function map_ident_(OutputDirectoryName)
+function map_ident_(OutputDirectoryName,opt_gsa)
% Copyright (C) 2012 Dynare Team
%
@@ -19,7 +19,7 @@ function map_ident_(OutputDirectoryName)
global bayestopt_ M_ options_ estim_params_ oo_
-opt_gsa = options_.opt_gsa;
+% opt_gsa = options_.opt_gsa;
fname_ = M_.fname;
nliv = opt_gsa.morris_nliv;
ntra = opt_gsa.morris_ntra;
diff --git a/matlab/gsa/redform_map.m b/matlab/gsa/redform_map.m
index 2c49e940f2e0206e611f1fb6d358c4c2ca95ea53..8fed92d98f312bc84792d51f86628607c3249d4f 100644
--- a/matlab/gsa/redform_map.m
+++ b/matlab/gsa/redform_map.m
@@ -1,4 +1,4 @@
-function redform_map(dirname)
+function redform_map(dirname,options_gsa_)
%function redform_map(dirname)
% inputs (from opt_gsa structure
% anamendo = options_gsa_.namendo;
@@ -39,7 +39,7 @@ function redform_map(dirname)
global M_ oo_ estim_params_ options_ bayestopt_
-options_gsa_ = options_.opt_gsa;
+% options_gsa_ = options_.opt_gsa;
anamendo = options_gsa_.namendo;
anamlagendo = options_gsa_.namlagendo;
@@ -48,8 +48,11 @@ iload = options_gsa_.load_redform;
pprior = options_gsa_.pprior;
ilog = options_gsa_.logtrans_redform;
threshold = options_gsa_.threshold_redform;
-ksstat = options_gsa_.ksstat_redform;
+% ksstat = options_gsa_.ksstat_redform;
alpha2 = options_gsa_.alpha2_redform;
+alpha2=0;
+pvalue_ks = options_gsa_.ksstat_redform;
+pvalue_corr = options_gsa_.alpha2_redform;
pnames = M_.param_names(estim_params_.param_vals(:,1),:);
if nargin==0,
@@ -57,14 +60,14 @@ if nargin==0,
end
if pprior
- load([dirname,'/',M_.fname,'_prior']);
+ load([dirname,'/',M_.fname,'_prior'],'lpmat', 'lpmat0', 'istable','T');
adir=[dirname '/redform_stab'];
else
- load([dirname,'/',M_.fname,'_mc']);
+ load([dirname,'/',M_.fname,'_mc'],'lpmat', 'lpmat0', 'istable','T');
adir=[dirname '/redform_mc'];
end
if ~exist('T')
- stab_map_(dirname);
+ stab_map_(dirname,options_gsa_);
if pprior
load([dirname,'/',M_.fname,'_prior'],'T');
else
@@ -79,6 +82,7 @@ adir0=pwd;
nspred=size(T,2)-M_.exo_nbr;
x0=lpmat(istable,:);
+xx0=lpmat0(istable,:);
[kn, np]=size(x0);
offset = length(bayestopt_.pshape)-np;
if options_gsa_.prior_range,
@@ -90,7 +94,7 @@ else
end
nsok = length(find(M_.lead_lag_incidence(M_.maximum_lag,:)));
-clear lpmat lpmat0 egg iunstable yys
+clear lpmat lpmat0
js=0;
for j=1:size(anamendo,1)
namendo=deblank(anamendo(j,:));
@@ -105,7 +109,7 @@ for j=1:size(anamendo,1)
%y0=squeeze(T(iendo,iexo+nspred,istable));
y0=squeeze(T(iendo,iexo+nspred,:));
if (max(y0)-min(y0))>1.e-10,
- if mod(iplo,9)==0,
+ if mod(iplo,9)==0 && isempty(threshold),
ifig=ifig+1;
hfig = dyn_figure(options_,'name',[namendo,' vs. shocks ',int2str(ifig)]);
iplo=0;
@@ -115,28 +119,46 @@ for j=1:size(anamendo,1)
xdir0 = [adir,'/',namendo,'_vs_', namexo];
if ilog==0,
if isempty(threshold)
- si(:,js) = redform_private(x0, y0, pshape, pd, iload, pnames, namendo, namexo, xdir0);
+ if isempty(dir(xdir0))
+ mkdir(xdir0)
+ end
+ si(:,js) = redform_private(x0, y0, pshape, pd, iload, pnames, namendo, namexo, xdir0, options_gsa_);
else
iy=find( (y0>threshold(1)) & (y0<threshold(2)));
iyc=find( (y0<=threshold(1)) | (y0>=threshold(2)));
xdir = [xdir0,'_cut'];
- if ~isempty(iy),
- si(:,js) = redform_private(x0(iy,:), y0(iy), pshape, pd, iload, pnames, namendo, namexo, xdir);
+ if isempty(dir(xdir))
+ mkdir(xdir)
end
- if ~isempty(iy) & ~isempty(iyc)
+% if ~isempty(iy),
+% si(:,js) = redform_private(x0(iy,:), y0(iy), pshape, pd, iload, pnames, namendo, namexo, xdir, options_gsa_);
+% else
+ si(:,js) = NaN(np,1);
+% end
+ if ~isempty(iy) && ~isempty(iyc)
delete([xdir, '/*cut*.*'])
[proba, dproba] = stab_map_1(x0, iy, iyc, 'cut',0);
- indsmirnov = find(dproba>ksstat);
+% indsmirnov = find(dproba>ksstat);
+ indsmirnov = find(proba<pvalue_ks);
stab_map_1(x0, iy, iyc, 'cut',1,indsmirnov,xdir);
- stab_map_2(x0(iy,:),alpha2,'cut',xdir)
- stab_map_2(x0(iyc,:),alpha2,'trim',xdir)
+ stab_map_2(x0(iy,:),alpha2,pvalue_corr,'cut',xdir)
+ stab_map_2(x0(iyc,:),alpha2,pvalue_corr,'trim',xdir)
+ lpmat=x0(iy,:);
+ lpmat0=xx0(iy,:);
+ istable=[1:length(iy)];
+ save([xdir,filesep,'threshold.mat'],'lpmat','lpmat0','istable')
+ clear lpmat lpmat0 istable
end
end
else
[yy, xdir] = log_trans_(y0,xdir0);
- silog(:,js) = redform_private(x0, yy, pshape, pd, iload, pnames, namendo, namexo, xdir);
+ if isempty(dir(xdir))
+ mkdir(xdir)
+ end
+ silog(:,js) = redform_private(x0, yy, pshape, pd, iload, pnames, namendo, namexo, xdir, options_gsa_);
end
+ if isempty(threshold)
figure(hfig)
subplot(3,3,iplo),
if ilog,
@@ -161,11 +183,12 @@ for j=1:size(anamendo,1)
close(hfig);
end
end
+ end
end
end
end
- if iplo<9 & iplo>0 & ifig,
+ if iplo<9 && iplo>0 && ifig,
dyn_saveas(hfig,[dirname,'/',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)],options_);
if ~options_.nodisplay
close(hfig);
@@ -181,7 +204,7 @@ for j=1:size(anamendo,1)
%y0=squeeze(T(iendo,ilagendo,istable));
y0=squeeze(T(iendo,ilagendo,:));
if (max(y0)-min(y0))>1.e-10,
- if mod(iplo,9)==0,
+ if mod(iplo,9)==0 && isempty(threshold),
ifig=ifig+1;
hfig = dyn_figure(options_,'name',[namendo,' vs. lags ',int2str(ifig)]);
iplo=0;
@@ -191,28 +214,45 @@ for j=1:size(anamendo,1)
xdir0 = [adir,'/',namendo,'_vs_', namlagendo];
if ilog==0,
if isempty(threshold)
- si(:,js) = redform_private(x0, y0, pshape, pd, iload, pnames, namendo, namlagendo, xdir0);
+ if isempty(dir(xdir0))
+ mkdir(xdir0)
+ end
+ si(:,js) = redform_private(x0, y0, pshape, pd, iload, pnames, namendo, namlagendo, xdir0, options_gsa_);
else
iy=find( (y0>threshold(1)) & (y0<threshold(2)));
iyc=find( (y0<=threshold(1)) | (y0>=threshold(2)));
xdir = [xdir0,'_cut'];
- if ~isempty(iy)
- si(:,js) = redform_private(x0(iy,:), y0(iy), pshape, pd, iload, pnames, namendo, namlagendo, xdir);
+ if isempty(dir(xdir))
+ mkdir(xdir)
end
- if ~isempty(iy) & ~isempty(iyc),
+% if ~isempty(iy)
+% si(:,js) = redform_private(x0(iy,:), y0(iy), pshape, pd, iload, pnames, namendo, namlagendo, xdir, options_gsa_);
+% end
+ if ~isempty(iy) && ~isempty(iyc),
delete([xdir, '/*cut*.*'])
[proba, dproba] = stab_map_1(x0, iy, iyc, 'cut',0);
- indsmirnov = find(dproba>ksstat);
+% indsmirnov = find(dproba>ksstat);
+ indsmirnov = find(proba<pvalue_ks);
stab_map_1(x0, iy, iyc, 'cut',1,indsmirnov,xdir);
- stab_map_2(x0(iy,:),alpha2,'cut',xdir)
- stab_map_2(x0(iyc,:),alpha2,'trim',xdir)
+ stab_map_2(x0(iy,:),alpha2,pvalue_corr,'cut',xdir)
+ stab_map_2(x0(iyc,:),alpha2,pvalue_corr,'trim',xdir)
+ lpmat=x0(iy,:);
+ lpmat0=xx0(iy,:);
+ istable=[1:length(iy)];
+ save([xdir,filesep,'threshold.mat'],'lpmat','lpmat0','istable')
+ clear lpmat lpmat0 istable
+
end
end
else
[yy, xdir] = log_trans_(y0,xdir0);
- silog(:,js) = redform_private(x0, yy, pshape, pd, iload, pnames, namendo, namlagendo, xdir);
+ if isempty(dir(xdir))
+ mkdir(xdir)
+ end
+ silog(:,js) = redform_private(x0, yy, pshape, pd, iload, pnames, namendo, namlagendo, xdir, options_gsa_);
end
+ if isempty(threshold)
figure(hfig),
subplot(3,3,iplo),
if ilog,
@@ -237,11 +277,12 @@ for j=1:size(anamendo,1)
close(hfig);
end
end
+ end
end
end
end
- if iplo<9 & iplo>0 & ifig,
+ if iplo<9 && iplo>0 && ifig,
dyn_saveas(hfig,[dirname,'/',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)],options_);
if ~options_.nodisplay
close(hfig);
@@ -249,6 +290,7 @@ for j=1:size(anamendo,1)
end
end
+if isempty(threshold),
if ilog==0,
hfig=dyn_figure(options_); %bar(si)
% boxplot(si','whis',10,'symbol','r.')
@@ -280,10 +322,12 @@ title('Reduced form GSA - Log-transformed elements')
dyn_saveas(hfig,[dirname,'/',M_.fname,'_redform_gsa_log'],options_);
end
-function si = redform_private(x0, y0, pshape, pd, iload, pnames, namy, namx, xdir)
+end
+
+function si = redform_private(x0, y0, pshape, pd, iload, pnames, namy, namx, xdir, opt_gsa)
global bayestopt_ options_
-opt_gsa=options_.opt_gsa;
+% opt_gsa=options_.opt_gsa;
np=size(x0,2);
x00=x0;
if opt_gsa.prior_range,
@@ -319,9 +363,6 @@ if iload==0,
hfig=gsa_sdp_plot(gsa_,fname,pnames,iii(1:min(12,np)));
close(hfig);
gsa_.x0=x0(1:nfit,:);
- if ~options_.nodisplay
- close(hfig);
- end
% copyfile([fname,'_est.mat'],[fname,'.mat'])
hfig=dyn_figure(options_);
plot(y0(1:nfit),[gsa_.fit y0(1:nfit)],'.'),
diff --git a/matlab/gsa/redform_screen.m b/matlab/gsa/redform_screen.m
index 5be8a871eb70bc97140364f03c413c5bbe3dec45..0df22054621fc71dced1a8a37195fe1db3680eb3 100644
--- a/matlab/gsa/redform_screen.m
+++ b/matlab/gsa/redform_screen.m
@@ -1,5 +1,5 @@
-function redform_screen(dirname)
-%function redform_map(dirname)
+function redform_screen(dirname, options_gsa_)
+%function redform_map(dirname, options_gsa_)
% inputs (from opt_gsa structure
% anamendo = options_gsa_.namendo;
% anamlagendo = options_gsa_.namlagendo;
@@ -33,7 +33,7 @@ function redform_screen(dirname)
global M_ oo_ estim_params_ options_ bayestopt_
-options_gsa_ = options_.opt_gsa;
+% options_gsa_ = options_.opt_gsa;
anamendo = options_gsa_.namendo;
anamlagendo = options_gsa_.namlagendo;
diff --git a/matlab/gsa/stab_map_.m b/matlab/gsa/stab_map_.m
index d042127c59672a0ac6f3f48ef4ad16f666fa6603..dc0f369ee5c64c0b807bbec49afff5c72ec94ba5 100644
--- a/matlab/gsa/stab_map_.m
+++ b/matlab/gsa/stab_map_.m
@@ -1,4 +1,4 @@
-function x0 = stab_map_(OutputDirectoryName)
+function x0 = stab_map_(OutputDirectoryName,opt_gsa)
%
% function x0 = stab_map_(OutputDirectoryName)
%
@@ -58,7 +58,7 @@ function x0 = stab_map_(OutputDirectoryName)
%global bayestopt_ estim_params_ dr_ options_ ys_ fname_
global bayestopt_ estim_params_ options_ oo_ M_
-opt_gsa=options_.opt_gsa;
+% opt_gsa=options_.opt_gsa;
Nsam = opt_gsa.Nsam;
fload = opt_gsa.load_stab;
@@ -568,6 +568,8 @@ else
end
+xparam1=x0;
+save prior_ok xparam1;
options_.periods=opt.periods;
if isfield(opt,'nomoments'),
diff --git a/matlab/gsa/stab_map_2.m b/matlab/gsa/stab_map_2.m
index a1a72c790760bb06dd6f3551bac32010faf26eef..aecfabf81e00c5daea1dd25b85f70a770e5c5426 100644
--- a/matlab/gsa/stab_map_2.m
+++ b/matlab/gsa/stab_map_2.m
@@ -1,5 +1,5 @@
function stab_map_2(x,alpha2, pvalue, fnam, dirname,xparam1)
-% function stab_map_2(x, alpha2, pvalue, fnam, dirname)
+% function stab_map_2(x, alpha2, pvalue, fnam, dirname,xparam1)
%
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
@@ -108,10 +108,13 @@ for j=1:npar,
end
if (j==(npar)) && j2>0,
dyn_saveas(hh,[dirname,'/',fig_nam_,int2str(ifig)],options_);
+ if ~options_.nodisplay
+ close(hh);
+ end
end
end
if ifig==0,
- disp(['No correlation term >', num2str(alpha2),' found for ',fnam])
+ disp(['No correlation term with pvalue <', num2str(pvalue),' found for ',fnam])
end
%close all
diff --git a/matlab/gsa/tcrit.m b/matlab/gsa/tcrit.m
index f401ecabe60d47bc0900f0c733efc322ab703a07..2efd947e646d58c355952ba05f46dc9bfd8c33f5 100644
--- a/matlab/gsa/tcrit.m
+++ b/matlab/gsa/tcrit.m
@@ -30,6 +30,14 @@ function t_crit = tcrit(n,pval0)
if nargin==1 || isempty(pval0),
pval0=0.05;
end
+if pval0==1,
+ t_crit=0;
+ return,
+end
+if pval0==0,
+ t_crit=inf;
+ return,
+end
pval = [ 0.10 0.05 0.025 0.01 0.005 0.001];
pval0=max(pval0,min(pval));
ncol=min(find(pval<=pval0))+1;