redform_map

function redform_map(dirname)
 inputs (from opt_gsa structure
 anamendo    = options_gsa_.namendo;
 anamlagendo = options_gsa_.namlagendo;
 anamexo     = options_gsa_.namexo;
 iload       = options_gsa_.load_redform;
 pprior      = options_gsa_.pprior;
 ilog        = options_gsa_.logtrans_redform;
 threshold   = options_gsa_.threshold_redform;
 ksstat      = options_gsa_.ksstat_redform;
 alpha2      = options_gsa_.alpha2_redform;

 Written by Marco Ratto
 Joint Research Centre, The European Commission,
 (http://eemc.jrc.ec.europa.eu/),
 marco.ratto@jrc.it 

 Reference:
 M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.

0001 function redform_map(dirname)
0002 %function redform_map(dirname)
0003 % inputs (from opt_gsa structure
0004 % anamendo    = options_gsa_.namendo;
0005 % anamlagendo = options_gsa_.namlagendo;
0006 % anamexo     = options_gsa_.namexo;
0007 % iload       = options_gsa_.load_redform;
0008 % pprior      = options_gsa_.pprior;
0009 % ilog        = options_gsa_.logtrans_redform;
0010 % threshold   = options_gsa_.threshold_redform;
0011 % ksstat      = options_gsa_.ksstat_redform;
0012 % alpha2      = options_gsa_.alpha2_redform;
0013 %
0014 % Written by Marco Ratto
0015 % Joint Research Centre, The European Commission,
0016 % (http://eemc.jrc.ec.europa.eu/),
0017 % marco.ratto@jrc.it
0018 %
0019 % Reference:
0020 % M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006.
0021 
0022 % Copyright (C) 2012 Dynare Team
0023 %
0024 % This file is part of Dynare.
0025 %
0026 % Dynare is free software: you can redistribute it and/or modify
0027 % it under the terms of the GNU General Public License as published by
0028 % the Free Software Foundation, either version 3 of the License, or
0029 % (at your option) any later version.
0030 %
0031 % Dynare is distributed in the hope that it will be useful,
0032 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0033 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0034 % GNU General Public License for more details.
0035 %
0036 % You should have received a copy of the GNU General Public License
0037 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0038 
0039 
0040 global M_ oo_ estim_params_ options_ bayestopt_
0041 
0042 options_gsa_ = options_.opt_gsa;
0043 
0044 anamendo = options_gsa_.namendo;
0045 anamlagendo = options_gsa_.namlagendo;
0046 anamexo = options_gsa_.namexo;
0047 iload = options_gsa_.load_redform;
0048 pprior = options_gsa_.pprior;
0049 ilog = options_gsa_.logtrans_redform;
0050 threshold = options_gsa_.threshold_redform;
0051 ksstat = options_gsa_.ksstat_redform;
0052 alpha2 = options_gsa_.alpha2_redform;
0053 
0054 pnames = M_.param_names(estim_params_.param_vals(:,1),:);
0055 if nargin==0,
0056   dirname='';
0057 end
0058 
0059 if pprior
0060   load([dirname,'/',M_.fname,'_prior']);
0061   adir=[dirname '/redform_stab'];
0062 else
0063   load([dirname,'/',M_.fname,'_mc']);
0064   adir=[dirname '/redform_mc'];
0065 end
0066 if ~exist('T')
0067   stab_map_(dirname);
0068 if pprior
0069   load([dirname,'/',M_.fname,'_prior'],'T');
0070 else
0071   load([dirname,'/',M_.fname,'_mc'],'T');
0072 end
0073 end
0074 if isempty(dir(adir))
0075   mkdir(adir)
0076 end
0077 adir0=pwd;
0078 %cd(adir)
0079 
0080 nspred=size(T,2)-M_.exo_nbr;
0081 x0=lpmat(istable,:);
0082 [kn, np]=size(x0);
0083 offset = length(bayestopt_.pshape)-np;
0084 if options_gsa_.prior_range,
0085   pshape=5*(ones(np,1));
0086   pd =  [NaN(np,1) NaN(np,1) bayestopt_.lb(offset+1:end) bayestopt_.ub(offset+1:end)];
0087 else
0088   pshape = bayestopt_.pshape(offset+1:end);
0089   pd =  [bayestopt_.p6(offset+1:end) bayestopt_.p7(offset+1:end) bayestopt_.p3(offset+1:end) bayestopt_.p4(offset+1:end)];
0090 end
0091 
0092 nsok = length(find(M_.lead_lag_incidence(M_.maximum_lag,:)));
0093 clear lpmat lpmat0 egg iunstable yys
0094 js=0;
0095 for j=1:size(anamendo,1)
0096   namendo=deblank(anamendo(j,:));
0097   iendo=strmatch(namendo,M_.endo_names(oo_.dr.order_var,:),'exact');
0098   ifig=0;
0099   iplo=0;
0100   for jx=1:size(anamexo,1)
0101     namexo=deblank(anamexo(jx,:));
0102     iexo=strmatch(namexo,M_.exo_names,'exact');
0103 
0104     if ~isempty(iexo),
0105       %y0=squeeze(T(iendo,iexo+nspred,istable));
0106       y0=squeeze(T(iendo,iexo+nspred,:));
0107       if (max(y0)-min(y0))>1.e-10,
0108         if mod(iplo,9)==0,
0109           ifig=ifig+1;
0110           hfig = dyn_figure(options_,'name',[namendo,' vs. shocks ',int2str(ifig)]);
0111           iplo=0;
0112         end
0113         iplo=iplo+1;
0114         js=js+1;
0115         xdir0 = [adir,'/',namendo,'_vs_', namexo];
0116         if ilog==0,
0117           if isempty(threshold)
0118             si(:,js) = redform_private(x0, y0, pshape, pd, iload, pnames, namendo, namexo, xdir0);
0119           else
0120             iy=find( (y0>threshold(1)) & (y0<threshold(2)));
0121             iyc=find( (y0<=threshold(1)) | (y0>=threshold(2)));
0122             xdir = [xdir0,'_cut'];
0123             if ~isempty(iy),
0124               si(:,js) = redform_private(x0(iy,:), y0(iy), pshape, pd, iload, pnames, namendo, namexo, xdir);
0125             end
0126             if ~isempty(iy) & ~isempty(iyc)
0127             delete([xdir, '/*cut*.*'])
0128             [proba, dproba] = stab_map_1(x0, iy, iyc, 'cut',0);
0129             indsmirnov = find(dproba>ksstat);
0130             stab_map_1(x0, iy, iyc, 'cut',1,indsmirnov,xdir);
0131             stab_map_2(x0(iy,:),alpha2,'cut',xdir)
0132             stab_map_2(x0(iyc,:),alpha2,'trim',xdir)
0133             end
0134           end
0135         else
0136           [yy, xdir] = log_trans_(y0,xdir0);
0137           silog(:,js) = redform_private(x0, yy, pshape, pd, iload, pnames, namendo, namexo, xdir);
0138         end
0139 
0140         figure(hfig)
0141         subplot(3,3,iplo),
0142         if ilog,
0143           [saso, iso] = sort(-silog(:,js));
0144           bar([silog(iso(1:min(np,10)),js)])
0145           logflag='log';
0146         else
0147           [saso, iso] = sort(-si(:,js));
0148           bar(si(iso(1:min(np,10)),js))
0149           logflag='';
0150         end
0151         %set(gca,'xticklabel',pnames(iso(1:min(np,10)),:),'fontsize',8)
0152         set(gca,'xticklabel',' ','fontsize',10)
0153         set(gca,'xlim',[0.5 10.5])
0154         for ip=1:min(np,10),
0155           text(ip,-0.02,deblank(pnames(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
0156         end
0157         title([logflag,' ',namendo,' vs. ',namexo],'interpreter','none')
0158         if iplo==9,
0159           dyn_saveas(hfig,[dirname,'/',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)],options_);  
0160           if ~options_.nodisplay
0161             close(hfig);
0162           end
0163         end
0164       
0165       end
0166     end
0167   end
0168   if iplo<9 & iplo>0 & ifig,
0169     dyn_saveas(hfig,[dirname,'/',M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)],options_);
0170     if ~options_.nodisplay
0171         close(hfig);
0172     end
0173   end
0174   ifig=0;
0175   iplo=0;
0176   for je=1:size(anamlagendo,1)
0177     namlagendo=deblank(anamlagendo(je,:));
0178     ilagendo=strmatch(namlagendo,M_.endo_names(oo_.dr.order_var(oo_.dr.nstatic+1:oo_.dr.nstatic+nsok),:),'exact');
0179 
0180     if ~isempty(ilagendo),
0181       %y0=squeeze(T(iendo,ilagendo,istable));
0182       y0=squeeze(T(iendo,ilagendo,:));
0183       if (max(y0)-min(y0))>1.e-10,
0184         if mod(iplo,9)==0,
0185           ifig=ifig+1;
0186           hfig = dyn_figure(options_,'name',[namendo,' vs. lags ',int2str(ifig)]);
0187           iplo=0;
0188         end
0189         iplo=iplo+1;
0190         js=js+1;
0191         xdir0 = [adir,'/',namendo,'_vs_', namlagendo];
0192         if ilog==0,
0193         if isempty(threshold)
0194           si(:,js) = redform_private(x0, y0, pshape, pd, iload, pnames, namendo, namlagendo, xdir0);
0195         else
0196           iy=find( (y0>threshold(1)) & (y0<threshold(2)));
0197           iyc=find( (y0<=threshold(1)) | (y0>=threshold(2)));
0198           xdir = [xdir0,'_cut'];
0199           if ~isempty(iy)
0200           si(:,js) = redform_private(x0(iy,:), y0(iy), pshape, pd, iload, pnames, namendo, namlagendo, xdir);
0201           end
0202           if ~isempty(iy) & ~isempty(iyc),
0203           delete([xdir, '/*cut*.*'])
0204           [proba, dproba] = stab_map_1(x0, iy, iyc, 'cut',0);
0205           indsmirnov = find(dproba>ksstat);
0206           stab_map_1(x0, iy, iyc, 'cut',1,indsmirnov,xdir);
0207           stab_map_2(x0(iy,:),alpha2,'cut',xdir)
0208           stab_map_2(x0(iyc,:),alpha2,'trim',xdir)
0209           end
0210         end
0211         else
0212           [yy, xdir] = log_trans_(y0,xdir0);
0213           silog(:,js) = redform_private(x0, yy, pshape, pd, iload, pnames, namendo, namlagendo, xdir);
0214         end
0215 
0216         figure(hfig),
0217         subplot(3,3,iplo),
0218         if ilog,
0219           [saso, iso] = sort(-silog(:,js));
0220           bar([silog(iso(1:min(np,10)),js)])
0221           logflag='log';
0222         else
0223           [saso, iso] = sort(-si(:,js));
0224           bar(si(iso(1:min(np,10)),js))
0225           logflag='';
0226         end
0227         %set(gca,'xticklabel',pnames(iso(1:min(np,10)),:),'fontsize',8)
0228         set(gca,'xticklabel',' ','fontsize',10)
0229         set(gca,'xlim',[0.5 10.5])
0230         for ip=1:min(np,10),
0231           text(ip,-0.02,deblank(pnames(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
0232         end
0233         title([logflag,' ',namendo,' vs. ',namlagendo,'(-1)'],'interpreter','none')
0234         if iplo==9,
0235           dyn_saveas(hfig,[dirname,'/',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)],options_);  
0236             if ~options_.nodisplay
0237                 close(hfig);
0238             end
0239         end
0240       
0241       end
0242     end
0243   end
0244   if iplo<9 & iplo>0 & ifig,
0245     dyn_saveas(hfig,[dirname,'/',M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)],options_);  
0246     if ~options_.nodisplay
0247       close(hfig);
0248     end
0249   end
0250 end
0251 
0252 if ilog==0,
0253 hfig=dyn_figure(options_); %bar(si)
0254 % boxplot(si','whis',10,'symbol','r.')
0255 myboxplot(si',[],'.',[],10)
0256 xlabel(' ')
0257 set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
0258 set(gca,'xlim',[0.5 np+0.5])
0259 set(gca,'ylim',[0 1])
0260 set(gca,'position',[0.13 0.2 0.775 0.7])
0261 for ip=1:np,
0262   text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
0263 end
0264 title('Reduced form GSA')
0265 dyn_saveas(hfig,[dirname,'/',M_.fname,'_redform_gsa'],options_);
0266 
0267 else
0268 hfig=dyn_figure(options_); %bar(silog)
0269 % boxplot(silog','whis',10,'symbol','r.')
0270 myboxplot(silog',[],'.',[],10)
0271 set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
0272 xlabel(' ')
0273 set(gca,'xlim',[0.5 np+0.5])
0274 set(gca,'ylim',[0 1])
0275 set(gca,'position',[0.13 0.2 0.775 0.7])
0276 for ip=1:np,
0277   text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
0278 end
0279 title('Reduced form GSA - Log-transformed elements')
0280 dyn_saveas(hfig,[dirname,'/',M_.fname,'_redform_gsa_log'],options_);
0281 
0282 end
0283 function si  = redform_private(x0, y0, pshape, pd, iload, pnames, namy, namx, xdir)
0284 global bayestopt_ options_
0285 
0286 opt_gsa=options_.opt_gsa;
0287 np=size(x0,2);
0288 x00=x0;
0289   if opt_gsa.prior_range,
0290     for j=1:np,
0291       x0(:,j)=(x0(:,j)-pd(j,3))./(pd(j,4)-pd(j,3));
0292     end
0293   else
0294     x0=priorcdf(x0,pshape, pd(:,1), pd(:,2), pd(:,3), pd(:,4));
0295   end
0296 
0297 fname=[xdir,'/map'];
0298 if iload==0,
0299   hfig=dyn_figure(options_); hist(y0,30), title([namy,' vs. ', namx])
0300   if isempty(dir(xdir))
0301     mkdir(xdir)
0302   end
0303   dyn_saveas(hfig,[xdir,'/', namy,'_vs_', namx],options_);
0304   if ~options_.nodisplay
0305     close(hfig);
0306   end
0307 %   gsa_ = gsa_sdp_dyn(y0, x0, -2, [],[],[],1,fname, pnames);
0308   nrun=length(y0);
0309   nest=min(250,nrun);
0310   nfit=min(1000,nrun);
0311 %   dotheplots = (nfit<=nest);
0312   gsa_ = gsa_sdp(y0(1:nest), x0(1:nest,:), 2, [],[-1 -1 -1 -1 -1 0],[],0,[fname,'_est'], pnames);
0313   if nfit>nest,
0314     gsa_ = gsa_sdp(y0(1:nfit), x0(1:nfit,:), -2, gsa_.nvr*nest^3/nfit^3,[-1 -1 -1 -1 -1 0],[],0,fname, pnames);
0315   end
0316   save([fname,'.mat'],'gsa_')
0317   [sidum, iii]=sort(-gsa_.si);
0318   gsa_.x0=x00(1:nfit,:);
0319   hfig=gsa_sdp_plot(gsa_,fname,pnames,iii(1:min(12,np)));
0320   close(hfig);
0321   gsa_.x0=x0(1:nfit,:);
0322   if ~options_.nodisplay
0323     close(hfig);
0324   end
0325 %   copyfile([fname,'_est.mat'],[fname,'.mat'])
0326   hfig=dyn_figure(options_); 
0327   plot(y0(1:nfit),[gsa_.fit y0(1:nfit)],'.'), 
0328   title([namy,' vs. ', namx,' fit'])
0329   dyn_saveas(hfig,[xdir,'/', namy,'_vs_', namx,'_fit'],options_);
0330   if ~options_.nodisplay
0331    close(hfig);
0332   end
0333   if nfit<nrun,
0334     npred=[nfit+1:nrun];
0335   yf = ss_anova_fcast(x0(npred,:), gsa_);
0336   hfig=dyn_figure(options_); 
0337   plot(y0(npred),[yf y0(npred)],'.'), 
0338   title([namy,' vs. ', namx,' pred'])
0339   dyn_saveas(hfig,[xdir,'/', namy,'_vs_', namx,'_pred'],options_);
0340   if ~options_.nodisplay
0341     close(hfig);
0342   end
0343   end
0344 else
0345 %   gsa_ = gsa_sdp_dyn(y0, x0, 0, [],[],[],0,fname, pnames);
0346   gsa_ = gsa_sdp(y0, x0, 0, [],[],[],0,fname, pnames);
0347   yf = ss_anova_fcast(x0, gsa_);
0348   hfig=dyn_figure(options_); 
0349   plot(y0,[yf y0],'.'), 
0350   title([namy,' vs. ', namx,' pred'])
0351   dyn_saveas(hfig,[xdir,'/', namy,'_vs_', namx,'_pred'],options_);  
0352   if ~options_.nodisplay
0353     close(hfig);
0354   end
0355 end
0356 % si = gsa_.multivariate.si;
0357 si = gsa_.si;
0358