2012-06-11 01:33:01 +02:00
function x0 = stab_map_ ( OutputDirectoryName,opt_gsa)
2011-12-09 21:13:16 +01:00
%
% 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
%
2011-12-26 16:12:33 +01:00
% USES qmc_sequence, stab_map_1, stab_map_2
2011-12-09 21:13:16 +01:00
%
2012-01-09 13:12:31 +01:00
% Written by Marco Ratto
2011-12-09 21:13:16 +01:00
% Joint Research Centre, The European Commission,
2017-05-16 15:10:20 +02:00
% marco.ratto@ec.europa.eu
2012-01-09 13:12:31 +01:00
2022-04-13 13:15:19 +02:00
% Copyright © 2012-2016 European Commission
% Copyright © 2012-2018 Dynare Team
2012-01-09 13:12:31 +01:00
%
% 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.
2011-12-09 21:13:16 +01:00
%
2012-01-09 13:12:31 +01:00
% You should have received a copy of the GNU General Public License
2021-06-09 17:33:48 +02:00
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
2011-12-09 21:13:16 +01:00
%global bayestopt_ estim_params_ dr_ options_ ys_ fname_
global bayestopt_ estim_params_ options_ oo_ M_
2012-06-11 01:33:01 +02:00
% opt_gsa=options_.opt_gsa;
2011-12-09 21:13:16 +01:00
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 ;
ys_ = oo_ . dr . ys ;
2012-11-16 20:05:13 +01:00
nspred = M_ . nspred ; %size(dr_.ghx,2);
nboth = M_ . nboth ;
nfwrd = M_ . nfwrd ;
2011-12-09 21:13:16 +01:00
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 ;
2015-06-05 19:22:03 +02:00
lpmat0 = zeros ( Nsam , 0 ) ;
2012-02-06 16:01:09 +01:00
xparam1 = [ ] ;
2011-12-09 21:13:16 +01:00
2018-11-13 17:58:42 +01:00
[ ~ , ~ , ~ , lb , ub , ~ ] = set_prior ( estim_params_ , M_ , options_ ) ; %Prepare bounds
2015-06-05 20:48:43 +02:00
if ~ isempty ( bayestopt_ ) && any ( bayestopt_ . pshape > 0 )
% Set prior bounds
2015-12-11 18:26:07 +01:00
bounds = prior_bounds ( bayestopt_ , options_ . prior_trunc ) ;
2015-06-05 20:48:43 +02:00
bounds . lb = max ( bounds . lb , lb ) ;
bounds . ub = min ( bounds . ub , ub ) ;
else % estimated parameters but no declared priors
2017-05-16 15:10:20 +02:00
% No priors are declared so Dynare will estimate the model by
% maximum likelihood with inequality constraints for the parameters.
2015-06-05 20:48:43 +02:00
bounds . lb = lb ;
bounds . ub = ub ;
2015-06-08 09:31:27 +02:00
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
2015-06-05 20:48:43 +02:00
end
2014-10-20 16:18:54 +02:00
2017-05-16 12:42:01 +02:00
if nargin == 0
2011-12-09 21:13:16 +01:00
OutputDirectoryName = ' ' ;
end
2014-11-19 19:35:44 +01:00
options_mcf . pvalue_ks = pvalue_ks ;
options_mcf . pvalue_corr = pvalue_corr ;
options_mcf . alpha2 = alpha2 ;
2016-05-31 13:30:33 +02:00
name = cell ( np , 1 ) ;
name_tex = cell ( np , 1 ) ;
2017-05-16 15:10:20 +02:00
for jj = 1 : np
2016-05-31 13:30:33 +02:00
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
2017-10-10 10:05:59 +02:00
options_mcf . param_names_tex = name_tex ;
2016-05-31 13:30:33 +02:00
end
2017-10-10 10:05:59 +02:00
options_mcf . param_names = name ;
2016-05-31 13:30:33 +02:00
2014-11-19 19:35:44 +01:00
options_mcf . fname_ = fname_ ;
options_mcf . OutputDirectoryName = OutputDirectoryName ;
2015-10-13 10:44:13 +02:00
options_mcf . xparam1 = [ ] ;
2014-11-19 19:35:44 +01:00
2011-12-09 21:13:16 +01:00
opt = options_ ;
options_ . periods = 0 ;
options_ . nomoments = 1 ;
options_ . irf = 0 ;
options_ . noprint = 1 ;
2017-05-16 12:42:01 +02:00
if fload == 0
if isfield ( dr_ , ' ghx' )
2011-12-09 21:13:16 +01:00
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 ) ;
else
2022-10-20 15:14:39 +02:00
if np < 1112 && ilptau > 0
2012-06-29 00:36:02 +02:00
[ lpmat ] = qmc_sequence ( np , int64 ( 1 ) , 0 , Nsam ) ' ;
2017-05-16 12:42:01 +02:00
if np > 30 || ilptau == 2 % scrambled lptau
for j = 1 : np
2011-12-09 21:13:16 +01:00
lpmat ( : , j ) = lpmat ( randperm ( Nsam ) , j ) ;
end
end
else %ilptau==0
2015-06-05 19:22:03 +02:00
[ lpmat ] = NaN ( Nsam , np ) ;
2017-05-16 12:42:01 +02:00
for j = 1 : np
2011-12-09 21:13:16 +01:00
lpmat ( : , j ) = randperm ( Nsam ) ' ./ ( Nsam + 1 ) ; %latin hypercube
end
end
end
2022-05-17 10:55:47 +02:00
dummy = prior_draw_gsa ( 1 ) ;
2017-05-16 12:42:01 +02:00
if pprior
for j = 1 : nshock
if opt_gsa . morris ~= 1
2011-12-09 21:13:16 +01:00
lpmat0 ( : , j ) = randperm ( Nsam ) ' ./ ( Nsam + 1 ) ; %latin hypercube
end
if opt_gsa . prior_range
2014-10-20 16:18:54 +02:00
lpmat0 ( : , j ) = lpmat0 ( : , j ) .* ( bounds . ub ( j ) - bounds . lb ( j ) ) + bounds . lb ( j ) ;
2011-12-09 21:13:16 +01:00
end
end
if opt_gsa . prior_range
2017-05-16 12:42:01 +02:00
for j = 1 : np
2022-05-17 18:29:42 +02:00
lower_bound = max ( - options_ . huge_number , bounds . lb ( j + nshock ) ) ;
upper_bound = min ( options_ . huge_number , bounds . ub ( j + nshock ) ) ;
lpmat ( : , j ) = lpmat ( : , j ) .* ( upper_bound - lower_bound ) + lower_bound ;
2011-12-09 21:13:16 +01:00
end
else
xx = prior_draw_gsa ( 0 , [ lpmat0 lpmat ] ) ;
lpmat0 = xx ( : , 1 : nshock ) ;
lpmat = xx ( : , nshock + 1 : end ) ;
clear xx ;
end
else
2017-05-16 12:42:01 +02:00
if neighborhood_width > 0 && isempty ( options_ . mode_file )
2016-05-05 16:18:17 +02:00
xparam1 = get_all_parameters ( estim_params_ , M_ ) ;
else
2017-05-16 15:10:20 +02:00
eval ( [ ' load ' options_ . mode_file ' .mat;' ] ) ;
2016-05-05 16:18:17 +02:00
end
2017-05-16 12:42:01 +02:00
if neighborhood_width > 0
for j = 1 : nshock
2016-06-07 12:22:39 +02:00
if opt_gsa . morris ~= 1
2017-05-16 15:10:20 +02:00
lpmat0 ( : , j ) = randperm ( Nsam ) ' ./ ( Nsam + 1 ) ; %latin hypercube
2016-06-07 12:22:39 +02:00
end
2014-10-20 16:18:54 +02:00
ub = min ( [ bounds . ub ( j ) xparam1 ( j ) * ( 1 + neighborhood_width ) ] ) ;
lb = max ( [ bounds . lb ( j ) xparam1 ( j ) * ( 1 - neighborhood_width ) ] ) ;
2011-12-09 21:13:16 +01:00
lpmat0 ( : , j ) = lpmat0 ( : , j ) .* ( ub - lb ) + lb ;
end
2017-05-16 12:42:01 +02:00
for j = 1 : np
2016-06-07 12:22:39 +02:00
ub = xparam1 ( j + nshock ) * ( 1 + sign ( xparam1 ( j + nshock ) ) * neighborhood_width ) ;
lb = xparam1 ( j + nshock ) * ( 1 - sign ( xparam1 ( j + nshock ) ) * neighborhood_width ) ;
2022-05-20 15:12:12 +02:00
if bounds . ub ( j + nshock ) > = xparam1 ( j + nshock ) && bounds . lb ( j + nshock ) < = xparam1 ( j + nshock )
2016-06-07 12:22:39 +02:00
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
2011-12-09 21:13:16 +01:00
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 ' ) ;
2015-06-05 19:22:03 +02:00
lnprior = zeros ( 1 , Nsam * 2 ) ;
2017-05-16 12:42:01 +02:00
for j = 1 : Nsam * 2
2014-10-20 16:18:54 +02:00
lnprior ( j ) = any ( lp ( j , : ) ' < = bounds . lb | lp ( j , : ) ' > = bounds . ub ) ;
2011-12-09 21:13:16 +01:00
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
%
2011-12-13 18:32:57 +01:00
h = dyn_waitbar ( 0 , ' Please wait...' ) ;
2011-12-09 21:13:16 +01:00
istable = [ 1 : Nsam ] ;
jstab = 0 ;
iunstable = [ 1 : Nsam ] ;
iindeterm = zeros ( 1 , Nsam ) ;
iwrong = zeros ( 1 , Nsam ) ;
2014-11-14 17:59:08 +01:00
inorestriction = zeros ( 1 , Nsam ) ;
irestriction = zeros ( 1 , Nsam ) ;
2016-05-05 16:53:18 +02:00
infox = zeros ( Nsam , 1 ) ;
2017-05-16 12:42:01 +02:00
for j = 1 : Nsam
2014-09-04 16:39:59 +02:00
M_ = set_all_parameters ( [ lpmat0 ( j , : ) lpmat ( j , : ) ] ' , estim_params_ , M_ ) ;
2011-12-09 21:13:16 +01:00
%try stoch_simul([]);
try
2015-04-14 12:24:49 +02:00
if ~ isempty ( options_ . endogenous_prior_restrictions . moment )
2021-06-22 14:28:46 +02:00
[ Tt , Rr , SteadyState , info , M_ , oo_ ] = dynare_resolve ( M_ , options_ , oo_ ) ;
2015-04-14 12:24:49 +02:00
else
2021-06-22 14:28:46 +02:00
[ Tt , Rr , SteadyState , info , M_ , oo_ ] = dynare_resolve ( M_ , options_ , oo_ , ' restrict' ) ;
2015-04-14 12:24:49 +02:00
end
2013-03-19 17:04:33 +01:00
infox ( j , 1 ) = info ( 1 ) ;
2017-05-16 12:42:01 +02:00
if infox ( j , 1 ) == 0 && ~ exist ( ' T' , ' var' )
2011-12-09 21:13:16 +01:00
dr_ = oo_ . dr ;
2017-05-16 12:42:01 +02:00
if prepSA
2013-03-19 17:04:33 +01:00
try
T = zeros ( size ( dr_ . ghx , 1 ) , size ( dr_ . ghx , 2 ) + size ( dr_ . ghu , 2 ) , Nsam ) ;
2013-11-08 16:35:52 +01:00
catch
ME = lasterror ( ) ;
2017-05-16 12:42:01 +02:00
if strcmp ( ' MATLAB:nomem' , ME . identifier )
2013-03-19 17:04:33 +01:00
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
2011-12-09 21:13:16 +01:00
egg = zeros ( length ( dr_ . eigval ) , Nsam ) ;
end
2017-05-16 12:42:01 +02:00
if infox ( j , 1 )
2017-05-16 15:10:20 +02:00
% disp('no solution'),
2017-05-16 12:42:01 +02:00
if isfield ( oo_ . dr , ' ghx' )
2011-12-09 21:13:16 +01:00
oo_ . dr = rmfield ( oo_ . dr , ' ghx' ) ;
end
2017-05-16 12:42:01 +02:00
if ( infox ( j , 1 ) < 3 || infox ( j , 1 ) > 5 ) && isfield ( oo_ . dr , ' eigval' )
2013-02-08 14:31:49 +01:00
oo_ . dr = rmfield ( oo_ . dr , ' eigval' ) ;
end
2011-12-09 21:13:16 +01:00
end
2013-03-19 17:04:33 +01:00
catch ME
2017-05-16 12:42:01 +02:00
if isfield ( oo_ . dr , ' eigval' )
2011-12-09 21:13:16 +01:00
oo_ . dr = rmfield ( oo_ . dr , ' eigval' ) ;
end
2017-05-16 12:42:01 +02:00
if isfield ( oo_ . dr , ' ghx' )
2011-12-09 21:13:16 +01:00
oo_ . dr = rmfield ( oo_ . dr , ' ghx' ) ;
end
2017-05-16 12:42:01 +02:00
disp ( ' No solution could be found' )
2011-12-09 21:13:16 +01:00
end
dr_ = oo_ . dr ;
2017-05-16 12:42:01 +02:00
if isfield ( dr_ , ' ghx' )
2011-12-09 21:13:16 +01:00
egg ( : , j ) = sort ( dr_ . eigval ) ;
if prepSA
jstab = jstab + 1 ;
T ( : , : , jstab ) = [ dr_ . ghx dr_ . ghu ] ;
end
2017-05-16 12:42:01 +02:00
if ~ exist ( ' nspred' , ' var' )
2011-12-09 21:13:16 +01:00
nspred = dr_ . nspred ; %size(dr_.ghx,2);
nboth = dr_ . nboth ;
nfwrd = dr_ . nfwrd ;
end
2014-11-14 17:59:08 +01:00
info = endogenous_prior_restrictions ( Tt , Rr , M_ , options_ , oo_ ) ;
infox ( j , 1 ) = info ( 1 ) ;
2017-05-16 12:42:01 +02:00
if info ( 1 )
2014-11-14 17:59:08 +01:00
inorestriction ( j ) = j ;
else
iunstable ( j ) = 0 ;
irestriction ( j ) = j ;
end
2011-12-09 21:13:16 +01:00
else
istable ( j ) = 0 ;
if isfield ( dr_ , ' eigval' )
egg ( : , j ) = sort ( dr_ . eigval ) ;
2015-06-05 19:22:03 +02:00
if exist ( ' nspred' , ' var' )
2011-12-09 21:13:16 +01:00
if any ( isnan ( egg ( 1 : nspred , j ) ) )
iwrong ( j ) = j ;
else
2017-05-16 12:42:01 +02:00
if ( nboth || nfwrd ) && abs ( egg ( nspred + 1 , j ) ) < = options_ . qz_criterium
2011-12-09 21:13:16 +01:00
iindeterm ( j ) = j ;
end
end
end
else
2017-05-16 12:42:01 +02:00
if exist ( ' egg' , ' var' )
2011-12-09 21:13:16 +01:00
egg ( : , j ) = ones ( size ( egg , 1 ) , 1 ) .* NaN ;
end
iwrong ( j ) = j ;
end
end
ys_ = real ( dr_ . ys ) ;
yys ( : , j ) = ys_ ;
ys_ = yys ( : , 1 ) ;
2013-04-27 16:53:08 +02:00
dyn_waitbar ( j / Nsam , h , [ ' MC iteration ' , int2str ( j ) , ' /' , int2str ( Nsam ) ] )
2011-12-09 21:13:16 +01:00
end
2011-12-13 18:32:57 +01:00
dyn_waitbar_close ( h ) ;
2017-05-16 12:42:01 +02:00
if prepSA && jstab
2011-12-09 21:13:16 +01:00
T = T ( : , : , 1 : jstab ) ;
2013-03-19 17:04:33 +01:00
else
T = [ ] ;
2011-12-09 21:13:16 +01:00
end
2014-11-14 17:59:08 +01:00
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)
2011-12-09 21:13:16 +01:00
iindeterm = iindeterm ( find ( iindeterm ) ) ; % indeterminacy
iwrong = iwrong ( find ( iwrong ) ) ; % dynare could not find solution
2016-05-05 16:53:18 +02:00
ixun = iunstable ( find ( ~ ismember ( iunstable , [ iindeterm , iwrong , inorestriction ] ) ) ) ; % explosive roots
2011-12-09 21:13:16 +01:00
bkpprior . pshape = bayestopt_ . pshape ;
bkpprior . p1 = bayestopt_ . p1 ;
bkpprior . p2 = bayestopt_ . p2 ;
bkpprior . p3 = bayestopt_ . p3 ;
bkpprior . p4 = bayestopt_ . p4 ;
2017-05-16 12:42:01 +02:00
if pprior
2011-12-09 21:13:16 +01:00
if ~ prepSA
2013-04-27 16:22:25 +02:00
save ( [ OutputDirectoryName filesep fname_ ' _prior.mat' ] , ...
2017-05-16 15:10:20 +02:00
' bkpprior' , ' lpmat' , ' lpmat0' , ' irestriction' , ' iunstable' , ' istable' , ' iindeterm' , ' iwrong' , ' ixun' , ...
' egg' , ' yys' , ' nspred' , ' nboth' , ' nfwrd' , ' infox' )
2011-12-09 21:13:16 +01:00
else
2013-04-27 16:22:25 +02:00
save ( [ OutputDirectoryName filesep fname_ ' _prior.mat' ] , ...
2017-05-16 15:10:20 +02:00
' bkpprior' , ' lpmat' , ' lpmat0' , ' irestriction' , ' iunstable' , ' istable' , ' iindeterm' , ' iwrong' , ' ixun' , ...
' egg' , ' yys' , ' T' , ' nspred' , ' nboth' , ' nfwrd' , ' infox' )
2011-12-09 21:13:16 +01:00
end
else
if ~ prepSA
2013-04-27 16:22:25 +02:00
save ( [ OutputDirectoryName filesep fname_ ' _mc.mat' ] , ...
2017-05-16 15:10:20 +02:00
' lpmat' , ' lpmat0' , ' irestriction' , ' iunstable' , ' istable' , ' iindeterm' , ' iwrong' , ' ixun' , ...
' egg' , ' yys' , ' nspred' , ' nboth' , ' nfwrd' , ' infox' )
2011-12-09 21:13:16 +01:00
else
2013-04-27 16:22:25 +02:00
save ( [ OutputDirectoryName filesep fname_ ' _mc.mat' ] , ...
2017-05-16 15:10:20 +02:00
' lpmat' , ' lpmat0' , ' irestriction' , ' iunstable' , ' istable' , ' iindeterm' , ' iwrong' , ' ixun' , ...
' egg' , ' yys' , ' T' , ' nspred' , ' nboth' , ' nfwrd' , ' infox' )
2011-12-09 21:13:16 +01:00
end
end
else
2017-05-16 12:42:01 +02:00
if pprior
2013-04-27 16:22:25 +02:00
filetoload = [ OutputDirectoryName filesep fname_ ' _prior.mat' ] ;
2011-12-09 21:13:16 +01:00
else
2013-04-27 16:22:25 +02:00
filetoload = [ OutputDirectoryName filesep fname_ ' _mc.mat' ] ;
2011-12-09 21:13:16 +01:00
end
2014-11-14 17:59:08 +01:00
load ( filetoload , ' lpmat' , ' lpmat0' , ' irestriction' , ' iunstable' , ' istable' , ' iindeterm' , ' iwrong' , ' ixun' , ' egg' , ' yys' , ' nspred' , ' nboth' , ' nfwrd' , ' infox' )
2011-12-09 21:13:16 +01:00
Nsam = size ( lpmat , 1 ) ;
2017-05-16 12:42:01 +02:00
if pprior == 0 && ~ isempty ( options_ . mode_file )
2012-02-06 16:01:09 +01:00
eval ( [ ' load ' options_ . mode_file ' .mat;' ] ) ;
end
2011-12-09 21:13:16 +01:00
2017-05-16 12:42:01 +02:00
if prepSA && isempty ( strmatch ( ' T' , who ( ' -file' , filetoload ) , ' exact' ) )
2011-12-13 18:32:57 +01:00
h = dyn_waitbar ( 0 , ' Please wait...' ) ;
2011-12-09 21:13:16 +01:00
options_ . periods = 0 ;
options_ . nomoments = 1 ;
options_ . irf = 0 ;
options_ . noprint = 1 ;
2019-09-10 17:02:20 +02:00
[ ~ , oo_ , options_ ] = stoch_simul ( M_ , options_ , oo_ , [ ] ) ;
2011-12-09 21:13:16 +01:00
ntrans = length ( istable ) ;
2015-06-05 19:22:03 +02:00
yys = NaN ( length ( ys_ ) , ntrans ) ;
2017-05-16 12:42:01 +02:00
for j = 1 : ntrans
2011-12-09 21:13:16 +01:00
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' ) ;
2015-06-05 19:22:03 +02:00
if ~ exist ( ' T' , ' var' )
2011-12-09 21:13:16 +01:00
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 ] ;
2015-06-05 19:22:03 +02:00
if ~ exist ( ' nspred' , ' var' )
2011-12-09 21:13:16 +01:00
nspred = dr_ . nspred ; %size(dr_.ghx,2);
nboth = dr_ . nboth ;
nfwrd = dr_ . nfwrd ;
end
ys_ = real ( dr_ . ys ) ;
yys ( : , j ) = ys_ ;
ys_ = yys ( : , 1 ) ;
2013-04-27 16:53:08 +02:00
dyn_waitbar ( j / ntrans , h , [ ' MC iteration ' , int2str ( j ) , ' /' , int2str ( ntrans ) ] )
2011-12-09 21:13:16 +01:00
end
2011-12-13 18:32:57 +01:00
dyn_waitbar_close ( h ) ;
2011-12-09 21:13:16 +01:00
save ( filetoload , ' T' , ' -append' )
elseif prepSA
load ( filetoload , ' T' )
end
end
if pprior
2013-12-10 14:32:09 +01:00
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' ;
2014-11-14 17:59:08 +01:00
acalibname = ' prior_calib' ; acalibtitle = ' Prior StabMap: IRF/moment restrictions' ;
2014-11-27 11:32:03 +01:00
asname = ' prior_stable' ; atitle = ' Prior StabMap: Parameter driving non-existence of unique stable solution (Unacceptable)' ;
2011-12-09 21:13:16 +01:00
else
2013-12-10 14:32:09 +01:00
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' ;
2014-11-14 17:59:08 +01:00
acalibname = ' mc_calib' ; acalibtitle = ' MC (around posterior mode) StabMap: IRF/moment restrictions' ;
2014-11-27 11:32:03 +01:00
asname = ' mc_stable' ; atitle = ' MC (around posterior mode) StabMap: Parameter driving non-existence of unique stable solution (Unacceptable)' ;
2011-12-09 21:13:16 +01:00
end
2014-11-27 11:32:03 +01:00
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 , ' .*' ] ) ;
2011-12-09 21:13:16 +01:00
2017-05-16 12:42:01 +02:00
if length ( iunstable ) > 0 || length ( iwrong ) > 0
2013-03-19 17:04:33 +01:00
fprintf ( [ ' %4.1f%% of the prior support gives unique saddle-path solution.\n' ] , length ( istable ) / Nsam * 100 )
2016-05-05 16:53:18 +02:00
fprintf ( [ ' %4.1f%% of the prior support gives explosive dynamics.\n' ] , ( length ( ixun ) ) / Nsam * 100 )
2017-05-16 12:42:01 +02:00
if ~ isempty ( iindeterm )
2016-05-05 16:53:18 +02:00
fprintf ( [ ' %4.1f%% of the prior support gives indeterminacy.\n' ] , length ( iindeterm ) / Nsam * 100 )
2011-12-09 21:13:16 +01:00
end
2016-05-05 16:53:18 +02:00
inorestriction = istable ( find ( ~ ismember ( istable , irestriction ) ) ) ; % violation of prior restrictions
2017-05-16 12:42:01 +02:00
if ~ isempty ( iwrong ) || ~ isempty ( inorestriction )
2013-07-10 17:12:34 +02:00
skipline ( )
2017-05-16 12:42:01 +02:00
if any ( infox == 49 )
2016-05-05 16:53:18 +02:00
fprintf ( [ ' %4.1f%% of the prior support violates prior restrictions.\n' ] , ( length ( inorestriction ) ) / Nsam * 100 )
end
2017-05-16 12:42:01 +02:00
if ~ isempty ( iwrong )
2016-05-05 16:53:18 +02:00
skipline ( )
disp ( [ ' For ' , num2str ( length ( iwrong ) / Nsam * 100 , ' %4.1f' ) , ' % of the prior support dynare could not find a solution.' ] )
skipline ( )
end
2017-05-16 12:42:01 +02:00
if any ( infox == 1 )
2014-09-04 16:39:59 +02:00
disp ( [ ' For ' , num2str ( length ( find ( infox == 1 ) ) / Nsam * 100 , ' %4.1f' ) , ' % The model doesn' ' t determine the current variables uniquely.' ] )
2013-03-19 17:04:33 +01:00
end
2017-05-16 12:42:01 +02:00
if any ( infox == 2 )
2014-09-04 16:39:59 +02:00
disp ( [ ' For ' , num2str ( length ( find ( infox == 2 ) ) / Nsam * 100 , ' %4.1f' ) , ' % MJDGGES returned an error code.' ] )
2013-03-19 17:04:33 +01:00
end
2017-05-16 12:42:01 +02:00
if any ( infox == 6 )
2014-09-04 16:39:59 +02:00
disp ( [ ' For ' , num2str ( length ( find ( infox == 6 ) ) / Nsam * 100 , ' %4.1f' ) , ' % The jacobian evaluated at the deterministic steady state is complex.' ] )
2013-03-19 17:04:33 +01:00
end
2017-05-16 12:42:01 +02:00
if any ( infox == 19 )
2017-11-28 10:27:43 +01:00
disp ( [ ' For ' , num2str ( length ( find ( infox == 19 ) ) / Nsam * 100 , ' %4.1f' ) , ' % The steadystate routine has thrown an exception (inconsistent deep parameters).' ] )
2013-03-19 17:04:33 +01:00
end
2017-05-16 12:42:01 +02:00
if any ( infox == 20 )
2014-09-04 16:39:59 +02:00
disp ( [ ' For ' , num2str ( length ( find ( infox == 20 ) ) / Nsam * 100 , ' %4.1f' ) , ' % Cannot find the steady state.' ] )
2013-03-19 17:04:33 +01:00
end
2017-05-16 12:42:01 +02:00
if any ( infox == 21 )
2014-09-04 16:39:59 +02:00
disp ( [ ' For ' , num2str ( length ( find ( infox == 21 ) ) / Nsam * 100 , ' %4.1f' ) , ' % The steady state is complex.' ] )
2013-03-19 17:04:33 +01:00
end
2017-05-16 12:42:01 +02:00
if any ( infox == 22 )
2014-09-04 16:39:59 +02:00
disp ( [ ' For ' , num2str ( length ( find ( infox == 22 ) ) / Nsam * 100 , ' %4.1f' ) , ' % The steady has NaNs.' ] )
2013-03-19 17:04:33 +01:00
end
2017-05-16 12:42:01 +02:00
if any ( infox == 23 )
2014-09-04 16:39:59 +02:00
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.' ] )
2013-03-19 17:04:33 +01:00
end
2017-05-16 12:42:01 +02:00
if any ( infox == 24 )
2014-09-04 16:39:59 +02:00
disp ( [ ' For ' , num2str ( length ( find ( infox == 24 ) ) / Nsam * 100 , ' %4.1f' ) , ' % M_.params has been updated in the steadystate routine and has some NaNs.' ] )
2013-03-19 17:04:33 +01:00
end
2017-05-16 12:42:01 +02:00
if any ( infox == 30 )
2014-09-04 16:39:59 +02:00
disp ( [ ' For ' , num2str ( length ( find ( infox == 30 ) ) / Nsam * 100 , ' %4.1f' ) , ' % Ergodic variance can' ' t be computed.' ] )
2013-03-19 17:04:33 +01:00
end
2011-12-09 21:13:16 +01:00
end
2013-07-10 17:12:34 +02:00
skipline ( )
2014-11-14 17:59:08 +01:00
if length ( iunstable ) < Nsam || length ( istable ) > 1
itot = [ 1 : Nsam ] ;
2016-05-05 16:53:18 +02:00
isolve = itot ( find ( ~ ismember ( itot , iwrong ) ) ) ; % dynare could find a solution
2017-05-16 15:10:20 +02:00
% Blanchard Kahn
2017-05-16 12:42:01 +02:00
if neighborhood_width
2015-10-13 10:44:13 +02:00
options_mcf . xparam1 = xparam1 ( nshock + 1 : end ) ;
end
2014-11-14 17:59:08 +01:00
itmp = itot ( find ( ~ ismember ( itot , istable ) ) ) ;
2014-11-19 19:35:44 +01:00
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' ;
2015-04-16 16:09:39 +02:00
mcf_analysis ( lpmat , istable , itmp , options_mcf , options_ ) ;
2014-11-19 19:35:44 +01:00
2017-05-16 12:42:01 +02:00
if ~ isempty ( iindeterm )
2016-05-05 16:53:18 +02:00
itmp = isolve ( find ( ~ ismember ( isolve , iindeterm ) ) ) ;
2014-11-19 19:35:44 +01:00
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' ;
2015-04-16 16:09:39 +02:00
mcf_analysis ( lpmat , itmp , iindeterm , options_mcf , options_ ) ;
2011-12-09 21:13:16 +01:00
end
2017-05-16 15:10:20 +02:00
2017-05-16 12:42:01 +02:00
if ~ isempty ( ixun )
2016-05-05 16:53:18 +02:00
itmp = isolve ( find ( ~ ismember ( isolve , ixun ) ) ) ;
2014-11-19 19:35:44 +01:00
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' ;
2015-04-16 16:09:39 +02:00
mcf_analysis ( lpmat , itmp , ixun , options_mcf , options_ ) ;
2011-12-09 21:13:16 +01:00
end
2017-05-16 15:10:20 +02:00
2016-05-05 16:53:18 +02:00
inorestriction = istable ( find ( ~ ismember ( istable , irestriction ) ) ) ; % violation of prior restrictions
iwrong = iwrong ( find ( ~ ismember ( iwrong , inorestriction ) ) ) ; % what went wrong beyond prior restrictions
2017-05-16 12:42:01 +02:00
if ~ isempty ( iwrong )
2014-11-14 17:59:08 +01:00
itmp = itot ( find ( ~ ismember ( itot , iwrong ) ) ) ;
2014-11-19 19:35:44 +01:00
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' ;
2015-04-16 16:09:39 +02:00
mcf_analysis ( lpmat , itmp , iwrong , options_mcf , options_ ) ;
2011-12-09 21:13:16 +01:00
end
2017-05-16 15:10:20 +02:00
2017-05-16 12:42:01 +02:00
if ~ isempty ( irestriction )
if neighborhood_width
2015-10-13 10:44:13 +02:00
options_mcf . xparam1 = xparam1 ;
end
2016-06-19 14:32:33 +02:00
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
2017-10-10 10:05:59 +02:00
options_mcf . param_names_tex = name_tex ;
2016-06-19 14:32:33 +02:00
end
2017-10-10 10:05:59 +02:00
options_mcf . param_names = name ;
2014-11-19 19:35:44 +01:00
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' ;
2015-04-16 16:09:39 +02:00
mcf_analysis ( [ lpmat0 lpmat ] , irestriction , inorestriction , options_mcf , options_ ) ;
2014-11-14 17:59:08 +01:00
iok = irestriction ( 1 ) ;
2014-11-19 19:35:44 +01:00
x0 = [ lpmat0 ( iok , : ) ' ; lpmat ( iok , : ) ' ] ;
2014-11-14 17:59:08 +01:00
else
iok = istable ( 1 ) ;
2014-11-19 19:35:44 +01:00
x0 = 0.5 .* ( bounds . ub ( 1 : nshock ) - bounds . lb ( 1 : nshock ) ) + bounds . lb ( 1 : nshock ) ;
x0 = [ x0 ; lpmat ( iok , : ) ' ] ;
2014-09-04 16:39:59 +02:00
end
2017-05-16 15:10:20 +02:00
2014-11-19 19:35:44 +01:00
M_ = set_all_parameters ( x0 , estim_params_ , M_ ) ;
2021-06-22 14:28:46 +02:00
[ oo_ . dr , info , M_ , oo_ ] = resol ( 0 , M_ , options_ , oo_ ) ;
2011-12-09 21:13:16 +01:00
else
disp ( ' All parameter values in the specified ranges are not acceptable!' )
x0 = [ ] ;
end
2014-09-04 16:39:59 +02:00
else
2014-11-14 17:59:08 +01:00
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 ) , : ) ' ] ;
2011-12-09 21:13:16 +01:00
end
2012-06-11 01:33:01 +02:00
xparam1 = x0 ;
2021-09-17 10:47:01 +02:00
save ( [ OutputDirectoryName filesep ' prior_ok.mat' ] , ' xparam1' ) ;
2011-12-09 21:13:16 +01:00
options_ . periods = opt . periods ;
2017-05-16 12:42:01 +02:00
if isfield ( opt , ' nomoments' )
2011-12-09 21:13:16 +01:00
options_ . nomoments = opt . nomoments ;
end
options_ . irf = opt . irf ;
options_ . noprint = opt . noprint ;