map_calibration.m 26.4 KB
Newer Older
1
function map_calibration(OutputDirectoryName, Model, DynareOptions, DynareResults, EstimatedParameters, BayesInfo)
2

3 4
% Written by Marco Ratto
% Joint Research Centre, The European Commission,
5
% marco.ratto@ec.europa.eu
6 7

% Copyright (C) 2014-2016 European Commission
8
% Copyright (C) 2014-2017 Dynare Team
9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
%
% This file is part of Dynare.
%
% 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/>.

fname_ = Model.fname;
26 27 28 29 30

np = EstimatedParameters.np;
nshock = EstimatedParameters.nvx + EstimatedParameters.nvn + EstimatedParameters.ncx + EstimatedParameters.ncn;
pnames=cell(np,1);
pnames_tex=cell(np,1);
31
for jj=1:np
32 33
    if DynareOptions.TeX
        [param_name_temp, param_name_tex_temp]= get_the_name(nshock+jj, DynareOptions.TeX, Model, EstimatedParameters, DynareOptions);
34 35 36
        pnames_tex{jj,1} = strrep(param_name_tex_temp,'$','');
        pnames{jj,1} = param_name_temp;
    else
37
        param_name_temp = get_the_name(nshock+jj, DynareOptions.TeX, Model, EstimatedParameters, DynareOptions);
38 39 40 41
        pnames{jj,1} = param_name_temp;
    end
end

42 43 44
pvalue_ks = DynareOptions.opt_gsa.pvalue_ks;
indx_irf = [];
indx_moment = [];
45 46 47 48 49
init = ~DynareOptions.opt_gsa.load_stab;

options_mcf.pvalue_ks = DynareOptions.opt_gsa.pvalue_ks;
options_mcf.pvalue_corr = DynareOptions.opt_gsa.pvalue_corr;
options_mcf.alpha2 = DynareOptions.opt_gsa.alpha2_stab;
50
options_mcf.param_names = char(pnames);
51
if DynareOptions.TeX
52 53
    options_mcf.param_names_tex=char(pnames_tex);
end
54 55 56
options_mcf.fname_ = fname_;
options_mcf.OutputDirectoryName = OutputDirectoryName;

57 58 59
skipline()
disp('Sensitivity analysis for calibration criteria')

60
if DynareOptions.opt_gsa.ppost
61 62
    filetoload=dir([Model.dname filesep 'metropolis' filesep fname_ '_param_irf*.mat']);
    lpmat=[];
63
    for j=1:length(filetoload)
64 65 66 67 68
        load([Model.dname filesep 'metropolis' filesep fname_ '_param_irf',int2str(j),'.mat'])
        lpmat = [lpmat; stock];
        clear stock
    end
    type = 'post';
69
else
70 71 72 73 74 75 76 77 78 79
    if DynareOptions.opt_gsa.pprior
        filetoload=[OutputDirectoryName '/' fname_ '_prior'];
        load(filetoload,'lpmat','lpmat0','istable','iunstable','iindeterm','iwrong' ,'infox')
        lpmat = [lpmat0 lpmat];
        type = 'prior';
    else
        filetoload=[OutputDirectoryName '/' fname_ '_mc'];
        load(filetoload,'lpmat','lpmat0','istable','iunstable','iindeterm','iwrong' ,'infox')
        lpmat = [lpmat0 lpmat];
        type = 'mc';
80
    end
81 82
end
[Nsam, np] = size(lpmat);
83 84
npar = size(pnames,1);
nshock = np - npar;
85 86

nbr_irf_restrictions = size(DynareOptions.endogenous_prior_restrictions.irf,1);
87 88 89
nbr_moment_restrictions = size(DynareOptions.endogenous_prior_restrictions.moment,1);

if init
90
    mat_irf=cell(nbr_irf_restrictions,1);
91
    for ij=1:nbr_irf_restrictions
92 93
        mat_irf{ij}=NaN(Nsam,length(DynareOptions.endogenous_prior_restrictions.irf{ij,3}));
    end
94

95
    mat_moment=cell(nbr_moment_restrictions,1);
96
    for ij=1:nbr_moment_restrictions
97 98
        mat_moment{ij}=NaN(Nsam,length(DynareOptions.endogenous_prior_restrictions.moment{ij,3}));
    end
99

100 101
    irestrictions = [1:Nsam];
    h = dyn_waitbar(0,'Please wait...');
102
    for j=1:Nsam
103
        Model = set_all_parameters(lpmat(j,:)',EstimatedParameters,Model);
104
        if nbr_moment_restrictions
105 106 107
            [Tt,Rr,SteadyState,info,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults);
        else
            [Tt,Rr,SteadyState,info,Model,DynareOptions,DynareResults] = dynare_resolve(Model,DynareOptions,DynareResults,'restrict');
108
        end
109
        if info(1)==0
110 111
            [info, info_irf, info_moment, data_irf, data_moment]=endogenous_prior_restrictions(Tt,Rr,Model,DynareOptions,DynareResults);
            if ~isempty(info_irf)
112
                for ij=1:nbr_irf_restrictions
113 114 115
                    mat_irf{ij}(j,:)=data_irf{ij}(:,2)';
                end
                indx_irf(j,:)=info_irf(:,1);
116
            end
117
            if ~isempty(info_moment)
118
                for ij=1:nbr_moment_restrictions
119 120 121 122 123 124
                    mat_moment{ij}(j,:)=data_moment{ij}(:,2)';
                end
                indx_moment(j,:)=info_moment(:,1);
            end
        else
            irestrictions(j)=0;
125
        end
126
        dyn_waitbar(j/Nsam,h,['MC iteration ',int2str(j),'/',int2str(Nsam)])
127
    end
128
    dyn_waitbar_close(h);
129

130 131 132 133 134
    irestrictions=irestrictions(find(irestrictions));
    xmat=lpmat(irestrictions,:);
    skipline()
    endo_prior_restrictions=DynareOptions.endogenous_prior_restrictions;
    save([OutputDirectoryName,filesep,fname_,'_',type,'_restrictions'],'xmat','mat_irf','mat_moment','irestrictions','indx_irf','indx_moment','endo_prior_restrictions');
135
else
136 137
    load([OutputDirectoryName,filesep,fname_,'_',type,'_restrictions'],'xmat','mat_irf','mat_moment','irestrictions','indx_irf','indx_moment','endo_prior_restrictions');
end
138
if ~isempty(indx_irf)
139 140 141
    skipline()
    disp('Deleting old IRF calibration plots ...')
    a=dir([OutputDirectoryName,filesep,fname_,'_',type,'_irf_calib*.eps']);
142
    for j=1:length(a)
143 144 145
        delete([OutputDirectoryName,filesep,a(j).name]);
    end
    a=dir([OutputDirectoryName,filesep,fname_,'_',type,'_irf_calib*.fig']);
146
    for j=1:length(a)
147 148 149
        delete([OutputDirectoryName,filesep,a(j).name]);
    end
    a=dir([OutputDirectoryName,filesep,fname_,'_',type,'_irf_calib*.pdf']);
150
    for j=1:length(a)
151 152 153
        delete([OutputDirectoryName,filesep,a(j).name]);
    end
    a=dir([OutputDirectoryName,filesep,fname_,'_',type,'_irf_restrictions.eps']);
154
    for j=1:length(a)
155 156 157
        delete([OutputDirectoryName,filesep,a(j).name]);
    end
    a=dir([OutputDirectoryName,filesep,fname_,'_',type,'_irf_restrictions.fig']);
158
    for j=1:length(a)
159 160 161
        delete([OutputDirectoryName,filesep,a(j).name]);
    end
    a=dir([OutputDirectoryName,filesep,fname_,'_',type,'_irf_restrictions.pdf']);
162
    for j=1:length(a)
163 164 165 166
        delete([OutputDirectoryName,filesep,a(j).name]);
    end
    disp('done !')
    skipline()
167

168
    % For single legend search which has maximum nbr of restrictions
169 170 171 172 173
    all_irf_couples = cellstr([char(endo_prior_restrictions.irf(:,1)) char(endo_prior_restrictions.irf(:,2))]);
    irf_couples = unique(all_irf_couples);
    nbr_irf_couples = size(irf_couples,1);
    plot_indx = NaN(nbr_irf_couples,1);
    time_matrix=cell(nbr_irf_couples,1);
174
    indx_irf_matrix=zeros(length(irestrictions),nbr_irf_couples);
175 176 177 178 179 180
    irf_matrix=cell(nbr_irf_couples,1);
    irf_mean=cell(nbr_irf_couples,1);
    irf_median=cell(nbr_irf_couples,1);
    irf_var=cell(nbr_irf_couples,1);
    irf_HPD=cell(nbr_irf_couples,1);
    irf_distrib=cell(nbr_irf_couples,1);
181 182
    maxijv=0;
    for ij=1:nbr_irf_restrictions
183 184
        if length(endo_prior_restrictions.irf{ij,3})>maxijv
            maxij=ij;maxijv=length(endo_prior_restrictions.irf{ij,3});
185
        end
186 187 188
        plot_indx(ij) = find(strcmp(irf_couples,all_irf_couples(ij,:)));
        time_matrix{plot_indx(ij)} = [time_matrix{plot_indx(ij)} endo_prior_restrictions.irf{ij,3}];
    end
189
    iplot_indx = ones(size(plot_indx));
190

191
    indx_irf = indx_irf(irestrictions,:);
192
    if ~DynareOptions.nograph
193
        h1=dyn_figure(DynareOptions.nodisplay,'name',[type ' evaluation of irf restrictions']);
194 195
        nrow=ceil(sqrt(nbr_irf_couples));
        ncol=nrow;
196
        if nrow*(nrow-1)>nbr_irf_couples
197 198
            ncol=nrow-1;
        end
199
    end
200
    for ij=1:nbr_irf_restrictions
201
        mat_irf{ij}=mat_irf{ij}(irestrictions,:);
202
        irf_matrix{plot_indx(ij)} = [irf_matrix{plot_indx(ij)} mat_irf{ij}];
203
        indx_irf_matrix(:,plot_indx(ij)) = indx_irf_matrix(:,plot_indx(ij)) + indx_irf(:,ij);
204
        for ik=1:size(mat_irf{ij},2)
205 206 207 208 209 210 211
            [Mean,Median,Var,HPD,Distrib] = ...
                posterior_moments(mat_irf{ij}(:,ik),0,DynareOptions.mh_conf_sig);
            irf_mean{plot_indx(ij)} = [irf_mean{plot_indx(ij)}; Mean];
            irf_median{plot_indx(ij)} = [irf_median{plot_indx(ij)}; Median];
            irf_var{plot_indx(ij)} = [irf_var{plot_indx(ij)}; Var];
            irf_HPD{plot_indx(ij)} = [irf_HPD{plot_indx(ij)}; HPD];
            irf_distrib{plot_indx(ij)} = [irf_distrib{plot_indx(ij)}; Distrib'];
212
        end
213 214
        leg = num2str(endo_prior_restrictions.irf{ij,3}(1));
        aleg = num2str(endo_prior_restrictions.irf{ij,3}(1));
215
        if size(mat_irf{ij},2)>1
216 217
            leg = [leg,':' ,num2str(endo_prior_restrictions.irf{ij,3}(end))];
            aleg = [aleg,'-' ,num2str(endo_prior_restrictions.irf{ij,3}(end))];
218 219
            iplot_indx(ij)=0;
        end
220
        if ~DynareOptions.nograph && length(time_matrix{plot_indx(ij)})==1
221
            set(0,'currentfigure',h1),
222 223 224
            subplot(nrow,ncol, plot_indx(ij)),
            hc = cumplot(mat_irf{ij}(:,ik));
            a=axis;
Marco Ratto's avatar
Marco Ratto committed
225
            delete(hc);
226 227 228
            x1val=max(endo_prior_restrictions.irf{ij,4}(1),a(1));
            x2val=min(endo_prior_restrictions.irf{ij,4}(2),a(2));
            hp = patch([x1val x2val x2val x1val],a([3 3 4 4]),'b');
Marco Ratto's avatar
Marco Ratto committed
229 230 231 232
            hold all,
            set(hp,'FaceColor', [0.7 0.8 1])
            hc = cumplot(mat_irf{ij}(:,ik));
            set(hc,'color','k','linewidth',2)
233 234 235 236 237 238 239 240 241 242 243
            hold off,
            %         hold off,
            title([endo_prior_restrictions.irf{ij,1},' vs ',endo_prior_restrictions.irf{ij,2}, '(', leg,')'],'interpreter','none'),
            %set(legend_h,'Xlim',[0 1]);
            %         if ij==maxij
            %             leg1 = num2str(endo_prior_restrictions.irf{ij,3}(:));
            %             [legend_h,object_h,plot_h,text_strings]=legend(leg1);
            %             Position=get(legend_h,'Position');Position(1:2)=[-0.055 0.95-Position(4)];
            %             set(legend_h,'Position',Position);
            %         end
        end
244 245 246
        % hc = get(h,'Children');
        %for i=2:2:length(hc)
        %end
247 248
        indx1 = find(indx_irf(:,ij)==0);
        indx2 = find(indx_irf(:,ij)~=0);
249
        atitle0=[endo_prior_restrictions.irf{ij,1},' vs ',endo_prior_restrictions.irf{ij,2}, '(', leg,')'];
250
        fprintf(['%4.1f%% of the ',type,' support matches IRF ',atitle0,' inside [%4.1f, %4.1f]\n'],length(indx1)/length(irestrictions)*100,endo_prior_restrictions.irf{ij,4})
251
        % aname=[type '_irf_calib_',int2str(ij)];
252 253 254 255
        aname=[type '_irf_calib_',endo_prior_restrictions.irf{ij,1},'_vs_',endo_prior_restrictions.irf{ij,2},'_',aleg];
        atitle=[type ' IRF Calib: Parameter(s) driving ',endo_prior_restrictions.irf{ij,1},' vs ',endo_prior_restrictions.irf{ij,2}, '(', leg,')'];
        options_mcf.amcf_name = aname;
        options_mcf.amcf_title = atitle;
256 257
        options_mcf.beha_title = 'IRF restriction';
        options_mcf.nobeha_title = 'NO IRF restriction';
258 259 260 261
        options_mcf.title = atitle0;
        if ~isempty(indx1) && ~isempty(indx2)
            mcf_analysis(xmat(:,nshock+1:end), indx1, indx2, options_mcf, DynareOptions);
        end
262

263 264 265 266 267
        %         [proba, dproba] = stab_map_1(xmat, indx1, indx2, aname, 0);
        %         indplot=find(proba<pvalue_ks);
        %         if ~isempty(indplot)
        %             stab_map_1(xmat, indx1, indx2, aname, 1, indplot, OutputDirectoryName,[],atitle);
        %         end
268
    end
269 270 271
    for ij=1:nbr_irf_couples
        if length(time_matrix{ij})>1
            if ~DynareOptions.nograph
272
                set(0,'currentfigure',h1);
273 274
                subplot(nrow,ncol, ij)
                itmp = (find(plot_indx==ij));
Marco Ratto's avatar
Marco Ratto committed
275
                htmp = plot(time_matrix{ij},[max(irf_matrix{ij})' min(irf_matrix{ij})'],'k--','linewidth',2);
276
                a=axis;
Marco Ratto's avatar
Marco Ratto committed
277
                delete(htmp);
278
                tmp=[];
279
                for ir=1:length(itmp)
280 281 282 283
                    for it=1:length(endo_prior_restrictions.irf{itmp(ir),3})
                        temp_index = find(time_matrix{ij}==endo_prior_restrictions.irf{itmp(ir),3}(it));
                        tmp(temp_index,:) = endo_prior_restrictions.irf{itmp(ir),4};
                    end
284
                end
285 286 287
                %             tmp = cell2mat(endo_prior_restrictions.irf(itmp,4));
                tmp(isinf(tmp(:,1)),1)=a(3);
                tmp(isinf(tmp(:,2)),2)=a(4);
Marco Ratto's avatar
Marco Ratto committed
288 289 290 291 292 293
                hp = patch([time_matrix{ij} time_matrix{ij}(end:-1:1)],[tmp(:,1); tmp(end:-1:1,2)],'c');
                set(hp,'FaceColor',[0.7 0.8 1])
                hold on,
                plot(time_matrix{ij},[max(irf_matrix{ij})' min(irf_matrix{ij})'],'k--','linewidth',2)
                plot(time_matrix{ij},irf_median{ij},'k','linewidth',2)
                plot(time_matrix{ij},[irf_distrib{ij}],'k-')
294
                plot(a(1:2),[0 0],'r')
295
                hold off
296
                axis([max(1,a(1)) a(2:4)])
297
                box on
298
                %set(gca,'xtick',sort(time_matrix{ij}))
299 300
                itmp = min(itmp);
                title([endo_prior_restrictions.irf{itmp,1},' vs ',endo_prior_restrictions.irf{itmp,2}],'interpreter','none'),
301
            end
302
            if any(iplot_indx.*plot_indx==ij)
303
                % MCF of the couples with logical AND
304
                itmp = min(find(plot_indx==ij));
305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323
                indx1 = find(indx_irf_matrix(:,ij)==0);
                indx2 = find(indx_irf_matrix(:,ij)~=0);
                leg = num2str(time_matrix{ij}(1));
                leg = [leg '...' num2str(time_matrix{ij}(end))];
                aleg = 'ALL';
                atitle0=[endo_prior_restrictions.irf{itmp,1},' vs ',endo_prior_restrictions.irf{itmp,2}, '(', leg,')'];
                fprintf(['%4.1f%% of the ',type,' support matches IRF restrictions ',atitle0,'\n'],length(indx1)/length(irestrictions)*100)
                % aname=[type '_irf_calib_',int2str(ij)];
                aname=[type '_irf_calib_',endo_prior_restrictions.irf{itmp,1},'_vs_',endo_prior_restrictions.irf{itmp,2},'_',aleg];
                atitle=[type ' IRF Calib: Parameter(s) driving ',endo_prior_restrictions.irf{itmp,1},' vs ',endo_prior_restrictions.irf{itmp,2}, '(', leg,')'];
                options_mcf.amcf_name = aname;
                options_mcf.amcf_title = atitle;
                options_mcf.beha_title = 'IRF restriction';
                options_mcf.nobeha_title = 'NO IRF restriction';
                options_mcf.title = atitle0;
                if ~isempty(indx1) && ~isempty(indx2)
                    mcf_analysis(xmat(:,nshock+1:end), indx1, indx2, options_mcf, DynareOptions);
                end
            end
324 325
        end
    end
326
    if ~DynareOptions.nograph
327
        dyn_saveas(h1,[OutputDirectoryName,filesep,fname_,'_',type,'_irf_restrictions'],DynareOptions.nodisplay,DynareOptions.graph_format);
328
        create_TeX_loader(DynareOptions,[OutputDirectoryName,filesep,fname_,'_',type,'_irf_restrictions'],[type ' evaluation of irf restrictions'],'irf_restrictions',type,DynareOptions.figures.textwidth*min(ij/ncol,1))
329
    end
330 331 332 333
    skipline()
end

if ~isempty(indx_moment)
334 335 336
    skipline()
    disp('Deleting old MOMENT calibration plots ...')
    a=dir([OutputDirectoryName,filesep,fname_,'_',type,'_moment_calib*.eps']);
337
    for j=1:length(a)
338 339 340
        delete([OutputDirectoryName,filesep,a(j).name]);
    end
    a=dir([OutputDirectoryName,filesep,fname_,'_',type,'_moment_calib*.fig']);
341
    for j=1:length(a)
342 343 344
        delete([OutputDirectoryName,filesep,a(j).name]);
    end
    a=dir([OutputDirectoryName,filesep,fname_,'_',type,'_moment_calib*.pdf']);
345
    for j=1:length(a)
346 347 348
        delete([OutputDirectoryName,filesep,a(j).name]);
    end
    a=dir([OutputDirectoryName,filesep,fname_,'_',type,'_moment_restrictions.eps']);
349
    for j=1:length(a)
350 351 352
        delete([OutputDirectoryName,filesep,a(j).name]);
    end
    a=dir([OutputDirectoryName,filesep,fname_,'_',type,'_moment_restrictions.fig']);
353
    for j=1:length(a)
354 355 356
        delete([OutputDirectoryName,filesep,a(j).name]);
    end
    a=dir([OutputDirectoryName,filesep,fname_,'_',type,'_moment_restrictions.pdf']);
357
    for j=1:length(a)
358 359 360 361
        delete([OutputDirectoryName,filesep,a(j).name]);
    end
    disp('done !')
    skipline()
362

363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380
    %get parameter names including standard deviations
    np=size(BayesInfo.name,1);
    name=cell(np,1);
    name_tex=cell(np,1);
    for jj=1:np
        if DynareOptions.TeX
            [param_name_temp, param_name_tex_temp]= get_the_name(jj,DynareOptions.TeX,Model,EstimatedParameters,DynareOptions);
            name_tex{jj,1} = strrep(param_name_tex_temp,'$','');
            name{jj,1} = param_name_temp;
        else
            param_name_temp = get_the_name(jj,DynareOptions.TeX,Model,EstimatedParameters,DynareOptions);
            name{jj,1} = param_name_temp;
        end
    end
    options_mcf.param_names = char(name);
    if DynareOptions.TeX
        options_mcf.param_names_tex = char(name_tex);
    end
381 382 383 384 385 386
    options_mcf.param_names = char(BayesInfo.name);
    all_moment_couples = cellstr([char(endo_prior_restrictions.moment(:,1)) char(endo_prior_restrictions.moment(:,2))]);
    moment_couples = unique(all_moment_couples);
    nbr_moment_couples = size(moment_couples,1);
    plot_indx = NaN(nbr_moment_couples,1);
    time_matrix=cell(nbr_moment_couples,1);
387
    indx_moment_matrix=zeros(length(irestrictions),nbr_moment_couples);
388 389 390 391 392 393
    moment_matrix=cell(nbr_moment_couples,1);
    moment_mean=cell(nbr_moment_couples,1);
    moment_median=cell(nbr_moment_couples,1);
    moment_var=cell(nbr_moment_couples,1);
    moment_HPD=cell(nbr_moment_couples,1);
    moment_distrib=cell(nbr_moment_couples,1);
394
    % For single legend search which has maximum nbr of restrictions
395 396 397 398 399 400 401 402
    maxijv=0;
    for ij=1:nbr_moment_restrictions
        if length(endo_prior_restrictions.moment{ij,3})>maxijv
            maxij=ij;maxijv=length(endo_prior_restrictions.moment{ij,3});
        end
        plot_indx(ij) = find(strcmp(moment_couples,all_moment_couples(ij,:)));
        time_matrix{plot_indx(ij)} = [time_matrix{plot_indx(ij)} endo_prior_restrictions.moment{ij,3}];
    end
403
    iplot_indx = ones(size(plot_indx));
404

405
    indx_moment = indx_moment(irestrictions,:);
406
    if ~DynareOptions.nograph
407
        h2=dyn_figure(DynareOptions.nodisplay,'name',[type ' evaluation of moment restrictions']);
408 409
        nrow=ceil(sqrt(nbr_moment_couples));
        ncol=nrow;
410
        if nrow*(nrow-1)>nbr_moment_couples
411 412
            ncol=nrow-1;
        end
413
    end
414

415
    for ij=1:nbr_moment_restrictions
416
        mat_moment{ij}=mat_moment{ij}(irestrictions,:);
417
        moment_matrix{plot_indx(ij)} = [moment_matrix{plot_indx(ij)} mat_moment{ij}];
418
        indx_moment_matrix(:,plot_indx(ij)) = indx_moment_matrix(:,plot_indx(ij)) + indx_moment(:,ij);
419
        for ik=1:size(mat_moment{ij},2)
420 421 422 423 424 425 426
            [Mean,Median,Var,HPD,Distrib] = ...
                posterior_moments(mat_moment{ij}(:,ik),0,DynareOptions.mh_conf_sig);
            moment_mean{plot_indx(ij)} = [moment_mean{plot_indx(ij)}; Mean];
            moment_median{plot_indx(ij)} = [moment_median{plot_indx(ij)}; Median];
            moment_var{plot_indx(ij)} = [moment_var{plot_indx(ij)}; Var];
            moment_HPD{plot_indx(ij)} = [moment_HPD{plot_indx(ij)}; HPD];
            moment_distrib{plot_indx(ij)} = [moment_distrib{plot_indx(ij)}; Distrib'];
427
        end
428 429
        leg = num2str(endo_prior_restrictions.moment{ij,3}(1));
        aleg = num2str(endo_prior_restrictions.moment{ij,3}(1));
430
        if size(mat_moment{ij},2)>1
431 432
            leg = [leg,':' ,num2str(endo_prior_restrictions.moment{ij,3}(end))];
            aleg = [aleg,'_' ,num2str(endo_prior_restrictions.moment{ij,3}(end))];
433
            iplot_indx(ij)=0;
434
        end
435
        if ~DynareOptions.nograph && length(time_matrix{plot_indx(ij)})==1
436
            set(0,'currentfigure',h2);
437 438
            subplot(nrow,ncol,plot_indx(ij)),
            hc = cumplot(mat_moment{ij}(:,ik));
Marco Ratto's avatar
Marco Ratto committed
439
            a=axis; delete(hc),
440 441 442 443
            %     hist(mat_moment{ij}),
            x1val=max(endo_prior_restrictions.moment{ij,4}(1),a(1));
            x2val=min(endo_prior_restrictions.moment{ij,4}(2),a(2));
            hp = patch([x1val x2val x2val x1val],a([3 3 4 4]),'b');
Marco Ratto's avatar
Marco Ratto committed
444
            set(hp,'FaceColor', [0.7 0.8 1])
445
            hold all
Marco Ratto's avatar
Marco Ratto committed
446 447
            hc = cumplot(mat_moment{ij}(:,ik));
            set(hc,'color','k','linewidth',2)
448
            hold off
449
            title([endo_prior_restrictions.moment{ij,1},' vs ',endo_prior_restrictions.moment{ij,2},'(',leg,')'],'interpreter','none'),
450 451 452 453 454 455
            %         if ij==maxij
            %             leg1 = num2str(endo_prior_restrictions.moment{ij,3}(:));
            %             [legend_h,object_h,plot_h,text_strings]=legend(leg1);
            %             Position=get(legend_h,'Position');Position(1:2)=[-0.055 0.95-Position(4)];
            %             set(legend_h,'Position',Position);
            %         end
456 457 458
        end
        indx1 = find(indx_moment(:,ij)==0);
        indx2 = find(indx_moment(:,ij)~=0);
459
        atitle0=[endo_prior_restrictions.moment{ij,1},' vs ',endo_prior_restrictions.moment{ij,2}, '(', leg,')'];
460
        fprintf(['%4.1f%% of the ',type,' support matches MOMENT ',atitle0,' inside [%4.1f, %4.1f]\n'],length(indx1)/length(irestrictions)*100,endo_prior_restrictions.moment{ij,4})
461
        % aname=[type '_moment_calib_',int2str(ij)];
462 463 464 465
        aname=[type '_moment_calib_',endo_prior_restrictions.moment{ij,1},'_vs_',endo_prior_restrictions.moment{ij,2},'_',aleg];
        atitle=[type ' MOMENT Calib: Parameter(s) driving ',endo_prior_restrictions.moment{ij,1},' vs ',endo_prior_restrictions.moment{ij,2}, '(', leg,')'];
        options_mcf.amcf_name = aname;
        options_mcf.amcf_title = atitle;
466 467
        options_mcf.beha_title = 'moment restriction';
        options_mcf.nobeha_title = 'NO moment restriction';
468 469 470 471
        options_mcf.title = atitle0;
        if ~isempty(indx1) && ~isempty(indx2)
            mcf_analysis(xmat, indx1, indx2, options_mcf, DynareOptions);
        end
472

473 474 475 476 477
        %         [proba, dproba] = stab_map_1(xmat, indx1, indx2, aname, 0);
        %         indplot=find(proba<pvalue_ks);
        %         if ~isempty(indplot)
        %             stab_map_1(xmat, indx1, indx2, aname, 1, indplot, OutputDirectoryName,[],atitle);
        %         end
478
    end
479 480
    for ij=1:nbr_moment_couples
        if length(time_matrix{ij})>1
481 482 483 484
            if ~DynareOptions.nograph
                itmp = (find(plot_indx==ij));
                set(0,'currentfigure',h2);
                subplot(nrow,ncol, ij)
Marco Ratto's avatar
Marco Ratto committed
485
                htmp = plot(time_matrix{ij},[max(moment_matrix{ij})' min(moment_matrix{ij})'],'k--','linewidth',2);
486
                a=axis;
Marco Ratto's avatar
Marco Ratto committed
487
                delete(htmp);
488
                tmp=[];
489
                for ir=1:length(itmp)
490 491 492 493
                    for it=1:length(endo_prior_restrictions.moment{itmp(ir),3})
                        temp_index = find(time_matrix{ij}==endo_prior_restrictions.moment{itmp(ir),3}(it));
                        tmp(temp_index,:) = endo_prior_restrictions.moment{itmp(ir),4};
                    end
494
                end
495 496 497 498
                %             tmp = cell2mat(endo_prior_restrictions.moment(itmp,4));
                tmp(isinf(tmp(:,1)),1)=a(3);
                tmp(isinf(tmp(:,2)),2)=a(4);
                hp = patch([time_matrix{ij} time_matrix{ij}(end:-1:1)],[tmp(:,1); tmp(end:-1:1,2)],'b');
Marco Ratto's avatar
Marco Ratto committed
499
                set(hp,'FaceColor',[0.7 0.8 1])
500
                hold on
Marco Ratto's avatar
Marco Ratto committed
501 502 503
                plot(time_matrix{ij},[max(moment_matrix{ij})' min(moment_matrix{ij})'],'k--','linewidth',2)
                plot(time_matrix{ij},moment_median{ij},'k','linewidth',2)
                plot(time_matrix{ij},[moment_distrib{ij}],'k-')
504
                plot(a(1:2),[0 0],'r')
505
                hold off
506
                axis(a)
507
                box on
508 509 510
                set(gca,'xtick',sort(time_matrix{ij}))
                itmp = min(itmp);
                title([endo_prior_restrictions.moment{itmp,1},' vs ',endo_prior_restrictions.moment{itmp,2}],'interpreter','none'),
511
            end
512
            if any(iplot_indx.*plot_indx==ij)
513
                % MCF of the couples with logical AND
514
                itmp = min(find(plot_indx==ij));
515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533
                indx1 = find(indx_moment_matrix(:,ij)==0);
                indx2 = find(indx_moment_matrix(:,ij)~=0);
                leg = num2str(time_matrix{ij}(1));
                leg = [leg '...' num2str(time_matrix{ij}(end))];
                aleg = 'ALL';
                atitle0=[endo_prior_restrictions.moment{itmp,1},' vs ',endo_prior_restrictions.moment{itmp,2}, '(', leg,')'];
                fprintf(['%4.1f%% of the ',type,' support matches MOMENT restrictions ',atitle0,'\n'],length(indx1)/length(irestrictions)*100)
                % aname=[type '_moment_calib_',int2str(ij)];
                aname=[type '_moment_calib_',endo_prior_restrictions.moment{itmp,1},'_vs_',endo_prior_restrictions.moment{itmp,2},'_',aleg];
                atitle=[type ' MOMENT Calib: Parameter(s) driving ',endo_prior_restrictions.moment{itmp,1},' vs ',endo_prior_restrictions.moment{itmp,2}, '(', leg,')'];
                options_mcf.amcf_name = aname;
                options_mcf.amcf_title = atitle;
                options_mcf.beha_title = 'moment restriction';
                options_mcf.nobeha_title = 'NO moment restriction';
                options_mcf.title = atitle0;
                if ~isempty(indx1) && ~isempty(indx2)
                    mcf_analysis(xmat, indx1, indx2, options_mcf, DynareOptions);
                end
            end
534 535
        end
    end
536
    if ~DynareOptions.nograph
537
        dyn_saveas(h2,[OutputDirectoryName,filesep,fname_,'_',type,'_moment_restrictions'],DynareOptions.nodisplay,DynareOptions.graph_format);
538
        create_TeX_loader(DynareOptions,[OutputDirectoryName,filesep,fname_,'_',type,'_moment_restrictions'],[type ' evaluation of moment restrictions'],'moment_restrictions',type,DynareOptions.figures.textwidth*min(ij/ncol,1))
539
    end
540

541 542 543 544
    skipline()
end
return

545
function []=create_TeX_loader(options,figpath,caption,label_name,label_type,scale_factor)
546
if options.TeX && any(strcmp('eps',cellstr(options.graph_format)))
547
    fidTeX = fopen([figpath '.tex'],'w');
548 549 550 551
    fprintf(fidTeX,'%% TeX eps-loader file generated by map_calibration.m (Dynare).\n');
    fprintf(fidTeX,['%% ' datestr(now,0) '\n\n']);
    fprintf(fidTeX,'\\begin{figure}[H]\n');
    fprintf(fidTeX,'\\centering \n');
552
    fprintf(fidTeX,'\\includegraphics[width=%2.2f\\textwidth]{%s}\n',scale_factor,strrep(figpath,'\','/'));
553 554 555 556 557 558
    fprintf(fidTeX,'\\caption{%s.}',caption);
    fprintf(fidTeX,'\\label{Fig:%s:%s}\n',label_name,label_type);
    fprintf(fidTeX,'\\end{figure}\n\n');
    fprintf(fidTeX,'%% End Of TeX file. \n');
    fclose(fidTeX);
end