1function redform_map(dirname,options_gsa_) 2%function redform_map(dirname) 3% inputs (from opt_gsa structure 4% anamendo = options_gsa_.namendo; 5% anamlagendo = options_gsa_.namlagendo; 6% anamexo = options_gsa_.namexo; 7% iload = options_gsa_.load_redform; 8% pprior = options_gsa_.pprior; 9% ilog = options_gsa_.logtrans_redform; 10% threshold = options_gsa_.threshold_redform; 11% ksstat = options_gsa_.ksstat_redform; 12% alpha2 = options_gsa_.alpha2_redform; 13% 14% Written by Marco Ratto 15% Joint Research Centre, The European Commission, 16% marco.ratto@ec.europa.eu 17 18% Copyright (C) 2012-2016 European Commission 19% Copyright (C) 2012-2018 Dynare Team 20% 21% This file is part of Dynare. 22% 23% Dynare is free software: you can redistribute it and/or modify 24% it under the terms of the GNU General Public License as published by 25% the Free Software Foundation, either version 3 of the License, or 26% (at your option) any later version. 27% 28% Dynare is distributed in the hope that it will be useful, 29% but WITHOUT ANY WARRANTY; without even the implied warranty of 30% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 31% GNU General Public License for more details. 32% 33% You should have received a copy of the GNU General Public License 34% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 35 36 37global M_ oo_ estim_params_ options_ bayestopt_ 38 39% options_gsa_ = options_.opt_gsa; 40 41anamendo = options_gsa_.namendo; 42anamlagendo = options_gsa_.namlagendo; 43anamexo = options_gsa_.namexo; 44iload = options_gsa_.load_redform; 45pprior = options_gsa_.pprior; 46ilog = options_gsa_.logtrans_redform; 47threshold = options_gsa_.threshold_redform; 48% ksstat = options_gsa_.ksstat_redform; 49alpha2 = options_gsa_.alpha2_redform; 50alpha2=0; 51pvalue_ks = options_gsa_.ksstat_redform; 52pvalue_corr = options_gsa_.alpha2_redform; 53 54np = estim_params_.np; 55nshock = estim_params_.nvx + estim_params_.nvn + estim_params_.ncx + estim_params_.ncn; 56pnames=cell(np,1); 57pnames_tex=cell(np,1); 58for jj=1:np 59 if options_.TeX 60 [param_name_temp, param_name_tex_temp]= get_the_name(nshock+jj,options_.TeX,M_,estim_params_,options_); 61 pnames_tex{jj,1} = strrep(param_name_tex_temp,'$',''); 62 pnames{jj,1} = param_name_temp; 63 else 64 param_name_temp = get_the_name(nshock+jj,options_.TeX,M_,estim_params_,options_); 65 pnames{jj,1} = param_name_temp; 66 end 67end 68 69fname_ = M_.fname; 70 71bounds = prior_bounds(bayestopt_, options_.prior_trunc); 72 73if nargin==0 74 dirname=''; 75end 76 77if pprior 78 load([dirname,filesep,M_.fname,'_prior'],'lpmat', 'lpmat0', 'istable','T'); 79 adir=[dirname filesep 'redform_prior']; 80 type = 'prior'; 81else 82 load([dirname,filesep,M_.fname,'_mc'],'lpmat', 'lpmat0', 'istable','T'); 83 adir=[dirname filesep 'redform_mc']; 84 type = 'mc'; 85end 86options_mcf.pvalue_ks = options_gsa_.ksstat_redform; 87options_mcf.pvalue_corr = options_gsa_.alpha2_redform; 88options_mcf.alpha2 = options_gsa_.alpha2_redform; 89options_mcf.param_names = pnames; 90if options_.TeX 91 options_mcf.param_names_tex = pnames_tex; 92end 93options_mcf.fname_ = M_.fname; 94options_mcf.OutputDirectoryName = adir; 95 96if ~exist('T') 97 stab_map_(dirname,options_gsa_); 98 if pprior 99 load([dirname,filesep,M_.fname,'_prior'],'T'); 100 else 101 load([dirname,filesep,M_.fname,'_mc'],'T'); 102 end 103 if ~exist('T') 104 disp('The model is too large!') 105 disp('Reduced form mapping stopped!') 106 return 107 end 108end 109if isempty(dir(adir)) 110 mkdir(adir) 111end 112adir0=pwd; 113%cd(adir) 114 115nspred=size(T,2)-M_.exo_nbr; 116x0=lpmat(istable,:); 117if isempty(lpmat0) 118 xx0=[]; 119 nshocks=0; 120else 121 xx0=lpmat0(istable,:); 122 nshocks=size(xx0,2); 123end 124[kn, np]=size(x0); 125offset = length(bayestopt_.pshape)-np; 126if options_gsa_.prior_range 127 pshape=5*(ones(np,1)); 128 pd = [NaN(np,1) NaN(np,1) bounds.lb(offset+1:end) bounds.ub(offset+1:end)]; 129else 130 pshape = bayestopt_.pshape(offset+1:end); 131 pd = [bayestopt_.p6(offset+1:end) bayestopt_.p7(offset+1:end) bayestopt_.p3(offset+1:end) bayestopt_.p4(offset+1:end)]; 132end 133 134options_map.param_names = pnames; 135if options_.TeX 136 options_map.param_names_tex = pnames_tex; 137end 138options_map.fname_ = M_.fname; 139options_map.OutputDirectoryName = adir; 140options_map.iload = iload; 141options_map.log_trans = ilog; 142options_map.prior_range = options_gsa_.prior_range; 143options_map.pshape = pshape; 144options_map.pd = pd; 145 146nsok = length(find(M_.lead_lag_incidence(M_.maximum_lag,:))); 147lpmat=[]; 148lpmat0=[]; 149js=0; 150for j = 1:length(anamendo) 151 namendo = anamendo{j}; 152 iendo = strmatch(namendo, M_.endo_names(oo_.dr.order_var), 'exact'); 153 ifig = 0; 154 iplo = 0; 155 for jx = 1:length(anamexo) 156 namexo = anamexo{jx}; 157 iexo=strmatch(namexo, M_.exo_names, 'exact'); 158 skipline() 159 disp(['[', namendo,' vs ',namexo,']']) 160 161 162 if ~isempty(iexo) 163 %y0=squeeze(T(iendo,iexo+nspred,istable)); 164 y0=squeeze(T(iendo,iexo+nspred,:)); 165 if (max(y0)-min(y0))>1.e-10 166 if mod(iplo,9)==0 && isempty(threshold) && ~options_.nograph 167 ifig=ifig+1; 168 hfig = dyn_figure(options_.nodisplay,'name',['Reduced Form Mapping: ', namendo,' vs shocks ',int2str(ifig)]); 169 iplo=0; 170 end 171 iplo=iplo+1; 172 js=js+1; 173 xdir0 = [adir,filesep,namendo,'_vs_', namexo]; 174 if ilog==0 || ~isempty(threshold) 175 if isempty(threshold) 176 if isempty(dir(xdir0)) 177 mkdir(xdir0) 178 end 179 atitle0=['Reduced Form Mapping (ANOVA) for ',namendo,' vs ', namexo]; 180 aname=[type '_' namendo '_vs_' namexo]; 181 atitle=[type ' Reduced Form Mapping (ANOVA): Parameter(s) driving ',namendo,' vs ',namexo]; 182 options_map.amap_name = aname; 183 options_map.amap_title = atitle; 184 options_map.figtitle = atitle0; 185 options_map.title = [namendo,' vs ', namexo]; 186 options_map.OutputDirectoryName = xdir0; 187 si(:,js) = redform_private(x0, y0, options_map, options_); 188 else 189 iy=find( (y0>threshold(1)) & (y0<threshold(2))); 190 iyc=find( (y0<=threshold(1)) | (y0>=threshold(2))); 191 xdir = [xdir0,'_threshold']; 192 if isempty(dir(xdir)) 193 mkdir(xdir) 194 end 195 if ~options_.nograph 196 hf=dyn_figure(options_.nodisplay,'name',['Reduced Form Mapping (Monte Carlo Filtering): ',namendo,' vs ', namexo]); 197 hc = cumplot(y0); 198 a=axis; delete(hc); 199 % hist(mat_moment{ij}), 200 x1val=max(threshold(1),a(1)); 201 x2val=min(threshold(2),a(2)); 202 hp = patch([x1val x2val x2val x1val],a([3 3 4 4]),'b'); 203 set(hp,'FaceColor', [0.7 0.8 1]) 204 hold all, 205 hc = cumplot(y0); 206 set(hc,'color','k','linewidth',2) 207 hold off, 208 title([namendo,' vs ', namexo ' - threshold [' num2str(threshold(1)) ' ' num2str(threshold(2)) ']'],'interpreter','none') 209 dyn_saveas(hf,[xdir,filesep, fname_ '_' type '_' namendo,'_vs_', namexo],options_.nodisplay,options_.graph_format); 210 create_TeX_loader(options_,[xdir,filesep, fname_ '_' type '_' namendo,'_vs_', namexo],['Reduced Form Mapping (Monte Carlo Filtering): ',strrep(namendo,'_','\_'),' vs ', strrep(namexo,'_','\_')],[type '_' namendo,'_vs_', namexo]) 211 end 212 si(:,js) = NaN(np,1); 213 delete([xdir, '/*threshold*.*']) 214 215 atitle0=['Reduced Form Mapping (Monte Carlo Filtering) for ',namendo,' vs ', namexo]; 216 aname=[type '_' namendo '_vs_' namexo '_threshold']; 217 atitle=[type ' Reduced Form Mapping (Monte Carlo Filtering): Parameter(s) driving ',namendo,' vs ',namexo]; 218 options_mcf.amcf_name = aname; 219 options_mcf.amcf_title = atitle; 220 options_mcf.beha_title = 'inside threshold'; 221 options_mcf.nobeha_title = 'outside threshold'; 222 options_mcf.title = atitle0; 223 options_mcf.OutputDirectoryName = xdir; 224 if ~isempty(iy) && ~isempty(iyc) 225 fprintf(['%4.1f%% of the ',type,' support matches ',atitle0,'\n'],length(iy)/length(y0)*100) 226 icheck = mcf_analysis(x0, iy, iyc, options_mcf, options_); 227 228 lpmat=x0(iy,:); 229 if nshocks 230 lpmat0=xx0(iy,:); 231 end 232 istable=[1:length(iy)]; 233 save([xdir,filesep, fname_ '_' type '_' namendo,'_vs_', namexo '_threshold' ],'lpmat','lpmat0','istable','y0','x0','xx0','iy','iyc') 234 lpmat=[]; lpmat0=[]; istable=[]; 235 if length(iy)<=10 || length(iyc)<=10 236 icheck = []; % do the generic plot in any case 237 end 238 else 239 icheck=[]; 240 end 241 if isempty(icheck) 242 if length(iy)<=10 243 if isempty(iy) 244 disp(['There are NO MC samples in the desired range [' num2str(threshold) ']!']) 245 else 246 disp(['There are TOO FEW (<=10) MC samples in the desired range [' num2str(threshold) ']!']) 247 end 248 elseif length(iyc)<=10 249 if isempty(iyc) 250 disp(['ALL MC samples are in the desired range [' num2str(threshold) ']!']) 251 else 252 disp(['Almost ALL MC samples are in the desired range [' num2str(threshold) ']!']) 253 disp('There are TOO FEW (<=10) MC samples OUTSIDE the desired range!') 254 end 255 end 256 atitle0=['Monte Carlo Filtering for ',namendo,' vs ', namexo]; 257 options_mcf.title = atitle0; 258 indmcf = redform_mcf(y0, x0, options_mcf, options_); 259 260 end 261 end 262 else 263 [yy, xdir] = log_trans_(y0,xdir0); 264 atitle0=['Reduced Form Mapping (ANOVA) for log-transformed ',namendo,' vs ', namexo]; 265 aname=[type '_' namendo '_vs_' namexo]; 266 atitle=[type ' Reduced Form Mapping (ANOVA): Parameter(s) driving ',namendo,' vs ',namexo]; 267 options_map.amap_name = aname; 268 options_map.amap_title = atitle; 269 options_map.figtitle = atitle0; 270 options_map.title = ['log(' namendo ' vs ' namexo ')']; 271 options_map.OutputDirectoryName = xdir0; 272 silog(:,js) = redform_private(x0, y0, options_map, options_); 273 end 274 275 if isempty(threshold) && ~options_.nograph 276 figure(hfig) 277 subplot(3,3,iplo), 278 if ilog 279 [saso, iso] = sort(-silog(:,js)); 280 bar([silog(iso(1:min(np,10)),js)]) 281 logflag='log'; 282 else 283 [saso, iso] = sort(-si(:,js)); 284 bar(si(iso(1:min(np,10)),js)) 285 logflag=''; 286 end 287 %set(gca,'xticklabel',pnames(iso(1:min(np,10)),:),'fontsize',8) 288 set(gca,'xticklabel',' ','fontsize',10) 289 set(gca,'xlim',[0.5 10.5]) 290 for ip=1:min(np,10) 291 text(ip,-0.02,deblank(pnames(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none') 292 end 293 title([logflag,' ',namendo,' vs ',namexo],'interpreter','none') 294 if iplo==9 295 dyn_saveas(hfig,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)],options_.nodisplay,options_.graph_format); 296 create_TeX_loader(options_,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)],[logflag,' ',strrep(namendo,'_','\_'),' vs ',strrep(namexo,'_','\_')],['redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)],1) 297 end 298 end 299 300 else 301 disp(['This entry in the shock matrix is CONSTANT = ' num2str(mean(y0),3)]) 302 end 303 end 304 end 305 if iplo<9 && iplo>0 && ifig && ~options_.nograph 306 dyn_saveas(hfig,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)],options_.nodisplay,options_.graph_format); 307 create_TeX_loader(options_,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)],[logflag,' ',strrep(namendo,'_','\_'),' vs ',strrep(namexo,'_','\_')],['redform_', namendo,'_vs_shocks_',logflag,num2str(ifig)],options_.figures.textwidth*min(iplo/3,1)) 308 end 309 ifig=0; 310 iplo=0; 311 for je=1:length(anamlagendo) 312 namlagendo = anamlagendo{je}; 313 ilagendo=strmatch(namlagendo, M_.endo_names(oo_.dr.order_var(M_.nstatic+1:M_.nstatic+nsok)), 'exact'); 314 skipline() 315 disp(['[', namendo,' vs lagged ',namlagendo,']']) 316 317 if ~isempty(ilagendo) 318 %y0=squeeze(T(iendo,ilagendo,istable)); 319 y0=squeeze(T(iendo,ilagendo,:)); 320 if (max(y0)-min(y0))>1.e-10 321 if mod(iplo,9)==0 && isempty(threshold) && ~options_.nograph 322 ifig=ifig+1; 323 hfig = dyn_figure(options_.nodisplay,'name',['Reduced Form Mapping: ' namendo,' vs lags ',int2str(ifig)]); 324 iplo=0; 325 end 326 iplo=iplo+1; 327 js=js+1; 328 xdir0 = [adir,filesep,namendo,'_vs_', namlagendo]; 329 if ilog==0 || ~isempty(threshold) 330 if isempty(threshold) 331 if isempty(dir(xdir0)) 332 mkdir(xdir0) 333 end 334 atitle0=['Reduced Form Mapping (ANOVA) for ',namendo,' vs ', namlagendo]; 335 aname=[type '_' namendo '_vs_' namlagendo]; 336 atitle=[type ' Reduced Form Mapping (ANOVA): Parameter(s) driving ',namendo,' vs ',namlagendo]; 337 options_map.amap_name = aname; 338 options_map.amap_title = atitle; 339 options_map.figtitle = atitle0; 340 options_map.title = [namendo,' vs ', namlagendo]; 341 options_map.OutputDirectoryName = xdir0; 342 si(:,js) = redform_private(x0, y0, options_map, options_); 343 else 344 iy=find( (y0>threshold(1)) & (y0<threshold(2))); 345 iyc=find( (y0<=threshold(1)) | (y0>=threshold(2))); 346 xdir = [xdir0,'_threshold']; 347 if isempty(dir(xdir)) 348 mkdir(xdir) 349 end 350 if ~options_.nograph 351 hf=dyn_figure(options_.nodisplay,'name',['Reduced Form Mapping (Monte Carlo Filtering): ',namendo,' vs lagged ', namlagendo]); 352 hc = cumplot(y0); 353 a=axis; delete(hc); 354 % hist(mat_moment{ij}), 355 x1val=max(threshold(1),a(1)); 356 x2val=min(threshold(2),a(2)); 357 hp = patch([x1val x2val x2val x1val],a([3 3 4 4]),'b'); 358 set(hp,'FaceColor', [0.7 0.8 1]) 359 hold all, 360 hc = cumplot(y0); 361 set(hc,'color','k','linewidth',2) 362 hold off, 363 title([namendo,' vs lagged ', namlagendo ' - threshold [' num2str(threshold(1)) ' ' num2str(threshold(2)) ']'],'interpreter','none') 364 dyn_saveas(hf,[xdir,filesep, fname_ '_' type '_' namendo,'_vs_', namlagendo],options_.nodisplay,options_.graph_format); 365 create_TeX_loader(options_,[xdir,filesep, fname_ '_' type '_' namendo,'_vs_', namlagendo],['Reduced Form Mapping (Monte Carlo Filtering): ',strrep(namendo,'_','\_'),' vs lagged ', strrep(namlagendo,'_','\_')],[type '_' namendo,'_vs_', namlagendo],1) 366 end 367 368 delete([xdir, '/*threshold*.*']) 369 370 atitle0=['Reduced Form Mapping (Monte Carlo Filtering) for ',namendo,' vs ', namlagendo]; 371 aname=[type '_' namendo '_vs_' namlagendo '_threshold']; 372 atitle=[type ' Reduced Form Mapping (Monte Carlo Filtering): Parameter(s) driving ',namendo,' vs ',namlagendo]; 373 options_mcf.amcf_name = aname; 374 options_mcf.amcf_title = atitle; 375 options_mcf.beha_title = 'inside threshold'; 376 options_mcf.nobeha_title = 'outside threshold'; 377 options_mcf.title = atitle0; 378 options_mcf.OutputDirectoryName = xdir; 379 if ~isempty(iy) && ~isempty(iyc) 380 381 fprintf(['%4.1f%% of the ',type,' support matches ',atitle0,'\n'],length(iy)/length(y0)*100) 382 icheck = mcf_analysis(x0, iy, iyc, options_mcf, options_); 383 384 lpmat=x0(iy,:); 385 if nshocks 386 lpmat0=xx0(iy,:); 387 end 388 istable=[1:length(iy)]; 389 save([xdir,filesep, fname_ '_' type '_' namendo,'_vs_', namlagendo '_threshold' ],'lpmat','lpmat0','istable','y0','x0','xx0','iy','iyc') 390 lpmat=[]; lpmat0=[]; istable=[]; 391 if length(iy)<=10 || length(iyc)<=10, 392 icheck = []; % do the generic plot in any case 393 end 394 395 else 396 icheck = []; 397 end 398 if isempty(icheck) 399 if length(iy)<=10 400 if isempty(iy) 401 disp(['There are NO MC samples in the desired range [' num2str(threshold) ']!']) 402 else 403 disp(['There are TOO FEW (<=10) MC samples in the desired range [' num2str(threshold) ']!']) 404 end 405 elseif length(iyc)<=10 406 if isempty(iyc) 407 disp(['ALL MC samples are in the desired range [' num2str(threshold) ']!']) 408 else 409 disp(['Almost ALL MC samples are in the desired range [' num2str(threshold) ']!']) 410 disp('There are TOO FEW (<=10) MC samples OUTSIDE the desired range!') 411 end 412 end 413 atitle0=['Monte Carlo Filtering for ',namendo,' vs ', namlagendo]; 414 options_mcf.title = atitle0; 415 indmcf = redform_mcf(y0, x0, options_mcf, options_); 416 end 417 end 418 else 419 [yy, xdir] = log_trans_(y0,xdir0); 420 atitle0=['Reduced Form Mapping (ANOVA) for log-transformed ',namendo,' vs ', namlagendo]; 421 aname=[type '_' namendo '_vs_' namlagendo]; 422 atitle=[type ' Reduced Form Mapping (ANOVA): Parameter(s) driving ',namendo,' vs ',namlagendo]; 423 options_map.amap_name = aname; 424 options_map.amap_title = atitle; 425 options_map.figtitle = atitle0; 426 options_map.title = ['log(' namendo ' vs ' namlagendo ')']; 427 options_map.OutputDirectoryName = xdir0; 428 silog(:,js) = redform_private(x0, y0, options_map, options_); 429 end 430 431 if isempty(threshold) && ~options_.nograph 432 figure(hfig), 433 subplot(3,3,iplo), 434 if ilog 435 [saso, iso] = sort(-silog(:,js)); 436 bar([silog(iso(1:min(np,10)),js)]) 437 logflag='log'; 438 else 439 [saso, iso] = sort(-si(:,js)); 440 bar(si(iso(1:min(np,10)),js)) 441 logflag=''; 442 end 443 %set(gca,'xticklabel',pnames(iso(1:min(np,10)),:),'fontsize',8) 444 set(gca,'xticklabel',' ','fontsize',10) 445 set(gca,'xlim',[0.5 10.5]) 446 for ip=1:min(np,10) 447 text(ip,-0.02,deblank(pnames(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none') 448 end 449 title([logflag,' ',namendo,' vs ',namlagendo,'(-1)'],'interpreter','none') 450 if iplo==9 451 dyn_saveas(hfig,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)],options_.nodisplay,options_.graph_format); 452 create_TeX_loader(options_,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)],[logflag,' ',strrep(namendo,'_','\_'),' vs ',strrep(namlagendo,'_','\_'),'(-1)'],['redform_', namendo,'_vs_lags_',logflag,':',num2str(ifig)],1) 453 end 454 end 455 456 else 457 disp(['This entry in the transition matrix is CONSTANT = ' num2str(mean(y0),3)]) 458 end 459 end 460 end 461 if iplo<9 && iplo>0 && ifig && ~options_.nograph 462 dyn_saveas(hfig,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)],options_.nodisplay,options_.graph_format); 463 create_TeX_loader(options_,[dirname,filesep,M_.fname,'_redform_', namendo,'_vs_lags_',logflag,num2str(ifig)],[logflag,' ',strrep(namendo,'_','\_'),' vs ',strrep(namlagendo,'_','\_'),'(-1)'],['redform_', namendo,'_vs_lags_',logflag,':',num2str(ifig)],options_.figures.textwidth*min(iplo/3,1)); 464 end 465end 466 467if isempty(threshold) && ~options_.nograph 468 if ilog==0 469 hfig=dyn_figure(options_.nodisplay,'name','Reduced Form GSA'); %bar(si) 470 % boxplot(si','whis',10,'symbol','r.') 471 myboxplot(si',[],'.',[],10) 472 xlabel(' ') 473 set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np]) 474 set(gca,'xlim',[0.5 np+0.5]) 475 set(gca,'ylim',[0 1]) 476 set(gca,'position',[0.13 0.2 0.775 0.7]) 477 for ip=1:np 478 text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none') 479 end 480 title('Reduced form GSA') 481 dyn_saveas(hfig,[dirname,filesep,M_.fname,'_redform_gsa'],options_.nodisplay,options_.graph_format); 482 create_TeX_loader(options_,[dirname,filesep,M_.fname,'_redform_gsa'],'Reduced Form GSA','redform_gsa') 483 484 else 485 hfig=dyn_figure(options_.nodisplay,'name','Reduced Form GSA'); %bar(silog) 486 % boxplot(silog','whis',10,'symbol','r.') 487 myboxplot(silog',[],'.',[],10) 488 set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np]) 489 xlabel(' ') 490 set(gca,'xlim',[0.5 np+0.5]) 491 set(gca,'ylim',[0 1]) 492 set(gca,'position',[0.13 0.2 0.775 0.7]) 493 for ip=1:np 494 text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none') 495 end 496 title('Reduced form GSA - Log-transformed elements') 497 dyn_saveas(hfig,[dirname,filesep,M_.fname,'_redform_gsa_log'],options_.nodisplay,options_.graph_format); 498 create_TeX_loader(options_,[dirname,filesep,M_.fname,'_redform_gsa_log'],'Reduced form GSA - Log-transformed elements','redform_gsa_log') 499 500 end 501end 502 503function si = redform_private(x0, y0, options_map, options_) 504 505np=size(x0,2); 506x00=x0; 507ilog = options_map.log_trans; 508iload = options_map.iload; 509pnames = options_map.param_names; 510pd = options_map.pd; 511pshape = options_map.pshape; 512xdir = options_map.OutputDirectoryName; 513if options_map.prior_range 514 for j=1:np 515 x0(:,j)=(x0(:,j)-pd(j,3))./(pd(j,4)-pd(j,3)); 516 end 517else 518 x0=priorcdf(x0,pshape, pd(:,1), pd(:,2), pd(:,3), pd(:,4)); 519end 520 521if ilog 522 fname=[xdir filesep options_map.fname_ '_' options_map.amap_name '_log']; 523else 524 fname=[xdir filesep options_map.fname_ '_' options_map.amap_name]; 525end 526if iload==0 527 if isempty(dir(xdir)) 528 mkdir(xdir) 529 end 530 nrun=length(y0); 531 nest=max(50,nrun/2); 532 nest=min(250,nest); 533 nfit=min(1000,nrun); 534 % dotheplots = (nfit<=nest); 535 % gsa_ = gsa_sdp(y0(1:nest), x0(1:nest,:), 2, [],[-1 -1 -1 -1 -1 0],[],0,[fname,'_est'], pnames); 536 [ys,is] = sort(y0); 537 istep = ceil(nrun/nest); 538 if istep>1 539 iest = is(floor(istep/2):istep:end); 540 nest = length(iest); 541 irest = is(setdiff([1:nrun],[floor(istep/2):istep:nrun])); 542 istep = ceil(length(irest)/(nfit-nest)); 543 ifit = union(iest, irest(1:istep:end)); 544 else 545 warning('the number of samples is too small for ANOVA estimation') 546 si=nan(np,1); 547 return 548 end 549 if ~ismember(irest(end),ifit) 550 ifit = union(ifit, irest(end)); 551 end 552 nfit=length(ifit); 553 % ifit = union(iest, irest(randperm(nrun-nest,nfit-nest))); 554 % ifit = iest; 555 % nfit=nest; 556 ipred = setdiff([1:nrun],ifit); 557 558 if ilog 559 [y1, tmp, isig, lam] = log_trans_(y0(iest)); 560 y1 = log(y0*isig+lam); 561 end 562 if ~options_.nograph 563 hfig=dyn_figure(options_.nodisplay,'name',options_map.figtitle); 564 subplot(221) 565 if ilog 566 hist(y1,30) 567 else 568 hist(y0,30) 569 end 570 title(options_map.title,'interpreter','none') 571 subplot(222) 572 if ilog 573 hc = cumplot(y1); 574 else 575 hc = cumplot(y0); 576 end 577 set(hc,'color','k','linewidth',2) 578 title([options_map.title ' CDF'],'interpreter','none') 579 end 580 581 gsa0 = ss_anova(y0(iest), x0(iest,:), 1); 582 if ilog 583 [gsa22, gsa1, gsax] = ss_anova_log(y1(iest), x0(iest,:), isig, lam, gsa0); 584 end 585 % if (gsa1.out.bic-gsa0.out.bic) < 10, 586 % y00=y0; 587 % gsa00=gsa0; 588 % gsa0=gsa1; 589 % y0=y1; 590 % ilog=1; 591 % end 592 if nfit>nest 593 % 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); 594 nvr = gsa0.nvr*nest^3/nfit^3; 595 nvr(gsa0.stat<2) = gsa0.nvr(gsa0.stat<2)*nest^5/nfit^5; 596 gsa_ = ss_anova(y0(ifit), x0(ifit,:), 1, 0, 2, nvr); 597 if ilog 598 gsa0 = gsa_; 599 nvr1 = gsa1.nvr*nest^3/nfit^3; 600 nvr1(gsa1.stat<2) = gsa1.nvr(gsa1.stat<2)*nest^5/nfit^5; 601 nvrx = gsax.nvr*nest^3/nfit^3; 602 nvrx(gsax.stat<2) = gsax.nvr(gsax.stat<2)*nest^5/nfit^5; 603 [gsa22, gsa1, gsax] = ss_anova_log(y1(ifit), x0(ifit,:), isig, lam, gsa0, [nvr1' nvrx']); 604 % gsa1 = ss_anova(y1(ifit), x0(ifit,:), 1, 0, 2, nvr); 605 % gsa2=gsa1; 606 % gsa2.y = gsa0.y; 607 % gsa2.fit = (exp(gsa1.fit)-lam)*isig; 608 % gsa2.f0 = mean(gsa2.fit); 609 % gsa2.out.SSE = sum((gsa2.fit-gsa2.y).^2); 610 % gsa2.out.bic = gsa2.out.bic-nest*log(gsa1.out.SSE)+nest*log(gsa2.out.SSE); 611 % gsa2.r2 = 1-cov(gsa2.fit-gsa2.y)/cov(gsa2.y); 612 % for j=1:np, 613 % gsa2.fs(:,j) = exp(gsa1.fs(:,j)).*mean(exp(gsa1.fit-gsa1.f(:,j)))*isig-lam*isig-gsa2.f0; 614 % gsa2.f(:,j) = exp(gsa1.f(:,j)).*mean(exp(gsa1.fit-gsa1.f(:,j)))*isig-lam*isig-gsa2.f0; 615 % gsa2.si(j) = var(gsa2.f(:,j))/var(gsa2.y); 616 % end 617 % nvr = gsax.nvr*nest^3/nfit^3; 618 % nvr(gsax.stat<2) = gsax.nvr(gsax.stat<2)*nest^5/nfit^5; 619 % gsax = ss_anova([gsa2.y-gsa2.fit], x0(ifit,:), 1, 0, 2, nvr); 620 % gsa22=gsa2; 621 % gsa22.fit = gsa2.fit+gsax.fit; 622 % gsa22.f0 = mean(gsa22.fit); 623 % gsa22.out.SSE = sum((gsa22.fit-gsa22.y).^2); 624 % gsa22.out.bic = nest*log(gsa22.out.SSE/nest) + (gsax.out.df+gsa2.out.df-1)*log(nest); 625 % gsa22.r2 = 1-sum((gsa22.fit-gsa22.y).^2)/sum((gsa22.y-mean(gsa22.y)).^2); 626 % for j=1:np, 627 % gsa22.fs(:,j) = gsa2.fs(:,j)+gsax.fs(:,j); 628 % gsa22.f(:,j) = gsa2.f(:,j)+gsax.f(:,j); 629 % gsa22.si(j) = var(gsa22.f(:,j))/var(gsa22.y); 630 % end 631 gsa_ = gsa22; 632 end 633 else 634 if ilog 635 gsa_ = gsa22; 636 else 637 gsa_ = gsa0; 638 end 639 end 640 save([fname,'_map.mat'],'gsa_') 641 [sidum, iii]=sort(-gsa_.si); 642 gsa_.x0=x00(ifit,:); 643 if ~options_.nograph 644 hmap=gsa_sdp_plot(gsa_,[fname '_map'],pnames,iii(1:min(12,np))); 645 set(hmap,'name',options_map.amap_title); 646 end 647 gsa_.x0=x0(ifit,:); 648 % copyfile([fname,'_est.mat'],[fname,'.mat']) 649 if ~options_.nograph 650 figure(hfig); 651 subplot(223), 652 plot(y0(ifit),[gsa_.fit y0(ifit)],'.'), 653 r2 = gsa_.r2; 654 % if ilog, 655 % plot(y00(ifit),[log_trans_(gsa_.fit,'',isig,lam) y00(ifit)],'.'), 656 % r2 = 1 - cov(log_trans_(gsa_.fit,'',isig,lam)-y00(ifit))/cov(y00(ifit)); 657 % else 658 % plot(y0(ifit),[gsa_.fit y0(ifit)],'.'), 659 % r2 = gsa_.r2; 660 % end 661 title(['Learning sample fit - R2=' num2str(r2,2)],'interpreter','none') 662 if nfit<nrun 663 if ilog 664 yf = ss_anova_fcast(x0(ipred,:), gsa1); 665 yf = log_trans_(yf,'',isig,lam)+ss_anova_fcast(x0(ipred,:), gsax); 666 else 667 yf = ss_anova_fcast(x0(ipred,:), gsa_); 668 end 669 yn = y0(ipred); 670 r2 = 1-cov(yf-yn)/cov(yn); 671 subplot(224), 672 plot(yn,[yf yn],'.'), 673 title(['Out-of-sample prediction - R2=' num2str(r2,2)],'interpreter','none') 674 end 675 dyn_saveas(hfig,fname,options_.nodisplay,options_.graph_format); 676 create_TeX_loader(options_,fname,['Out-of-sample prediction - R2=' num2str(r2,2)],'redform_gsa_log') 677 678 if options_.nodisplay 679 close(hmap); 680 end 681 end 682else 683 % gsa_ = gsa_sdp_dyn(y0, x0, 0, [],[],[],0,fname, pnames); 684 % gsa_ = gsa_sdp(y0, x0, 0, [],[],[],0,fname, pnames); 685 load([fname,'_map.mat'],'gsa_') 686 if ~options_.nograph 687 yf = ss_anova_fcast(x0, gsa_); 688 hfig=dyn_figure(options_.nodisplay,'name',options_map.title); 689 plot(y0,[yf y0],'.'), 690 title([namy,' vs ', namx,' pred'],'interpreter','none') 691 dyn_saveas(hfig,[fname '_pred'],options_.nodisplay,options_.graph_format); 692 create_TeX_loader(options_,[fname '_pred'],options_map.title,[namy,' vs ', namx,' pred']) 693 694 end 695end 696% si = gsa_.multivariate.si; 697si = gsa_.si; 698 699return 700 701function gsa2 = log2level_map(gsa1, isig, lam) 702 703nest=length(gsa1.y); 704np = size(gsa1.x0,2); 705gsa2=gsa1; 706gsa2.y = log_trans_(gsa1.y,'',isig,lam); 707gsa2.fit = (exp(gsa1.fit)-lam)*isig; 708gsa2.f0 = mean(gsa2.fit); 709gsa2.out.SSE = sum((gsa2.fit-gsa2.y).^2); 710gsa2.out.bic = gsa2.out.bic-nest*log(gsa1.out.SSE)+nest*log(gsa2.out.SSE); 711gsa2.r2 = 1-cov(gsa2.fit-gsa2.y)/cov(gsa2.y); 712for j=1:np 713 gsa2.fs(:,j) = exp(gsa1.fs(:,j)).*mean(exp(gsa1.fit-gsa1.f(:,j)))*isig-lam*isig-gsa2.f0; 714 gsa2.fses(:,j) = exp(gsa1.fs(:,j)+gsa1.fses(:,j)).*mean(exp(gsa1.fit-gsa1.f(:,j)))*isig-lam*isig-gsa2.f0-gsa2.fs(:,j); 715 gsa2.f(:,j) = exp(gsa1.f(:,j)).*mean(exp(gsa1.fit-gsa1.f(:,j)))*isig-lam*isig-gsa2.f0; 716 gsa2.si(j) = var(gsa2.f(:,j))/var(gsa2.y); 717end 718 719return 720 721 722function [gsa22, gsa1, gsax] = ss_anova_log(y,x,isig,lam,gsa0,nvrs) 723 724[nest, np]=size(x); 725 726if nargin==6 727 gsa1 = ss_anova(y, x, 1, 0, 2, nvrs(:,1)); 728else 729 gsa1 = ss_anova(y, x, 1); 730end 731gsa2 = log2level_map(gsa1, isig, lam); 732if nargin >=5 && ~isempty(gsa0) 733 for j=1:np 734 nvr2(j) = var(diff(gsa2.fs(:,j),2)); 735 nvr0(j) = var(diff(gsa0.fs(:,j),2)); 736 end 737 inda = find((gsa0.stat<2)&(gsa1.stat>2)); 738 inda = inda(log10(nvr0(inda)./nvr2(inda))/2<0); 739 gsa1.nvr(inda)=gsa1.nvr(inda).*10.^(log10(nvr0(inda)./nvr2(inda))); 740 gsa1 = ss_anova(y, x, 1, 0, 2, gsa1.nvr); 741 gsa2 = log2level_map(gsa1, isig, lam); 742end 743if nargin==6 744 gsax = ss_anova(gsa2.y-gsa2.fit, x, 1, 0, 2, nvrs(:,2)); 745else 746 gsax = ss_anova(gsa2.y-gsa2.fit, x, 1); 747end 748gsa22=gsa2; 749gsa22.fit = gsa2.fit+gsax.fit; 750gsa22.f0 = mean(gsa22.fit); 751gsa22.out.SSE = sum((gsa22.fit-gsa22.y).^2); 752gsa22.out.bic = nest*log(gsa22.out.SSE/nest) + (gsax.out.df+gsa2.out.df-1)*log(nest); 753gsa22.r2 = 1-sum((gsa22.fit-gsa22.y).^2)/sum((gsa22.y-mean(gsa22.y)).^2); 754for j=1:np 755 gsa22.fs(:,j) = gsa2.fs(:,j)+gsax.fs(:,j); 756 gsa22.fses(:,j) = gsax.fses(:,j); 757 gsa22.f(:,j) = gsa2.f(:,j)+gsax.f(:,j); 758 gsa22.si(j) = var(gsa22.f(:,j))/var(gsa22.y); 759end 760 761return 762 763function indmcf = redform_mcf(y0, x0, options_mcf, options_) 764 765hfig=dyn_figure(options_.nodisplay,'name',options_mcf.amcf_title); 766 767[post_mean, post_median, post_var, hpd_interval, post_deciles, ... 768 density] = posterior_moments(y0,1,0.9); 769post_deciles = [-inf; post_deciles; inf]; 770 771for jt=1:10 772 indy{jt}=find( (y0>post_deciles(jt)) & (y0<=post_deciles(jt+1))); 773 leg{jt}=[int2str(jt) '-dec']; 774end 775[proba, dproba] = stab_map_1(x0, indy{1}, indy{end}, [],0); 776indmcf=find(proba<options_mcf.pvalue_ks); 777if isempty(indmcf) 778 [tmp,jtmp] = sort(proba,2,'ascend'); 779 indmcf = jtmp(1); 780% indmcf = jtmp(1:min(2,length(proba))); 781end 782[tmp,jtmp] = sort(proba(indmcf),2,'ascend'); 783indmcf = indmcf(jtmp); 784nbr_par = length(indmcf); 785nrow=ceil(sqrt(nbr_par+1)); 786ncol=nrow; 787if nrow*(nrow-1)>nbr_par 788 ncol=nrow-1; 789end 790 791cmap = colormap(jet(10)); 792for jx=1:nbr_par 793 subplot(nrow,ncol,jx) 794 hold off 795 for jt=1:10 796 h=cumplot(x0(indy{jt},indmcf(jx))); 797 set(h,'color', cmap(jt,:), 'linewidth', 2) 798 hold all 799 end 800 title(options_mcf.param_names(indmcf(jx),:),'interpreter','none') 801end 802hleg = legend(leg); 803aa=get(hleg,'Position'); 804aa(1)=1-aa(3)-0.02; 805aa(2)=0.02; 806set(hleg,'Position',aa); 807if ~isoctave 808 annotation('textbox', [0.25,0.01,0.5,0.05], ... 809 'String', options_mcf.title, ... 810 'Color','black',... 811 'FontWeight','bold',... 812 'interpreter','none',... 813 'horizontalalignment','center'); 814end 815 816dyn_saveas(hfig,[options_mcf.OutputDirectoryName filesep options_mcf.fname_,'_',options_mcf.amcf_name],options_.nodisplay,options_.graph_format); 817create_TeX_loader(options_,[options_mcf.OutputDirectoryName filesep options_mcf.fname_,'_',options_mcf.amcf_name],strrep(options_mcf.amcf_title,'_','\_'),[options_mcf.fname_,'_',options_mcf.amcf_name]) 818 819return 820 821function []=create_TeX_loader(options_,figpath,caption,label_name,scale_factor) 822if nargin<5 823 scale_factor=1; 824end 825if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format))) 826 fidTeX = fopen([figpath '.tex'],'w'); 827 fprintf(fidTeX,'%% TeX eps-loader file generated by redform_map.m (Dynare).\n'); 828 fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']); 829 fprintf(fidTeX,'\\begin{figure}[H]\n'); 830 fprintf(fidTeX,'\\centering \n'); 831 fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s}\n',scale_factor,strrep(figpath,'\','/')); 832 fprintf(fidTeX,'\\caption{%s.}',caption); 833 fprintf(fidTeX,'\\label{Fig:%s}\n',label_name); 834 fprintf(fidTeX,'\\end{figure}\n\n'); 835 fprintf(fidTeX,'%% End Of TeX file. \n'); 836 fclose(fidTeX); 837end 838