redform_map.m 36.3 KB
Newer Older
1
function redform_map(dirname,options_gsa_)
2 3 4 5 6 7 8 9 10 11 12 13
%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;
%
14
% Written by Marco Ratto
15
% Joint Research Centre, The European Commission,
16
% marco.ratto@ec.europa.eu
17

Marco Ratto's avatar
Marco Ratto committed
18
% Copyright (C) 2012-2016 European Commission
19
% Copyright (C) 2012-2018 Dynare Team
20 21
%
% This file is part of Dynare.
22
%
23 24 25 26 27 28 29 30 31 32 33 34 35
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.

36 37 38

global M_ oo_ estim_params_ options_ bayestopt_

39
% options_gsa_ = options_.opt_gsa;
40 41 42 43 44 45 46 47

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;
48
% ksstat = options_gsa_.ksstat_redform;
49
alpha2 = options_gsa_.alpha2_redform;
50 51 52
alpha2=0;
pvalue_ks = options_gsa_.ksstat_redform;
pvalue_corr = options_gsa_.alpha2_redform;
53 54 55 56 57

np = estim_params_.np;
nshock = estim_params_.nvx + estim_params_.nvn + estim_params_.ncx + estim_params_.ncn;
pnames=cell(np,1);
pnames_tex=cell(np,1);
58
for jj=1:np
59 60 61 62 63 64 65 66 67 68
    if options_.TeX
        [param_name_temp, param_name_tex_temp]= get_the_name(nshock+jj,options_.TeX,M_,estim_params_,options_);
        pnames_tex{jj,1} = strrep(param_name_tex_temp,'$','');
        pnames{jj,1} = param_name_temp;
    else
        param_name_temp = get_the_name(nshock+jj,options_.TeX,M_,estim_params_,options_);
        pnames{jj,1} = param_name_temp;
    end
end

69
fname_ = M_.fname;
70

71
bounds = prior_bounds(bayestopt_, options_.prior_trunc);
72

73
if nargin==0
74
    dirname='';
75 76 77
end

if pprior
78
    load([dirname,filesep,M_.fname,'_prior'],'lpmat', 'lpmat0', 'istable','T');
79 80
    adir=[dirname filesep 'redform_prior'];
    type = 'prior';
81
else
82 83
    load([dirname,filesep,M_.fname,'_mc'],'lpmat', 'lpmat0', 'istable','T');
    adir=[dirname filesep 'redform_mc'];
84
    type = 'mc';
85
end
86 87 88
options_mcf.pvalue_ks = options_gsa_.ksstat_redform;
options_mcf.pvalue_corr = options_gsa_.alpha2_redform;
options_mcf.alpha2 = options_gsa_.alpha2_redform;
89
options_mcf.param_names = pnames;
90
if options_.TeX
91
    options_mcf.param_names_tex = pnames_tex;
92
end
93 94 95
options_mcf.fname_ = M_.fname;
options_mcf.OutputDirectoryName = adir;

96
if ~exist('T')
97 98
    stab_map_(dirname,options_gsa_);
    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 105 106 107
        disp('The model is too large!')
        disp('Reduced form mapping stopped!')
        return
    end
108 109
end
if isempty(dir(adir))
110
    mkdir(adir)
111 112 113 114 115 116
end
adir0=pwd;
%cd(adir)

nspred=size(T,2)-M_.exo_nbr;
x0=lpmat(istable,:);
117
if isempty(lpmat0)
118 119 120 121 122 123
    xx0=[];
    nshocks=0;
else
    xx0=lpmat0(istable,:);
    nshocks=size(xx0,2);
end
124 125
[kn, np]=size(x0);
offset = length(bayestopt_.pshape)-np;
126
if 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)];
129
else
130 131
    pshape = bayestopt_.pshape(offset+1:end);
    pd =  [bayestopt_.p6(offset+1:end) bayestopt_.p7(offset+1:end) bayestopt_.p3(offset+1:end) bayestopt_.p4(offset+1:end)];
132 133
end

134
options_map.param_names = pnames;
135 136 137
if options_.TeX
    options_map.param_names_tex = pnames_tex;
end
138 139 140 141 142 143 144 145
options_map.fname_ = M_.fname;
options_map.OutputDirectoryName = adir;
options_map.iload = iload;
options_map.log_trans = ilog;
options_map.prior_range = options_gsa_.prior_range;
options_map.pshape = pshape;
options_map.pd = pd;

146
nsok = length(find(M_.lead_lag_incidence(M_.maximum_lag,:)));
147 148
lpmat=[];
lpmat0=[];
149
js=0;
150 151 152 153 154 155 156 157
for j = 1:length(anamendo)
    namendo = anamendo{j};
    iendo = strmatch(namendo, M_.endo_names(oo_.dr.order_var), 'exact');
    ifig = 0;
    iplo = 0;
    for jx = 1:length(anamexo)
        namexo = anamexo{jx};
        iexo=strmatch(namexo, M_.exo_names, 'exact');
158
        skipline()
159
        disp(['[', namendo,' vs ',namexo,']'])
160

161

162
        if ~isempty(iexo)
163 164
            %y0=squeeze(T(iendo,iexo+nspred,istable));
            y0=squeeze(T(iendo,iexo+nspred,:));
165 166
            if (max(y0)-min(y0))>1.e-10
                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 170 171 172
                    iplo=0;
                end
                iplo=iplo+1;
                js=js+1;
173
                xdir0 = [adir,filesep,namendo,'_vs_', namexo];
174
                if ilog==0 || ~isempty(threshold)
175 176 177 178
                    if isempty(threshold)
                        if isempty(dir(xdir0))
                            mkdir(xdir0)
                        end
179 180 181 182 183 184 185 186 187
                        atitle0=['Reduced Form Mapping (ANOVA) for ',namendo,' vs ', namexo];
                        aname=[type '_' namendo '_vs_' namexo];
                        atitle=[type ' Reduced Form Mapping (ANOVA): Parameter(s) driving ',namendo,' vs ',namexo];
                        options_map.amap_name = aname;
                        options_map.amap_title = atitle;
                        options_map.figtitle = atitle0;
                        options_map.title = [namendo,' vs ', namexo];
                        options_map.OutputDirectoryName = xdir0;
                        si(:,js) = redform_private(x0, y0, options_map, options_);
188 189 190
                    else
                        iy=find( (y0>threshold(1)) & (y0<threshold(2)));
                        iyc=find( (y0<=threshold(1)) | (y0>=threshold(2)));
191
                        xdir = [xdir0,'_threshold'];
192 193 194
                        if isempty(dir(xdir))
                            mkdir(xdir)
                        end
195
                        if ~options_.nograph
196
                            hf=dyn_figure(options_.nodisplay,'name',['Reduced Form Mapping (Monte Carlo Filtering): ',namendo,' vs ', namexo]);
197 198 199 200 201 202 203 204 205 206 207 208
                            hc = cumplot(y0);
                            a=axis; delete(hc);
                            %     hist(mat_moment{ij}),
                            x1val=max(threshold(1),a(1));
                            x2val=min(threshold(2),a(2));
                            hp = patch([x1val x2val x2val x1val],a([3 3 4 4]),'b');
                            set(hp,'FaceColor', [0.7 0.8 1])
                            hold all,
                            hc = cumplot(y0);
                            set(hc,'color','k','linewidth',2)
                            hold off,
                            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 216 217 218 219 220 221 222 223 224 225 226
                        atitle0=['Reduced Form Mapping (Monte Carlo Filtering) for ',namendo,' vs ', namexo];
                        aname=[type '_' namendo '_vs_' namexo '_threshold'];
                        atitle=[type ' Reduced Form Mapping (Monte Carlo Filtering): Parameter(s) driving ',namendo,' vs ',namexo];
                        options_mcf.amcf_name = aname;
                        options_mcf.amcf_title = atitle;
                        options_mcf.beha_title = 'inside threshold';
                        options_mcf.nobeha_title = 'outside threshold';
                        options_mcf.title = atitle0;
                        options_mcf.OutputDirectoryName = xdir;
                        if ~isempty(iy) && ~isempty(iyc)
                            fprintf(['%4.1f%% of the ',type,' support matches ',atitle0,'\n'],length(iy)/length(y0)*100)
                            icheck = mcf_analysis(x0, iy, iyc, options_mcf, options_);
227

228
                            lpmat=x0(iy,:);
229
                            if nshocks
230 231
                                lpmat0=xx0(iy,:);
                            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 236 237
                            if length(iy)<=10 || length(iyc)<=10
                                icheck = [];  % do the generic plot in any case
                            end
238 239 240
                        else
                            icheck=[];
                        end
241
                        if isempty(icheck)
242 243 244 245 246 247 248 249 250 251 252 253 254 255
                            if length(iy)<=10 
                                if isempty(iy)
                                    disp(['There are NO MC samples in the desired range [' num2str(threshold) ']!'])
                                else
                                    disp(['There are TOO FEW (<=10) MC samples  in the desired range [' num2str(threshold) ']!'])
                                end
                            elseif length(iyc)<=10
                                if isempty(iyc)
                                    disp(['ALL MC samples are in the desired range [' num2str(threshold) ']!'])
                                else
                                    disp(['Almost ALL MC samples are in the desired range [' num2str(threshold) ']!'])
                                    disp('There are TOO FEW (<=10) MC samples OUTSIDE the desired range!')
                                end
                            end
256 257 258
                            atitle0=['Monte Carlo Filtering for ',namendo,' vs ', namexo];
                            options_mcf.title = atitle0;
                            indmcf = redform_mcf(y0, x0, options_mcf, options_);
259

260 261 262 263
                        end
                    end
                else
                    [yy, xdir] = log_trans_(y0,xdir0);
264 265 266 267 268 269 270 271 272
                    atitle0=['Reduced Form Mapping (ANOVA) for log-transformed ',namendo,' vs ', namexo];
                    aname=[type '_' namendo '_vs_' namexo];
                    atitle=[type ' Reduced Form Mapping (ANOVA): Parameter(s) driving ',namendo,' vs ',namexo];
                    options_map.amap_name = aname;
                    options_map.amap_title = atitle;
                    options_map.figtitle = atitle0;
                    options_map.title = ['log(' namendo ' vs ' namexo ')'];
                    options_map.OutputDirectoryName = xdir0;
                    silog(:,js) = redform_private(x0, y0, options_map, options_);
273
                end
274

275
                if isempty(threshold) && ~options_.nograph
276 277
                    figure(hfig)
                    subplot(3,3,iplo),
278
                    if ilog
279 280 281 282 283 284 285 286 287 288 289
                        [saso, iso] = sort(-silog(:,js));
                        bar([silog(iso(1:min(np,10)),js)])
                        logflag='log';
                    else
                        [saso, iso] = sort(-si(:,js));
                        bar(si(iso(1:min(np,10)),js))
                        logflag='';
                    end
                    %set(gca,'xticklabel',pnames(iso(1:min(np,10)),:),'fontsize',8)
                    set(gca,'xticklabel',' ','fontsize',10)
                    set(gca,'xlim',[0.5 10.5])
290
                    for ip=1:min(np,10)
291 292
                        text(ip,-0.02,deblank(pnames(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
                    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 298
                    end
                end
299

300 301
            else
               disp(['This entry in the shock matrix is CONSTANT = ' num2str(mean(y0),3)])
302 303
            end
        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 310
    ifig=0;
    iplo=0;
311 312 313
    for je=1:length(anamlagendo)
        namlagendo = anamlagendo{je};
        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 319
            %y0=squeeze(T(iendo,ilagendo,istable));
            y0=squeeze(T(iendo,ilagendo,:));
320 321
            if (max(y0)-min(y0))>1.e-10
                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 325 326 327
                    iplo=0;
                end
                iplo=iplo+1;
                js=js+1;
328
                xdir0 = [adir,filesep,namendo,'_vs_', namlagendo];
329
                if ilog==0 || ~isempty(threshold)
330 331 332 333
                    if isempty(threshold)
                        if isempty(dir(xdir0))
                            mkdir(xdir0)
                        end
334 335 336 337 338 339 340 341 342
                        atitle0=['Reduced Form Mapping (ANOVA) for ',namendo,' vs ', namlagendo];
                        aname=[type '_' namendo '_vs_' namlagendo];
                        atitle=[type ' Reduced Form Mapping (ANOVA): Parameter(s) driving ',namendo,' vs ',namlagendo];
                        options_map.amap_name = aname;
                        options_map.amap_title = atitle;
                        options_map.figtitle = atitle0;
                        options_map.title = [namendo,' vs ', namlagendo];
                        options_map.OutputDirectoryName = xdir0;
                        si(:,js) = redform_private(x0, y0, options_map, options_);
343 344 345
                    else
                        iy=find( (y0>threshold(1)) & (y0<threshold(2)));
                        iyc=find( (y0<=threshold(1)) | (y0>=threshold(2)));
346
                        xdir = [xdir0,'_threshold'];
347 348 349
                        if isempty(dir(xdir))
                            mkdir(xdir)
                        end
350
                        if ~options_.nograph
351
                            hf=dyn_figure(options_.nodisplay,'name',['Reduced Form Mapping (Monte Carlo Filtering): ',namendo,' vs lagged ', namlagendo]);
352 353 354 355 356 357 358 359 360 361 362 363
                            hc = cumplot(y0);
                            a=axis; delete(hc);
                            %     hist(mat_moment{ij}),
                            x1val=max(threshold(1),a(1));
                            x2val=min(threshold(2),a(2));
                            hp = patch([x1val x2val x2val x1val],a([3 3 4 4]),'b');
                            set(hp,'FaceColor', [0.7 0.8 1])
                            hold all,
                            hc = cumplot(y0);
                            set(hc,'color','k','linewidth',2)
                            hold off,
                            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 371 372 373 374 375 376 377 378 379 380 381 382
                        atitle0=['Reduced Form Mapping (Monte Carlo Filtering) for ',namendo,' vs ', namlagendo];
                        aname=[type '_' namendo '_vs_' namlagendo '_threshold'];
                        atitle=[type ' Reduced Form Mapping (Monte Carlo Filtering): Parameter(s) driving ',namendo,' vs ',namlagendo];
                        options_mcf.amcf_name = aname;
                        options_mcf.amcf_title = atitle;
                        options_mcf.beha_title = 'inside threshold';
                        options_mcf.nobeha_title = 'outside threshold';
                        options_mcf.title = atitle0;
                        options_mcf.OutputDirectoryName = xdir;
                        if ~isempty(iy) && ~isempty(iyc)

                            fprintf(['%4.1f%% of the ',type,' support matches ',atitle0,'\n'],length(iy)/length(y0)*100)
                            icheck = mcf_analysis(x0, iy, iyc, options_mcf, options_);
383

384
                            lpmat=x0(iy,:);
385
                            if nshocks
386 387
                                lpmat0=xx0(iy,:);
                            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 392 393
                            if length(iy)<=10 || length(iyc)<=10,
                                icheck = [];  % do the generic plot in any case
                            end
394

395 396 397
                        else
                            icheck = [];
                        end
398
                        if isempty(icheck)
399 400 401 402 403 404 405 406 407 408 409 410 411 412
                            if length(iy)<=10 
                                if isempty(iy)
                                    disp(['There are NO MC samples in the desired range [' num2str(threshold) ']!'])
                                else
                                    disp(['There are TOO FEW (<=10) MC samples  in the desired range [' num2str(threshold) ']!'])
                                end
                            elseif length(iyc)<=10
                                if isempty(iyc)
                                    disp(['ALL MC samples are in the desired range [' num2str(threshold) ']!'])
                                else
                                    disp(['Almost ALL MC samples are in the desired range [' num2str(threshold) ']!'])
                                    disp('There are TOO FEW (<=10) MC samples OUTSIDE the desired range!')
                                end
                            end
413 414 415
                            atitle0=['Monte Carlo Filtering for ',namendo,' vs ', namlagendo];
                            options_mcf.title = atitle0;
                            indmcf = redform_mcf(y0, x0, options_mcf, options_);
416 417 418 419
                        end
                    end
                else
                    [yy, xdir] = log_trans_(y0,xdir0);
420 421 422 423 424 425 426 427 428
                    atitle0=['Reduced Form Mapping (ANOVA) for log-transformed ',namendo,' vs ', namlagendo];
                    aname=[type '_' namendo '_vs_' namlagendo];
                    atitle=[type ' Reduced Form Mapping (ANOVA): Parameter(s) driving ',namendo,' vs ',namlagendo];
                    options_map.amap_name = aname;
                    options_map.amap_title = atitle;
                    options_map.figtitle = atitle0;
                    options_map.title = ['log(' namendo ' vs ' namlagendo ')'];
                    options_map.OutputDirectoryName = xdir0;
                    silog(:,js) = redform_private(x0, y0, options_map, options_);
429
                end
430

431 432 433
                if isempty(threshold) && ~options_.nograph
                    figure(hfig),
                    subplot(3,3,iplo),
434
                    if ilog
435 436 437 438 439 440 441 442 443 444 445
                        [saso, iso] = sort(-silog(:,js));
                        bar([silog(iso(1:min(np,10)),js)])
                        logflag='log';
                    else
                        [saso, iso] = sort(-si(:,js));
                        bar(si(iso(1:min(np,10)),js))
                        logflag='';
                    end
                    %set(gca,'xticklabel',pnames(iso(1:min(np,10)),:),'fontsize',8)
                    set(gca,'xticklabel',' ','fontsize',10)
                    set(gca,'xlim',[0.5 10.5])
446
                    for ip=1:min(np,10)
447 448
                        text(ip,-0.02,deblank(pnames(iso(ip),:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
                    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 454
                    end
                end
455

456 457
            else
               disp(['This entry in the transition matrix is CONSTANT = ' num2str(mean(y0),3)])
458
            end
459 460
        end
    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
465 466
end

467 468
if isempty(threshold) && ~options_.nograph
    if ilog==0
469
        hfig=dyn_figure(options_.nodisplay,'name','Reduced Form GSA'); %bar(si)
470
                                                                       % boxplot(si','whis',10,'symbol','r.')
471 472 473 474 475 476
        myboxplot(si',[],'.',[],10)
        xlabel(' ')
        set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
        set(gca,'xlim',[0.5 np+0.5])
        set(gca,'ylim',[0 1])
        set(gca,'position',[0.13 0.2 0.775 0.7])
477
        for ip=1:np
478 479 480
            text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
        end
        title('Reduced form GSA')
481
        dyn_saveas(hfig,[dirname,filesep,M_.fname,'_redform_gsa'],options_.nodisplay,options_.graph_format);
482 483
        create_TeX_loader(options_,[dirname,filesep,M_.fname,'_redform_gsa'],'Reduced Form GSA','redform_gsa')

484
    else
485
        hfig=dyn_figure(options_.nodisplay,'name','Reduced Form GSA'); %bar(silog)
486
                                                                       % boxplot(silog','whis',10,'symbol','r.')
487 488 489 490 491 492
        myboxplot(silog',[],'.',[],10)
        set(gca,'xticklabel',' ','fontsize',10,'xtick',[1:np])
        xlabel(' ')
        set(gca,'xlim',[0.5 np+0.5])
        set(gca,'ylim',[0 1])
        set(gca,'position',[0.13 0.2 0.775 0.7])
493
        for ip=1:np
494 495 496
            text(ip,-0.02,deblank(pnames(ip,:)),'rotation',90,'HorizontalAlignment','right','interpreter','none')
        end
        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
501 502
end

503
function si  = redform_private(x0, y0, options_map, options_)
504 505 506

np=size(x0,2);
x00=x0;
507 508 509 510 511 512
ilog = options_map.log_trans;
iload = options_map.iload;
pnames = options_map.param_names;
pd = options_map.pd;
pshape = options_map.pshape;
xdir = options_map.OutputDirectoryName;
513 514
if options_map.prior_range
    for j=1:np
515
        x0(:,j)=(x0(:,j)-pd(j,3))./(pd(j,4)-pd(j,3));
516
    end
517
else
518
    x0=priorcdf(x0,pshape, pd(:,1), pd(:,2), pd(:,3), pd(:,4));
519
end
520

521
if ilog
522 523 524 525
    fname=[xdir filesep options_map.fname_ '_' options_map.amap_name '_log'];
else
    fname=[xdir filesep options_map.fname_ '_' options_map.amap_name];
end
526
if iload==0
527 528 529 530
    if isempty(dir(xdir))
        mkdir(xdir)
    end
    nrun=length(y0);
531 532
    nest=max(50,nrun/2);
    nest=min(250,nest);
533 534
    nfit=min(1000,nrun);
    %   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 537
    [ys,is] = sort(y0);
    istep = ceil(nrun/nest);
538
    if istep>1
539 540 541 542 543
        iest = is(floor(istep/2):istep:end);
        nest = length(iest);
        irest = is(setdiff([1:nrun],[floor(istep/2):istep:nrun]));
        istep = ceil(length(irest)/(nfit-nest));
        ifit = union(iest, irest(1:istep:end));
544 545 546 547 548
    else
        warning('the number of samples is too small for ANOVA estimation')
        si=nan(np,1);
        return
    end
549
    if ~ismember(irest(end),ifit)
550 551 552
        ifit = union(ifit, irest(end));
    end
    nfit=length(ifit);
553 554 555
    %     ifit = union(iest, irest(randperm(nrun-nest,nfit-nest)));
    %     ifit = iest;
    %     nfit=nest;
556 557
    ipred = setdiff([1:nrun],ifit);

558
    if ilog
559 560
        [y1, tmp, isig, lam] = log_trans_(y0(iest));
        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 566
        if ilog
            hist(y1,30)
567
        else
568
            hist(y0,30)
569 570 571
        end
        title(options_map.title,'interpreter','none')
        subplot(222)
572
        if ilog
573 574 575
            hc = cumplot(y1);
        else
            hc = cumplot(y0);
576
        end
577 578 579 580 581
        set(hc,'color','k','linewidth',2)
        title([options_map.title ' CDF'],'interpreter','none')
    end

    gsa0 = ss_anova(y0(iest), x0(iest,:), 1);
582
    if ilog
583 584
        [gsa22, gsa1, gsax] = ss_anova_log(y1(iest), x0(iest,:), isig, lam, gsa0);
    end
585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632
    %     if (gsa1.out.bic-gsa0.out.bic) < 10,
    %         y00=y0;
    %         gsa00=gsa0;
    %         gsa0=gsa1;
    %         y0=y1;
    %         ilog=1;
    %     end
    if nfit>nest
        %         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);
        nvr =  gsa0.nvr*nest^3/nfit^3;
        nvr(gsa0.stat<2) = gsa0.nvr(gsa0.stat<2)*nest^5/nfit^5;
        gsa_ = ss_anova(y0(ifit), x0(ifit,:), 1, 0, 2, nvr);
        if ilog
            gsa0 = gsa_;
            nvr1 =  gsa1.nvr*nest^3/nfit^3;
            nvr1(gsa1.stat<2) = gsa1.nvr(gsa1.stat<2)*nest^5/nfit^5;
            nvrx =  gsax.nvr*nest^3/nfit^3;
            nvrx(gsax.stat<2) = gsax.nvr(gsax.stat<2)*nest^5/nfit^5;
            [gsa22, gsa1, gsax] = ss_anova_log(y1(ifit), x0(ifit,:), isig, lam, gsa0, [nvr1' nvrx']);
            %         gsa1 = ss_anova(y1(ifit), x0(ifit,:), 1, 0, 2, nvr);
            %         gsa2=gsa1;
            %         gsa2.y = gsa0.y;
            %         gsa2.fit = (exp(gsa1.fit)-lam)*isig;
            %         gsa2.f0 = mean(gsa2.fit);
            %         gsa2.out.SSE = sum((gsa2.fit-gsa2.y).^2);
            %         gsa2.out.bic = gsa2.out.bic-nest*log(gsa1.out.SSE)+nest*log(gsa2.out.SSE);
            %         gsa2.r2 = 1-cov(gsa2.fit-gsa2.y)/cov(gsa2.y);
            %         for j=1:np,
            %             gsa2.fs(:,j) = exp(gsa1.fs(:,j)).*mean(exp(gsa1.fit-gsa1.f(:,j)))*isig-lam*isig-gsa2.f0;
            %             gsa2.f(:,j) = exp(gsa1.f(:,j)).*mean(exp(gsa1.fit-gsa1.f(:,j)))*isig-lam*isig-gsa2.f0;
            %             gsa2.si(j) = var(gsa2.f(:,j))/var(gsa2.y);
            %         end
            %         nvr =  gsax.nvr*nest^3/nfit^3;
            %         nvr(gsax.stat<2) = gsax.nvr(gsax.stat<2)*nest^5/nfit^5;
            %         gsax = ss_anova([gsa2.y-gsa2.fit], x0(ifit,:), 1, 0, 2, nvr);
            %         gsa22=gsa2;
            %         gsa22.fit = gsa2.fit+gsax.fit;
            %         gsa22.f0 = mean(gsa22.fit);
            %         gsa22.out.SSE = sum((gsa22.fit-gsa22.y).^2);
            %         gsa22.out.bic = nest*log(gsa22.out.SSE/nest) + (gsax.out.df+gsa2.out.df-1)*log(nest);
            %         gsa22.r2 = 1-sum((gsa22.fit-gsa22.y).^2)/sum((gsa22.y-mean(gsa22.y)).^2);
            %         for j=1:np,
            %             gsa22.fs(:,j) = gsa2.fs(:,j)+gsax.fs(:,j);
            %             gsa22.f(:,j) = gsa2.f(:,j)+gsax.f(:,j);
            %             gsa22.si(j) = var(gsa22.f(:,j))/var(gsa22.y);
            %         end
            gsa_ = gsa22;
        end
633
    else
634 635 636 637 638
        if ilog
            gsa_ = gsa22;
        else
            gsa_ = gsa0;
        end
639 640 641 642
    end
    save([fname,'_map.mat'],'gsa_')
    [sidum, iii]=sort(-gsa_.si);
    gsa_.x0=x00(ifit,:);
643
    if ~options_.nograph
644 645
        hmap=gsa_sdp_plot(gsa_,[fname '_map'],pnames,iii(1:min(12,np)));
        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 651 652 653
        figure(hfig);
        subplot(223),
        plot(y0(ifit),[gsa_.fit y0(ifit)],'.'),
        r2 = gsa_.r2;
654 655 656 657 658 659 660
        %         if ilog,
        %             plot(y00(ifit),[log_trans_(gsa_.fit,'',isig,lam) y00(ifit)],'.'),
        %             r2 = 1 - cov(log_trans_(gsa_.fit,'',isig,lam)-y00(ifit))/cov(y00(ifit));
        %         else
        %             plot(y0(ifit),[gsa_.fit y0(ifit)],'.'),
        %             r2 = gsa_.r2;
        %         end
661
        title(['Learning sample fit - R2=' num2str(r2,2)],'interpreter','none')
662 663
        if nfit<nrun
            if ilog
664 665 666 667 668 669 670 671 672 673
                yf = ss_anova_fcast(x0(ipred,:), gsa1);
                yf = log_trans_(yf,'',isig,lam)+ss_anova_fcast(x0(ipred,:), gsax);
            else
                yf = ss_anova_fcast(x0(ipred,:), gsa_);
            end
            yn = y0(ipred);
            r2  = 1-cov(yf-yn)/cov(yn);
            subplot(224),
            plot(yn,[yf yn],'.'),
            title(['Out-of-sample prediction - R2=' num2str(r2,2)],'interpreter','none')
674
        end
675
        dyn_saveas(hfig,fname,options_.nodisplay,options_.graph_format);
676 677
        create_TeX_loader(options_,fname,['Out-of-sample prediction - R2=' num2str(r2,2)],'redform_gsa_log')

678 679 680
        if options_.nodisplay
            close(hmap);
        end
681
    end
682
else
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 693
        create_TeX_loader(options_,[fname '_pred'],options_map.title,[namy,' vs ', namx,' pred'])

694
    end
695 696 697 698
end
% si = gsa_.multivariate.si;
si = gsa_.si;

699 700 701 702 703 704 705 706 707 708 709 710 711
return

function gsa2 = log2level_map(gsa1, isig, lam)

nest=length(gsa1.y);
np = size(gsa1.x0,2);
gsa2=gsa1;
gsa2.y = log_trans_(gsa1.y,'',isig,lam);
gsa2.fit = (exp(gsa1.fit)-lam)*isig;
gsa2.f0 = mean(gsa2.fit);
gsa2.out.SSE = sum((gsa2.fit-gsa2.y).^2);
gsa2.out.bic = gsa2.out.bic-nest*log(gsa1.out.SSE)+nest*log(gsa2.out.SSE);
gsa2.r2 = 1-cov(gsa2.fit-gsa2.y)/cov(gsa2.y);
712
for j=1:np
713 714 715 716 717 718 719 720 721 722 723 724 725
    gsa2.fs(:,j) = exp(gsa1.fs(:,j)).*mean(exp(gsa1.fit-gsa1.f(:,j)))*isig-lam*isig-gsa2.f0;
    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);
    gsa2.f(:,j) = exp(gsa1.f(:,j)).*mean(exp(gsa1.fit-gsa1.f(:,j)))*isig-lam*isig-gsa2.f0;
    gsa2.si(j) = var(gsa2.f(:,j))/var(gsa2.y);
end

return


function [gsa22, gsa1, gsax] = ss_anova_log(y,x,isig,lam,gsa0,nvrs)

[nest, np]=size(x);

726
if nargin==6
727 728 729 730 731
    gsa1 = ss_anova(y, x, 1, 0, 2, nvrs(:,1));
else
    gsa1 = ss_anova(y, x, 1);
end
gsa2 = log2level_map(gsa1, isig, lam);
732 733
if nargin >=5 && ~isempty(gsa0)
    for j=1:np
734 735 736 737 738 739 740 741 742
        nvr2(j) = var(diff(gsa2.fs(:,j),2));
        nvr0(j) = var(diff(gsa0.fs(:,j),2));
    end
    inda = find((gsa0.stat<2)&(gsa1.stat>2));
    inda = inda(log10(nvr0(inda)./nvr2(inda))/2<0);
    gsa1.nvr(inda)=gsa1.nvr(inda).*10.^(log10(nvr0(inda)./nvr2(inda)));
    gsa1 = ss_anova(y, x, 1, 0, 2, gsa1.nvr);
    gsa2 = log2level_map(gsa1, isig, lam);
end
743
if nargin==6
744 745 746 747 748 749 750 751 752 753
    gsax = ss_anova(gsa2.y-gsa2.fit, x, 1, 0, 2, nvrs(:,2));
else
    gsax = ss_anova(gsa2.y-gsa2.fit, x, 1);
end
gsa22=gsa2;
gsa22.fit = gsa2.fit+gsax.fit;
gsa22.f0 = mean(gsa22.fit);
gsa22.out.SSE = sum((gsa22.fit-gsa22.y).^2);
gsa22.out.bic = nest*log(gsa22.out.SSE/nest) + (gsax.out.df+gsa2.out.df-1)*log(nest);
gsa22.r2 = 1-sum((gsa22.fit-gsa22.y).^2)/sum((gsa22.y-mean(gsa22.y)).^2);
754
for j=1:np
755 756 757 758 759 760 761 762 763 764
    gsa22.fs(:,j) = gsa2.fs(:,j)+gsax.fs(:,j);
    gsa22.fses(:,j) = gsax.fses(:,j);
    gsa22.f(:,j) = gsa2.f(:,j)+gsax.f(:,j);
    gsa22.si(j) = var(gsa22.f(:,j))/var(gsa22.y);
end

return

function indmcf = redform_mcf(y0, x0, options_mcf, options_)

765
hfig=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);
769 770
post_deciles = [-inf; post_deciles; inf];

771
for jt=1:10
772 773 774 775 776
    indy{jt}=find( (y0>post_deciles(jt)) & (y0<=post_deciles(jt+1)));
    leg{jt}=[int2str(jt) '-dec'];
end
[proba, dproba] = stab_map_1(x0, indy{1}, indy{end}, [],0);
indmcf=find(proba<options_mcf.pvalue_ks);
777 778 779 780 781
if isempty(indmcf)
    [tmp,jtmp] = sort(proba,2,'ascend');
    indmcf = jtmp(1);
%     indmcf = jtmp(1:min(2,length(proba)));
end
782 783 784 785 786
[tmp,jtmp] = sort(proba(indmcf),2,'ascend');
indmcf = indmcf(jtmp);
nbr_par = length(indmcf);
nrow=ceil(sqrt(nbr_par+1));
ncol=nrow;
787
if nrow*(nrow-1)>nbr_par
788 789 790 791
    ncol=nrow-1;
end

cmap = colormap(jet(10));
792
for jx=1:nbr_par
793 794
    subplot(nrow,ncol,jx)
    hold off
795
    for jt=1:10
796 797
        h=cumplot(x0(indy{jt},indmcf(jx)));
        set(h,'color', cmap(jt,:), 'linewidth', 2)
798
        hold all
799 800 801 802 803 804 805 806 807 808
    end
    title(options_mcf.param_names(indmcf(jx),:),'interpreter','none')
end
hleg = legend(leg);
aa=get(hleg,'Position');
aa(1)=1-aa(3)-0.02;
aa(2)=0.02;
set(hleg,'Position',aa);
if ~isoctave
    annotation('textbox', [0.25,0.01,0.5,0.05], ...
809 810 811 812 813
               'String', options_mcf.title, ...
               'Color','black',...
               'FontWeight','bold',...
               'interpreter','none',...
               'horizontalalignment','center');
814 815
end

816
dyn_saveas(hfig,[options_mcf.OutputDirectoryName filesep options_mcf.fname_,'_',options_mcf.amcf_name],options_.nodisplay,options_.graph_format);
817
create_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

819
return
820

821 822 823 824
function []=create_TeX_loader(options_,figpath,caption,label_name,scale_factor)
if nargin<5
    scale_factor=1;
end
825
if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format)))
826
    fidTeX = fopen([figpath '.tex'],'w');
827 828 829 830
    fprintf(fidTeX,'%% TeX eps-loader file generated by redform_map.m (Dynare).\n');
    fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
    fprintf(fidTeX,'\\begin{figure}[H]\n');
    fprintf(fidTeX,'\\centering \n');
831
    fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s}\n',scale_factor,strrep(figpath,'\','/'));
832 833 834 835 836 837
    fprintf(fidTeX,'\\caption{%s.}',caption);
    fprintf(fidTeX,'\\label{Fig:%s}\n',label_name);
    fprintf(fidTeX,'\\end{figure}\n\n');
    fprintf(fidTeX,'%% End Of TeX file. \n');
    fclose(fidTeX);
end