stab_map_.m 28.7 KB
Newer Older
1
function x0 = stab_map_(OutputDirectoryName,opt_gsa)
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
%
% function x0 = stab_map_(OutputDirectoryName)
%
% Mapping of stability regions in the prior ranges applying
% Monte Carlo filtering techniques.
%
% INPUTS (from opt_gsa structure)
% Nsam = MC sample size
% fload = 0 to run new MC; 1 to load prevoiusly generated analysis
% alpha2 =  significance level for bivariate sensitivity analysis
% [abs(corrcoef) > alpha2]
% prepSA = 1: save transition matrices for mapping reduced form
%        = 0: no transition matrix saved (default)
% pprior = 1: sample from prior ranges (default): sample saved in
%            _prior.mat   file
%        = 0: sample from posterior ranges: sample saved in
%            _mc.mat file
% OUTPUT:
% x0: one parameter vector for which the model is stable.
%
% GRAPHS
% 1) Pdf's of marginal distributions under the stability (dotted
%     lines) and unstability (solid lines) regions
% 2) Cumulative distributions of:
%   - stable subset (dotted lines)
%   - unacceptable subset (solid lines)
% 3) Bivariate plots of significant correlation patterns
%  ( abs(corrcoef) > alpha2) under the stable and unacceptable subsets
%
31
% USES qmc_sequence, stab_map_1, stab_map_2
32
%
33
% Written by Marco Ratto
34
% Joint Research Centre, The European Commission,
35
% marco.ratto@ec.europa.eu
36

Marco Ratto's avatar
Marco Ratto committed
37
% Copyright (C) 2012-2016 European Commission
Stéphane Adjemian's avatar
Stéphane Adjemian committed
38
% Copyright (C) 2012-2017 Dynare Team
39
40
41
42
43
44
45
46
47
48
49
50
%
% 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.
51
%
52
53
% You should have received a copy of the GNU General Public License
% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
54
55
56
57

%global bayestopt_ estim_params_ dr_ options_ ys_ fname_
global bayestopt_ estim_params_ options_ oo_ M_

58
% opt_gsa=options_.opt_gsa;
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74

Nsam   = opt_gsa.Nsam;
fload  = opt_gsa.load_stab;
alpha2 = opt_gsa.alpha2_stab;
pvalue_ks = opt_gsa.pvalue_ks;
pvalue_corr = opt_gsa.pvalue_corr;
prepSA = (opt_gsa.redform | opt_gsa.identification);
pprior = opt_gsa.pprior;
neighborhood_width = opt_gsa.neighborhood_width;
ilptau = opt_gsa.ilptau;
nliv   = opt_gsa.morris_nliv;
ntra   = opt_gsa.morris_ntra;

dr_ = oo_.dr;
%if isfield(dr_,'ghx'),
ys_ = oo_.dr.ys;
75
76
77
nspred = M_.nspred; %size(dr_.ghx,2);
nboth = M_.nboth;
nfwrd = M_.nfwrd;
78
79
80
81
82
83
84
85
%end
fname_ = M_.fname;

np = estim_params_.np;
nshock = estim_params_.nvx;
nshock = nshock + estim_params_.nvn;
nshock = nshock + estim_params_.ncx;
nshock = nshock + estim_params_.ncn;
86
lpmat0=zeros(Nsam,0);
87
xparam1=[];
88
89
90
91
92
93
94

pshape = bayestopt_.pshape(nshock+1:end);
p1 = bayestopt_.p1(nshock+1:end);
p2 = bayestopt_.p2(nshock+1:end);
p3 = bayestopt_.p3(nshock+1:end);
p4 = bayestopt_.p4(nshock+1:end);

95
96
97
[junk1,junk2,junk3,lb,ub,junk4] = set_prior(estim_params_,M_,options_); %Prepare bounds
if ~isempty(bayestopt_) && any(bayestopt_.pshape > 0)
    % Set prior bounds
98
    bounds = prior_bounds(bayestopt_, options_.prior_trunc);
99
100
101
    bounds.lb = max(bounds.lb,lb);
    bounds.ub = min(bounds.ub,ub);
else  % estimated parameters but no declared priors
102
103
      % No priors are declared so Dynare will estimate the model by
      % maximum likelihood with inequality constraints for the parameters.
104
105
    bounds.lb = lb;
    bounds.ub = ub;
106
107
108
109
    if opt_gsa.prior_range==0
        warning('GSA:: When using ML, sampling from the prior is not possible. Setting prior_range=1')
        opt_gsa.prior_range=1;
    end
110
end
111

112
if nargin==0
113
114
115
    OutputDirectoryName='';
end

116
117
118
options_mcf.pvalue_ks = pvalue_ks;
options_mcf.pvalue_corr = pvalue_corr;
options_mcf.alpha2 = alpha2;
119
120
121

name=cell(np,1);
name_tex=cell(np,1);
122
for jj=1:np
123
124
125
126
127
128
129
130
131
132
133
134
135
136
    if options_.TeX
        [param_name_temp, param_name_tex_temp]= get_the_name(nshock+jj,options_.TeX,M_,estim_params_,options_);
        name_tex{jj,1} = strrep(param_name_tex_temp,'$','');
        name{jj,1} = param_name_temp;
    else
        param_name_temp = get_the_name(nshock+jj,options_.TeX,M_,estim_params_,options_);
        name{jj,1} = param_name_temp;
    end
end
if options_.TeX
    options_mcf.param_names_tex=char(name_tex);
end
options_mcf.param_names = char(name);

137
138
options_mcf.fname_ = fname_;
options_mcf.OutputDirectoryName = OutputDirectoryName;
139
options_mcf.xparam1 = [];
140

141
142
143
144
145
opt=options_;
options_.periods=0;
options_.nomoments=1;
options_.irf=0;
options_.noprint=1;
146
if fload==0
147
148
149
150
    %   if prepSA
    %     T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),Nsam/2);
    %   end

151
    if isfield(dr_,'ghx')
152
153
154
155
156
157
158
159
160
161
        egg=zeros(length(dr_.eigval),Nsam);
    end
    yys=zeros(length(dr_.ys),Nsam);

    if opt_gsa.morris == 1
        [lpmat, OutFact] = Sampling_Function_2(nliv, np+nshock, ntra, ones(np+nshock, 1), zeros(np+nshock,1), []);
        lpmat = lpmat.*(nliv-1)/nliv+1/nliv/2;
        Nsam=size(lpmat,1);
        lpmat0 = lpmat(:,1:nshock);
        lpmat = lpmat(:,nshock+1:end);
162
163
164
        %     elseif opt_gsa.morris==3,
        %         lpmat = prep_ide(Nsam,np,5);
        %         Nsam=size(lpmat,1);
165
    else
166
        if np<52 && ilptau>0
167
            [lpmat] = qmc_sequence(np, int64(1), 0, Nsam)';
168
169
            if np>30 || ilptau==2 % scrambled lptau
                for j=1:np
170
171
172
173
                    lpmat(:,j)=lpmat(randperm(Nsam),j);
                end
            end
        else %ilptau==0
174
            [lpmat] = NaN(Nsam,np);
175
            for j=1:np
176
177
178
179
180
181
                lpmat(:,j) = randperm(Nsam)'./(Nsam+1); %latin hypercube
            end

        end
    end
    %   try
182
    dummy=prior_draw_gsa(1); %initialize persistent variables
183
184
185
186
187
188
189
190
                             %   catch
                             %     if pprior,
                             %       if opt_gsa.prior_range==0;
                             %         error('Some unknown prior is specified or ML estimation,: use prior_range=1 option!!');
                             %       end
                             %     end
                             %
                             %   end
191
192
193
    if pprior
        for j=1:nshock
            if opt_gsa.morris~=1
194
195
196
                lpmat0(:,j) = randperm(Nsam)'./(Nsam+1); %latin hypercube
            end
            if opt_gsa.prior_range
197
                lpmat0(:,j)=lpmat0(:,j).*(bounds.ub(j)-bounds.lb(j))+bounds.lb(j);
198
199
200
201
202
203
204
205
206
            end
        end
        if opt_gsa.prior_range
            %       if opt_gsa.identification,
            %         deltx=min(0.001, 1/Nsam/2);
            %         for j=1:np,
            %           xdelt(:,:,j)=prior_draw_gsa(0,[lpmat0 lpmat]+deltx);
            %         end
            %       end
207
            for j=1:np
208
                lpmat(:,j)=lpmat(:,j).*(bounds.ub(j+nshock)-bounds.lb(j+nshock))+bounds.lb(j+nshock);
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
            end
        else
            xx=prior_draw_gsa(0,[lpmat0 lpmat]);
            %       if opt_gsa.identification,
            %         deltx=min(0.001, 1/Nsam/2);
            %         ldum=[lpmat0 lpmat];
            %         ldum = prior_draw_gsa(0,ldum+deltx);
            %         for j=1:nshock+np,
            %           xdelt(:,:,j)=xx;
            %           xdelt(:,j,j)=ldum(:,j);
            %         end
            %         clear ldum
            %       end
            lpmat0=xx(:,1:nshock);
            lpmat=xx(:,nshock+1:end);
            clear xx;
        end
    else
        %         for j=1:nshock,
        %             xparam1(j) = oo_.posterior_mode.shocks_std.(bayestopt_.name{j});
        %             sd(j) = oo_.posterior_std.shocks_std.(bayestopt_.name{j});
        %             lpmat0(:,j) = randperm(Nsam)'./(Nsam+1); %latin hypercube
        %             lb = max(bayestopt_.lb(j), xparam1(j)-2*sd(j));
        %             ub1=xparam1(j)+(xparam1(j) - lb); % define symmetric range around the mode!
        %             ub = min(bayestopt_.ub(j),ub1);
        %             if ub<ub1,
        %                 lb=xparam1(j)-(ub-xparam1(j)); % define symmetric range around the mode!
        %             end
        %             lpmat0(:,j) = lpmat0(:,j).*(ub-lb)+lb;
        %         end
        %         %
        %         for j=1:np,
        %             xparam1(j+nshock) = oo_.posterior_mode.parameters.(bayestopt_.name{j+nshock});
        %             sd(j+nshock) = oo_.posterior_std.parameters.(bayestopt_.name{j+nshock});
        %             lb = max(bayestopt_.lb(j+nshock),xparam1(j+nshock)-2*sd(j+nshock));
        %             ub1=xparam1(j+nshock)+(xparam1(j+nshock) - lb); % define symmetric range around the mode!
        %             ub = min(bayestopt_.ub(j+nshock),ub1);
        %             if ub<ub1,
        %                 lb=xparam1(j+nshock)-(ub-xparam1(j+nshock)); % define symmetric range around the mode!
        %             end
        %             %ub = min(bayestopt_.ub(j+nshock),xparam1(j+nshock)+2*sd(j+nshock));
        %             if np>30 & np<52
        %                 lpmat(:,j) = lpmat(randperm(Nsam),j).*(ub-lb)+lb;
        %             else
        %                 lpmat(:,j) = lpmat(:,j).*(ub-lb)+lb;
        %             end
        %         end
        %load([fname_,'_mode'])
257
        if neighborhood_width>0 && isempty(options_.mode_file)
258
259
            xparam1 = get_all_parameters(estim_params_,M_);
        else
260
            eval(['load ' options_.mode_file '.mat;']);
261
        end
262
263
        if neighborhood_width>0
            for j=1:nshock
264
                if opt_gsa.morris ~= 1
265
                    lpmat0(:,j) = randperm(Nsam)'./(Nsam+1); %latin hypercube
266
                end
267
268
                ub=min([bounds.ub(j) xparam1(j)*(1+neighborhood_width)]);
                lb=max([bounds.lb(j) xparam1(j)*(1-neighborhood_width)]);
269
270
                lpmat0(:,j)=lpmat0(:,j).*(ub-lb)+lb;
            end
271
            for j=1:np
272
273
                ub=xparam1(j+nshock)*(1+sign(xparam1(j+nshock))*neighborhood_width);
                lb=xparam1(j+nshock)*(1-sign(xparam1(j+nshock))*neighborhood_width);
274
                if bounds.ub(j+nshock)>=xparam1(j) && bounds.lb(j)<=xparam1(j+nshock)
275
276
277
278
279
                    ub=min([bounds.ub(j+nshock) ub]);
                    lb=max([bounds.lb(j+nshock) lb]);
                else
                    fprintf('\nstab_map_:: the calibrated value of param %s for neighborhood_width sampling is outside prior bounds.\nWe allow violation of bounds for this parameter, but if this was not done on purpose, please change calibration before running neighborhood_width sampling\n', bayestopt_.name{j+nshock})
                end
280
281
282
283
284
                lpmat(:,j)=lpmat(:,j).*(ub-lb)+lb;
            end
        else
            d = chol(inv(hh));
            lp=randn(Nsam*2,nshock+np)*d+kron(ones(Nsam*2,1),xparam1');
285
            lnprior=zeros(1,Nsam*2);
286
            for j=1:Nsam*2
287
                lnprior(j) = any(lp(j,:)'<=bounds.lb | lp(j,:)'>=bounds.ub);
288
289
290
291
292
293
294
295
296
297
298
            end
            ireal=[1:2*Nsam];
            ireal=ireal(find(lnprior==0));
            lp=lp(ireal,:);
            Nsam=min(Nsam, length(ireal));
            lpmat0=lp(1:Nsam,1:nshock);
            lpmat=lp(1:Nsam,nshock+1:end);
            clear lp lnprior ireal;
        end
    end
    %
299
    h = dyn_waitbar(0,'Please wait...');
300
301
302
303
304
    istable=[1:Nsam];
    jstab=0;
    iunstable=[1:Nsam];
    iindeterm=zeros(1,Nsam);
    iwrong=zeros(1,Nsam);
305
306
    inorestriction=zeros(1,Nsam);
    irestriction=zeros(1,Nsam);
307
    infox=zeros(Nsam,1);
308
    for j=1:Nsam
309
        M_ = set_all_parameters([lpmat0(j,:) lpmat(j,:)]',estim_params_,M_);
310
311
        %try stoch_simul([]);
        try
312
313
314
315
316
            if ~ isempty(options_.endogenous_prior_restrictions.moment)
                [Tt,Rr,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_);
            else
                [Tt,Rr,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_,'restrict');
            end
317
            infox(j,1)=info(1);
318
            if infox(j,1)==0 && ~exist('T','var')
319
                dr_=oo_.dr;
320
                if prepSA
321
322
                    try
                        T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),Nsam);
323
324
                    catch
                        ME = lasterror();
325
                        if strcmp('MATLAB:nomem',ME.identifier)
326
327
328
329
330
331
332
333
334
                            prepSA=0;
                            disp('The model is too large for storing state space matrices ...')
                            disp('for mapping reduced form or for identification')
                        end
                        T=[];
                    end
                else
                    T=[];
                end
335
336
                egg=zeros(length(dr_.eigval),Nsam);
            end
337
            if infox(j,1)
338
                %                 disp('no solution'),
339
                if isfield(oo_.dr,'ghx')
340
341
                    oo_.dr=rmfield(oo_.dr,'ghx');
                end
342
                if (infox(j,1)<3 || infox(j,1)>5) && isfield(oo_.dr,'eigval')
343
344
                    oo_.dr=rmfield(oo_.dr,'eigval');
                end
345
            end
346
        catch ME
347
            if isfield(oo_.dr,'eigval')
348
349
                oo_.dr=rmfield(oo_.dr,'eigval');
            end
350
            if isfield(oo_.dr,'ghx')
351
352
                oo_.dr=rmfield(oo_.dr,'ghx');
            end
353
            disp('No solution could be found')
354
355
        end
        dr_ = oo_.dr;
356
        if isfield(dr_,'ghx')
357
358
359
360
361
362
363
364
365
            egg(:,j) = sort(dr_.eigval);
            if prepSA
                jstab=jstab+1;
                T(:,:,jstab) = [dr_.ghx dr_.ghu];
                %         [A,B] = ghx2transition(squeeze(T(:,:,jstab)), ...
                %           bayestopt_.restrict_var_list, ...
                %           bayestopt_.restrict_columns, ...
                %           bayestopt_.restrict_aux);
            end
366
            if ~exist('nspred','var')
367
368
369
370
                nspred = dr_.nspred; %size(dr_.ghx,2);
                nboth = dr_.nboth;
                nfwrd = dr_.nfwrd;
            end
371
372
            info=endogenous_prior_restrictions(Tt,Rr,M_,options_,oo_);
            infox(j,1)=info(1);
373
            if info(1)
374
375
376
377
378
                inorestriction(j)=j;
            else
                iunstable(j)=0;
                irestriction(j)=j;
            end
379
380
381
382
        else
            istable(j)=0;
            if isfield(dr_,'eigval')
                egg(:,j) = sort(dr_.eigval);
383
                if exist('nspred','var')
384
385
386
                    if any(isnan(egg(1:nspred,j)))
                        iwrong(j)=j;
                    else
387
                        if (nboth || nfwrd) && abs(egg(nspred+1,j))<=options_.qz_criterium
388
389
390
391
392
                            iindeterm(j)=j;
                        end
                    end
                end
            else
393
                if exist('egg','var')
394
395
396
397
398
399
400
401
                    egg(:,j)=ones(size(egg,1),1).*NaN;
                end
                iwrong(j)=j;
            end
        end
        ys_=real(dr_.ys);
        yys(:,j) = ys_;
        ys_=yys(:,1);
402
        dyn_waitbar(j/Nsam,h,['MC iteration ',int2str(j),'/',int2str(Nsam)])
403
    end
404
    dyn_waitbar_close(h);
405
    if prepSA && jstab
406
        T=T(:,:,1:jstab);
407
408
    else
        T=[];
409
    end
410
411
412
413
    istable=istable(find(istable));  % stable params ignoring restrictions
    irestriction=irestriction(find(irestriction));  % stable params & restrictions OK
    inorestriction=inorestriction(find(inorestriction));  % stable params violating restrictions
    iunstable=iunstable(find(iunstable));   % violation of BK & restrictions & solution could not be found (whatever goes wrong)
414
415
    iindeterm=iindeterm(find(iindeterm));  % indeterminacy
    iwrong=iwrong(find(iwrong));  % dynare could not find solution
416
    ixun=iunstable(find(~ismember(iunstable,[iindeterm,iwrong,inorestriction]))); % explosive roots
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456

    %     % map stable samples
    %     istable=[1:Nsam];
    %     for j=1:Nsam,
    %         if any(isnan(egg(1:nspred,j)))
    %             istable(j)=0;
    %         else
    %             if abs(egg(nspred,j))>=options_.qz_criterium; %(1-(options_.qz_criterium-1)); %1-1.e-5;
    %                 istable(j)=0;
    %                 %elseif (dr_.nboth | dr_.nfwrd) & abs(egg(nspred+1,j))<=options_.qz_criterium; %1+1.e-5;
    %             elseif (nboth | nfwrd) & abs(egg(nspred+1,j))<=options_.qz_criterium; %1+1.e-5;
    %                 istable(j)=0;
    %             end
    %         end
    %     end
    %     istable=istable(find(istable));  % stable params
    %
    %     % map unstable samples
    %     iunstable=[1:Nsam];
    %     for j=1:Nsam,
    %         %if abs(egg(dr_.npred+1,j))>1+1.e-5 & abs(egg(dr_.npred,j))<1-1.e-5;
    %         %if (dr_.nboth | dr_.nfwrd),
    %         if ~any(isnan(egg(1:5,j)))
    %             if (nboth | nfwrd),
    %                 if abs(egg(nspred+1,j))>options_.qz_criterium & abs(egg(nspred,j))<options_.qz_criterium; %(1-(options_.qz_criterium-1));
    %                     iunstable(j)=0;
    %                 end
    %             else
    %                 if abs(egg(nspred,j))<options_.qz_criterium; %(1-(options_.qz_criterium-1));
    %                     iunstable(j)=0;
    %                 end
    %             end
    %         end
    %     end
    %     iunstable=iunstable(find(iunstable));   % unstable params
    bkpprior.pshape=bayestopt_.pshape;
    bkpprior.p1=bayestopt_.p1;
    bkpprior.p2=bayestopt_.p2;
    bkpprior.p3=bayestopt_.p3;
    bkpprior.p4=bayestopt_.p4;
457
    if pprior
458
        if ~prepSA
459
            save([OutputDirectoryName filesep fname_ '_prior.mat'], ...
460
461
                 'bkpprior','lpmat','lpmat0','irestriction','iunstable','istable','iindeterm','iwrong','ixun', ...
                 'egg','yys','nspred','nboth','nfwrd','infox')
462
        else
463
            save([OutputDirectoryName filesep fname_ '_prior.mat'], ...
464
465
                 'bkpprior','lpmat','lpmat0','irestriction','iunstable','istable','iindeterm','iwrong','ixun', ...
                 'egg','yys','T','nspred','nboth','nfwrd','infox')
466
467
468
469
        end

    else
        if ~prepSA
470
            save([OutputDirectoryName filesep fname_ '_mc.mat'], ...
471
472
                 'lpmat','lpmat0','irestriction','iunstable','istable','iindeterm','iwrong','ixun', ...
                 'egg','yys','nspred','nboth','nfwrd','infox')
473
        else
474
            save([OutputDirectoryName filesep fname_ '_mc.mat'], ...
475
476
                 'lpmat','lpmat0','irestriction','iunstable','istable','iindeterm','iwrong','ixun', ...
                 'egg','yys','T','nspred','nboth','nfwrd','infox')
477
478
479
        end
    end
else
480
    if pprior
481
        filetoload=[OutputDirectoryName filesep fname_ '_prior.mat'];
482
    else
483
        filetoload=[OutputDirectoryName filesep fname_ '_mc.mat'];
484
    end
485
    load(filetoload,'lpmat','lpmat0','irestriction','iunstable','istable','iindeterm','iwrong','ixun','egg','yys','nspred','nboth','nfwrd','infox')
486
    Nsam = size(lpmat,1);
487
    if pprior==0 && ~isempty(options_.mode_file)
488
489
        eval(['load ' options_.mode_file '.mat;']);
    end
490
491


492
    if prepSA && isempty(strmatch('T',who('-file', filetoload),'exact'))
493
        h = dyn_waitbar(0,'Please wait...');
494
495
496
497
498
499
500
        options_.periods=0;
        options_.nomoments=1;
        options_.irf=0;
        options_.noprint=1;
        stoch_simul([]);
        %T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),length(istable));
        ntrans=length(istable);
501
        yys=NaN(length(ys_),ntrans);
502
        for j=1:ntrans
503
504
505
506
507
508
509
            M_.params(estim_params_.param_vals(:,1)) = lpmat(istable(j),:)';
            %stoch_simul([]);
            [Tt,Rr,SteadyState,info,M_,options_,oo_] = dynare_resolve(M_,options_,oo_,'restrict');
            % This syntax is not compatible with the current version of dynare_resolve [stepan].
            %[Tt,Rr,SteadyState,info] = dynare_resolve(bayestopt_.restrict_var_list,...
            %    bayestopt_.restrict_columns,...
            %    bayestopt_.restrict_aux);
510
            if ~exist('T','var')
511
512
513
514
                T=zeros(size(dr_.ghx,1),size(dr_.ghx,2)+size(dr_.ghu,2),ntrans);
            end
            dr_ = oo_.dr;
            T(:,:,j) = [dr_.ghx dr_.ghu];
515
            if ~exist('nspred','var')
516
517
518
519
520
521
522
                nspred = dr_.nspred; %size(dr_.ghx,2);
                nboth = dr_.nboth;
                nfwrd = dr_.nfwrd;
            end
            ys_=real(dr_.ys);
            yys(:,j) = ys_;
            ys_=yys(:,1);
523
            dyn_waitbar(j/ntrans,h,['MC iteration ',int2str(j),'/',int2str(ntrans)])
524
        end
525
        dyn_waitbar_close(h);
526
527
528
529
530
531
532
        save(filetoload,'T','-append')
    elseif prepSA
        load(filetoload,'T')
    end
end

if pprior
533
534
535
    aunstname='prior_unstable'; aunsttitle='Prior StabMap: explosiveness of solution';
    aindname='prior_indeterm'; aindtitle='Prior StabMap: Indeterminacy';
    awrongname='prior_wrong'; awrongtitle='Prior StabMap: inability to find solution';
536
    acalibname='prior_calib'; acalibtitle='Prior StabMap: IRF/moment restrictions';
537
    asname='prior_stable'; atitle='Prior StabMap: Parameter driving non-existence of unique stable solution (Unacceptable)';
538
else
539
540
541
    aunstname='mc_unstable'; aunsttitle='MC (around posterior mode) StabMap: explosiveness of solution';
    aindname='mc_indeterm';  aindtitle='MC (around posterior mode) StabMap: Indeterminacy';
    awrongname='mc_wrong'; awrongtitle='MC (around posterior mode) StabMap: inability to find solution';
542
    acalibname='mc_calib'; acalibtitle='MC (around posterior mode) StabMap: IRF/moment restrictions';
543
    asname='mc_stable'; atitle='MC (around posterior mode) StabMap: Parameter driving non-existence of unique stable solution (Unacceptable)';
544
end
545
546
547
548
549
delete([OutputDirectoryName,filesep,fname_,'_',asname,'.*']);
delete([OutputDirectoryName,filesep,fname_,'_',acalibname,'.*']);
delete([OutputDirectoryName,filesep,fname_,'_',aindname,'.*']);
delete([OutputDirectoryName,filesep,fname_,'_',aunstname,'.*']);
delete([OutputDirectoryName,filesep,fname_,'_',awrongname,'.*']);
550

551
if length(iunstable)>0 || length(iwrong)>0
552
    fprintf(['%4.1f%% of the prior support gives unique saddle-path solution.\n'],length(istable)/Nsam*100)
553
    fprintf(['%4.1f%% of the prior support gives explosive dynamics.\n'],(length(ixun) )/Nsam*100)
554
    if ~isempty(iindeterm)
555
        fprintf(['%4.1f%% of the prior support gives indeterminacy.\n'],length(iindeterm)/Nsam*100)
556
    end
557
    inorestriction = istable(find(~ismember(istable,irestriction))); % violation of prior restrictions
558
    if ~isempty(iwrong) || ~isempty(inorestriction)
559
        skipline()
560
        if any(infox==49)
561
562
            fprintf(['%4.1f%% of the prior support violates prior restrictions.\n'],(length(inorestriction) )/Nsam*100)
        end
563
        if ~isempty(iwrong)
564
565
566
567
            skipline()
            disp(['For ',num2str(length(iwrong)/Nsam*100,'%4.1f'),'% of the prior support dynare could not find a solution.'])
            skipline()
        end
568
        if any(infox==1)
569
            disp(['    For ',num2str(length(find(infox==1))/Nsam*100,'%4.1f'),'% The model doesn''t determine the current variables uniquely.'])
570
        end
571
        if any(infox==2)
572
            disp(['    For ',num2str(length(find(infox==2))/Nsam*100,'%4.1f'),'% MJDGGES returned an error code.'])
573
        end
574
        if any(infox==6)
575
            disp(['    For ',num2str(length(find(infox==6))/Nsam*100,'%4.1f'),'% The jacobian evaluated at the deterministic steady state is complex.'])
576
        end
577
        if any(infox==19)
578
            disp(['    For ',num2str(length(find(infox==19))/Nsam*100,'%4.1f'),'% The steadystate routine has thrown an exception (inconsistent deep parameters).'])
579
        end
580
        if any(infox==20)
581
            disp(['    For ',num2str(length(find(infox==20))/Nsam*100,'%4.1f'),'% Cannot find the steady state.'])
582
        end
583
        if any(infox==21)
584
            disp(['    For ',num2str(length(find(infox==21))/Nsam*100,'%4.1f'),'% The steady state is complex.'])
585
        end
586
        if any(infox==22)
587
            disp(['    For ',num2str(length(find(infox==22))/Nsam*100,'%4.1f'),'% The steady has NaNs.'])
588
        end
589
        if any(infox==23)
590
            disp(['    For ',num2str(length(find(infox==23))/Nsam*100,'%4.1f'),'% M_.params has been updated in the steadystate routine and has complex valued scalars.'])
591
        end
592
        if any(infox==24)
593
            disp(['    For ',num2str(length(find(infox==24))/Nsam*100,'%4.1f'),'% M_.params has been updated in the steadystate routine and has some NaNs.'])
594
        end
595
        if any(infox==30)
596
            disp(['    For ',num2str(length(find(infox==30))/Nsam*100,'%4.1f'),'% Ergodic variance can''t be computed.'])
597
598
        end

599
    end
600
    skipline()
601
602
    if length(iunstable)<Nsam || length(istable)>1
        itot = [1:Nsam];
603
        isolve = itot(find(~ismember(itot,iwrong))); % dynare could find a solution
604
                                                     % Blanchard Kahn
605
        if neighborhood_width
606
607
            options_mcf.xparam1 = xparam1(nshock+1:end);
        end
608
        itmp = itot(find(~ismember(itot,istable)));
609
610
611
612
613
        options_mcf.amcf_name = asname;
        options_mcf.amcf_title = atitle;
        options_mcf.beha_title = 'unique Stable Saddle-Path';
        options_mcf.nobeha_title = 'NO unique Stable Saddle-Path';
        options_mcf.title = 'unique solution';
614
        mcf_analysis(lpmat, istable, itmp, options_mcf, options_);
615

616
        if ~isempty(iindeterm)
617
            itmp = isolve(find(~ismember(isolve,iindeterm)));
618
619
620
621
622
            options_mcf.amcf_name = aindname;
            options_mcf.amcf_title = aindtitle;
            options_mcf.beha_title = 'NO indeterminacy';
            options_mcf.nobeha_title = 'indeterminacy';
            options_mcf.title = 'indeterminacy';
623
            mcf_analysis(lpmat, itmp, iindeterm, options_mcf, options_);
624
        end
625

626
        if ~isempty(ixun)
627
            itmp = isolve(find(~ismember(isolve,ixun)));
628
629
630
631
632
            options_mcf.amcf_name = aunstname;
            options_mcf.amcf_title = aunsttitle;
            options_mcf.beha_title = 'NO explosive solution';
            options_mcf.nobeha_title = 'explosive solution';
            options_mcf.title = 'instability';
633
            mcf_analysis(lpmat, itmp, ixun, options_mcf, options_);
634
        end
635

636
637
        inorestriction = istable(find(~ismember(istable,irestriction))); % violation of prior restrictions
        iwrong = iwrong(find(~ismember(iwrong,inorestriction))); % what went wrong beyond prior restrictions
638
        if ~isempty(iwrong)
639
            itmp = itot(find(~ismember(itot,iwrong)));
640
641
642
643
644
            options_mcf.amcf_name = awrongname;
            options_mcf.amcf_title = awrongtitle;
            options_mcf.beha_title = 'NO inability to find a solution';
            options_mcf.nobeha_title = 'inability to find a solution';
            options_mcf.title = 'inability to find a solution';
645
            mcf_analysis(lpmat, itmp, iwrong, options_mcf, options_);
646
        end
647

648
649
        if ~isempty(irestriction)
            if neighborhood_width
650
651
                options_mcf.xparam1 = xparam1;
            end
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
            np=size(bayestopt_.name,1);
            name=cell(np,1);
            name_tex=cell(np,1);
            for jj=1:np
                if options_.TeX
                    [param_name_temp, param_name_tex_temp]= get_the_name(jj,options_.TeX,M_,estim_params_,options_);
                    name_tex{jj,1} = strrep(param_name_tex_temp,'$','');
                    name{jj,1} = param_name_temp;
                else
                    param_name_temp = get_the_name(jj,options_.TeX,M_,estim_params_,options_);
                    name{jj,1} = param_name_temp;
                end
            end
            if options_.TeX
                options_mcf.param_names_tex = char(name_tex);
            end
            options_mcf.param_names = char(name);
669
670
671
672
673
            options_mcf.amcf_name = acalibname;
            options_mcf.amcf_title = acalibtitle;
            options_mcf.beha_title = 'prior IRF/moment calibration';
            options_mcf.nobeha_title = 'NO prior IRF/moment calibration';
            options_mcf.title = 'prior restrictions';
674
            mcf_analysis([lpmat0 lpmat], irestriction, inorestriction, options_mcf, options_);
675
            iok = irestriction(1);
676
            x0 = [lpmat0(iok,:)'; lpmat(iok,:)'];
677
678
        else
            iok = istable(1);
679
680
            x0=0.5.*(bounds.ub(1:nshock)-bounds.lb(1:nshock))+bounds.lb(1:nshock);
            x0 = [x0; lpmat(iok,:)'];
681
        end
682

683
        M_ = set_all_parameters(x0,estim_params_,M_);
684
685
        [oo_.dr,info,M_,options_,oo_] = resol(0,M_,options_,oo_);
        %     stoch_simul([]);
686
687
688
689
    else
        disp('All parameter values in the specified ranges are not acceptable!')
        x0=[];
    end
690
else
691
692
693
694
    disp('All parameter values in the specified ranges give unique saddle-path solution,')
    disp('and match prior IRF/moment restriction(s) if any!')
    x0=0.5.*(bounds.ub(1:nshock)-bounds.lb(1:nshock))+bounds.lb(1:nshock);
    x0 = [x0; lpmat(istable(1),:)'];
695
696
697

end

698
xparam1=x0;
699
save prior_ok.mat xparam1;
700
701

options_.periods=opt.periods;
702
if isfield(opt,'nomoments')
703
704
705
706
    options_.nomoments=opt.nomoments;
end
options_.irf=opt.irf;
options_.noprint=opt.noprint;