Skip to content

GitLab

Commit 2aaa88ca authored by Marco Ratto's avatar Marco Ratto
Browse files

Merge from the latest gsa routines for 4.2.5. 

1) options_gsa passed as function argument
2) use pvalues always to trigger Smirnov plots and correlation plots;
3) eliminated density plots in rmse analysis;
4) updated tex documentation
parent 7a795cfd
Hide whitespace changes
Inline Side-by-side
doc/gsa/gsa.tex
View file @ 2aaa88ca
......@@ -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}
......
matlab/dynare_sensitivity.m
View file @ 2aaa88ca
......@@ -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,
......
matlab/gsa/filt_mc_.m
View file @ 2aaa88ca
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):')
......
matlab/gsa/map_ident_.m
View file @ 2aaa88ca
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;
......
matlab/gsa/redform_map.m
View file @ 2aaa88ca
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)