2010-02-10 18:51:07 +01:00
function myoutput= prior_posterior_statistics_core ( myinputs,fpar,B,whoiam, ThisMatlab)
2010-05-31 11:54:22 +02:00
% PARALLEL CONTEXT
% Core functionality for prior_posterior.m function, which can be parallelized.
2016-05-19 17:13:39 +02:00
% See also the comment in posterior_sampler_core.m funtion.
2010-05-31 11:54:22 +02:00
%
2011-02-02 14:13:11 +01:00
% INPUTS
2016-05-19 17:13:39 +02:00
% See the comment in posterior_sampler_core.m funtion.
2011-02-04 17:27:33 +01:00
%
2010-05-31 11:54:22 +02:00
% OUTPUTS
% o myoutput [struc]
2011-02-02 14:13:11 +01:00
% Contained OutputFileName_smooth;
% _update;
% _inno;
% _error;
% _filter_step_ahead;
2010-05-31 11:54:22 +02:00
% _param;
% _forc_mean;
2016-06-03 17:09:09 +02:00
% _forc_point;
2016-06-18 14:33:11 +02:00
% _forc_point_ME;
2016-06-03 17:09:09 +02:00
% _filter_covar;
% _trend_coeff;
% _smoothed_trend;
% _smoothed_constant;
2016-11-02 15:55:23 +01:00
% _state_uncert;
2010-02-15 16:52:36 +01:00
%
2011-02-02 14:13:11 +01:00
% ALGORITHM
% Portion of prior_posterior.m function.
2010-02-15 16:52:36 +01:00
% This file is part of Dynare.
%
2010-05-31 11:54:22 +02:00
% SPECIAL REQUIREMENTS.
% None.
2011-02-04 17:27:33 +01:00
2022-04-13 13:15:19 +02:00
% Copyright © 2005-2020 Dynare Team
2011-02-04 17:27:33 +01:00
%
% This file is part of Dynare.
%
2010-02-15 16:52:36 +01:00
% 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
2021-06-09 17:33:48 +02:00
% along with Dynare. If not, see <https://www.gnu.org/licenses/>.
2010-02-10 18:51:07 +01:00
2010-04-02 21:37:43 +02:00
global options_ oo_ M_ bayestopt_ estim_params_
2017-05-16 12:42:01 +02:00
if nargin < 4
2010-02-10 18:51:07 +01:00
whoiam = 0 ;
end
2010-05-31 11:54:22 +02:00
% Reshape 'myinputs' for local computation.
% In order to avoid confusion in the name space, the instruction struct2local(myinputs) is replaced by:
type = myinputs . type ;
2010-04-02 21:37:43 +02:00
run_smoother = myinputs . run_smoother ;
2016-06-02 18:28:36 +02:00
filter_covariance = myinputs . filter_covariance ;
2016-11-02 15:55:23 +01:00
smoothed_state_uncertainty = myinputs . smoothed_state_uncertainty ;
2010-04-02 21:37:43 +02:00
gend = myinputs . gend ;
Y = myinputs . Y ;
data_index = myinputs . data_index ;
missing_value = myinputs . missing_value ;
varobs = myinputs . varobs ;
2015-03-02 10:45:12 +01:00
mean_varobs = myinputs . mean_varobs ;
2010-04-02 21:37:43 +02:00
irun = myinputs . irun ;
endo_nbr = myinputs . endo_nbr ;
nvn = myinputs . nvn ;
naK = myinputs . naK ;
horizon = myinputs . horizon ;
iendo = myinputs . iendo ;
2015-03-11 18:15:28 +01:00
IdObs = myinputs . IdObs ; %index of observables
2010-04-02 21:37:43 +02:00
if horizon
i_last_obs = myinputs . i_last_obs ;
2010-05-31 11:54:22 +02:00
MAX_nforc1 = myinputs . MAX_nforc1 ;
MAX_nforc2 = myinputs . MAX_nforc2 ;
2016-06-18 14:33:11 +02:00
if ~ isequal ( M_ . H , 0 )
MAX_nforc_ME = myinputs . MAX_nforc_ME ;
end
2010-04-02 21:37:43 +02:00
end
2010-05-31 11:54:22 +02:00
if naK
MAX_naK = myinputs . MAX_naK ;
end
2016-06-02 18:28:36 +02:00
if filter_covariance
2016-11-02 15:55:23 +01:00
MAX_filter_covariance = myinputs . MAX_filter_covariance ;
2016-06-02 18:28:36 +02:00
end
2016-11-02 15:55:23 +01:00
if smoothed_state_uncertainty
MAX_n_smoothed_state_uncertainty = myinputs . MAX_n_smoothed_state_uncertainty ;
end
2010-04-02 21:37:43 +02:00
exo_nbr = myinputs . exo_nbr ;
maxlag = myinputs . maxlag ;
MAX_nsmoo = myinputs . MAX_nsmoo ;
MAX_ninno = myinputs . MAX_ninno ;
2016-06-03 17:09:09 +02:00
MAX_n_smoothed_constant = myinputs . MAX_n_smoothed_constant ;
MAX_n_smoothed_trend = myinputs . MAX_n_smoothed_trend ;
MAX_n_trend_coeff = myinputs . MAX_n_trend_coeff ;
2010-05-31 11:54:22 +02:00
MAX_nerro = myinputs . MAX_nerro ;
2010-04-02 21:37:43 +02:00
MAX_nruns = myinputs . MAX_nruns ;
2010-05-31 11:54:22 +02:00
MAX_momentsno = myinputs . MAX_momentsno ;
2010-04-02 21:37:43 +02:00
ifil = myinputs . ifil ;
2010-02-10 18:51:07 +01:00
2017-05-16 12:42:01 +02:00
if ~ strcmpi ( type , ' prior' )
2010-04-02 21:37:43 +02:00
x = myinputs . x ;
2017-05-16 12:42:01 +02:00
if strcmpi ( type , ' posterior' )
2011-02-25 09:27:33 +01:00
logpost = myinputs . logpost ;
end
2021-05-27 16:13:58 +02:00
else
x = ( prior_draw ( bayestopt_ , options_ . prior_trunc ) ) ' ;
options_ = select_qz_criterium_value ( options_ ) ;
2010-04-02 21:37:43 +02:00
end
2010-05-31 11:54:22 +02:00
if whoiam
Parallel = myinputs . Parallel ;
end
2010-02-10 18:51:07 +01:00
2011-04-06 14:24:04 +02:00
% DirectoryName = myinputs.DirectoryName;
if strcmpi ( type , ' posterior' )
2011-12-15 11:56:23 +01:00
DirectoryName = CheckPath ( ' metropolis' , M_ . dname ) ;
2011-04-06 14:24:04 +02:00
elseif strcmpi ( type , ' gsa' )
if options_ . opt_gsa . pprior
2011-12-15 11:56:23 +01:00
DirectoryName = CheckPath ( [ ' gsa' , filesep , ' prior' ] , M_ . dname ) ;
2011-04-06 14:24:04 +02:00
else
2011-12-15 11:56:23 +01:00
DirectoryName = CheckPath ( [ ' gsa' , filesep , ' mc' ] , M_ . dname ) ;
2011-04-06 14:24:04 +02:00
end
elseif strcmpi ( type , ' prior' )
2011-12-15 11:56:23 +01:00
DirectoryName = CheckPath ( ' prior' , M_ . dname ) ;
2011-04-06 14:24:04 +02:00
end
2010-02-10 18:51:07 +01:00
RemoteFlag = 0 ;
if whoiam
2017-05-16 12:42:01 +02:00
if Parallel ( ThisMatlab ) . Local == 0
2011-12-13 18:32:57 +01:00
RemoteFlag = 1 ;
2010-10-12 08:07:57 +02:00
end
2010-02-10 18:51:07 +01:00
ifil = ifil ( : , whoiam ) ;
2011-12-13 18:32:57 +01:00
prct0 = { 0 , whoiam , Parallel ( ThisMatlab ) } ;
2010-02-12 17:53:43 +01:00
else
2011-12-13 18:32:57 +01:00
prct0 = 0 ;
2010-02-10 18:51:07 +01:00
end
2011-12-13 18:32:57 +01:00
h = dyn_waitbar ( prct0 , [ ' Taking ' , type , ' subdraws...' ] ) ;
2010-02-10 18:51:07 +01:00
2017-05-16 12:42:01 +02:00
if RemoteFlag == 1
2011-02-02 14:13:11 +01:00
OutputFileName_smooth = { } ;
OutputFileName_update = { } ;
OutputFileName_inno = { } ;
OutputFileName_error = { } ;
OutputFileName_filter_step_ahead = { } ;
OutputFileName_param = { } ;
OutputFileName_forc_mean = { } ;
OutputFileName_forc_point = { } ;
2016-06-18 14:33:11 +02:00
OutputFileName_forc_point_ME = { } ;
2016-06-02 18:28:36 +02:00
OutputFileName_filter_covar = { } ;
2016-06-03 17:09:09 +02:00
OutputFileName_trend_coeff = { } ;
OutputFileName_smoothed_trend = { } ;
OutputFileName_smoothed_constant = { } ;
2011-02-02 14:13:11 +01:00
% OutputFileName_moments = {};
2010-02-10 18:51:07 +01:00
end
2013-03-01 17:52:11 +01:00
%initialize arrays
if run_smoother
2017-05-16 15:10:20 +02:00
stock_smooth = NaN ( endo_nbr , gend , MAX_nsmoo ) ;
stock_update = NaN ( endo_nbr , gend , MAX_nsmoo ) ;
stock_innov = NaN ( M_ . exo_nbr , gend , MAX_ninno ) ;
stock_smoothed_constant = NaN ( endo_nbr , gend , MAX_n_smoothed_constant ) ;
stock_smoothed_trend = NaN ( endo_nbr , gend , MAX_n_smoothed_trend ) ;
stock_trend_coeff = zeros ( endo_nbr , MAX_n_trend_coeff ) ;
if horizon
stock_forcst_mean = NaN ( endo_nbr , horizon , MAX_nforc1 ) ;
stock_forcst_point = NaN ( endo_nbr , horizon , MAX_nforc2 ) ;
if ~ isequal ( M_ . H , 0 )
stock_forcst_point_ME = NaN ( length ( varobs ) , horizon , MAX_nforc_ME ) ;
end
end
2013-03-01 17:52:11 +01:00
end
if nvn
2017-05-16 15:10:20 +02:00
stock_error = NaN ( length ( varobs ) , gend , MAX_nerro ) ;
2013-03-01 17:52:11 +01:00
end
if naK
stock_filter_step_ahead = NaN ( length ( options_ . filter_step_ahead ) , endo_nbr , gend + max ( options_ . filter_step_ahead ) , MAX_naK ) ;
end
2021-05-27 16:13:58 +02:00
stock_param = NaN ( MAX_nruns , size ( x , 2 ) ) ;
2013-03-01 17:52:11 +01:00
stock_logpo = NaN ( MAX_nruns , 1 ) ;
stock_ys = NaN ( MAX_nruns , endo_nbr ) ;
2016-06-02 18:28:36 +02:00
if filter_covariance
stock_filter_covariance = zeros ( endo_nbr , endo_nbr , gend + 1 , MAX_filter_covariance ) ;
end
2016-11-02 15:55:23 +01:00
if smoothed_state_uncertainty
stock_smoothed_uncert = zeros ( endo_nbr , endo_nbr , gend , MAX_n_smoothed_state_uncertainty ) ;
end
2013-03-01 17:52:11 +01:00
2021-05-27 16:13:58 +02:00
opts_local = options_ ;
2010-02-10 18:51:07 +01:00
for b = fpar : B
2010-05-31 11:54:22 +02:00
if strcmpi ( type , ' prior' )
2021-05-27 16:13:58 +02:00
iter = 1 ;
logpo = [ ] ;
while ( isempty ( logpo ) || isinf ( logpo ) ) && iter < 1000
[ deep , logpo ] = GetOneDraw ( type , M_ , estim_params_ , oo_ , options_ , bayestopt_ ) ;
iter = iter + 1 ;
end
if iter == 1000
fprintf ( ' \nprior_posterior_statistics: unable to get a valid draw in 1000 tries. Giving up.' )
end
2010-02-10 18:51:07 +01:00
else
deep = x ( b , : ) ;
2011-02-25 09:27:33 +01:00
if strcmpi ( type , ' posterior' )
logpo = logpost ( b ) ;
else
2017-10-03 15:32:17 +02:00
logpo = evaluate_posterior_kernel ( deep ' , M_ , estim_params_ , oo_ , options_ , bayestopt_ ) ;
2011-02-25 09:27:33 +01:00
end
2010-02-10 18:51:07 +01:00
end
2012-06-07 14:15:59 +02:00
M_ = set_all_parameters ( deep , estim_params_ , M_ ) ;
2011-02-02 14:13:11 +01:00
2010-02-10 18:51:07 +01:00
if run_smoother
2021-06-22 14:28:46 +02:00
[ dr , info , M_ , oo_ ] = compute_decision_rules ( M_ , opts_local , oo_ ) ;
2021-05-27 16:13:58 +02:00
if ismember ( info ( 1 ) , [ 3 , 4 ] )
opts_local . qz_criterium = 1 + ( opts_local . qz_criterium - 1 ) * 10 ; %loosen tolerance, useful for big models where BK conditions were tested on restricted state space
2021-06-22 14:28:46 +02:00
[ dr , info , M_ , oo_ ] = compute_decision_rules ( M_ , opts_local , oo_ ) ;
2021-05-27 16:13:58 +02:00
end
if info ( 1 )
message = get_error_message ( info , opts_local ) ;
fprintf ( ' \nprior_posterior_statistics: One of the draws failed with the error:\n%s\n' , message )
fprintf ( ' prior_posterior_statistics: This should not happen. Please contact the developers.\n' , message )
end
2021-05-17 21:24:25 +02:00
if options_ . occbin . smoother . status
2022-01-06 14:57:57 +01:00
opts_local . occbin . simul . waitbar = 0 ;
opts_local . occbin . smoother . waitbar = 0 ;
2023-02-06 14:58:42 +01:00
opts_local . occbin . smoother . linear_smoother = false ; % speed-up
2021-05-17 21:24:25 +02:00
[ alphahat , etahat , epsilonhat , alphatilde , SteadyState , trend_coeff , aK , ~ , ~ , P , ~ , ~ , trend_addition , state_uncertainty , M_ , oo_ , bayestopt_ ] = ...
occbin . DSGE_smoother ( deep , gend , Y , data_index , missing_value , M_ , oo_ , opts_local , bayestopt_ , estim_params_ ) ;
2023-02-06 14:58:42 +01:00
if oo_ . occbin . smoother . error_flag ( 1 )
message = get_error_message ( oo_ . occbin . smoother . error_flag , opts_local ) ;
fprintf ( ' \nprior_posterior_statistics: One of the draws failed with the error:\n%s\n' , message )
end
2021-05-17 21:24:25 +02:00
else
[ alphahat , etahat , epsilonhat , alphatilde , SteadyState , trend_coeff , aK , ~ , ~ , P , ~ , ~ , trend_addition , state_uncertainty , M_ , oo_ , bayestopt_ ] = ...
DsgeSmoother ( deep , gend , Y , data_index , missing_value , M_ , oo_ , opts_local , bayestopt_ , estim_params_ ) ;
end
2017-05-16 15:10:20 +02:00
2016-06-03 17:09:09 +02:00
stock_trend_coeff ( options_ . varobs_id , irun ( 9 ) ) = trend_coeff ;
stock_smoothed_trend ( IdObs , : , irun ( 11 ) ) = trend_addition ;
2015-03-11 18:15:28 +01:00
if options_ . loglinear %reads values from smoother results, which are in dr-order and put them into declaration order
2016-06-03 17:09:09 +02:00
constant_part = repmat ( log ( SteadyState ( dr . order_var ) ) , 1 , gend ) ;
2010-02-10 18:51:07 +01:00
stock_smooth ( dr . order_var , : , irun ( 1 ) ) = alphahat ( 1 : endo_nbr , : ) + ...
2016-06-03 17:09:09 +02:00
constant_part ;
2010-02-10 18:51:07 +01:00
stock_update ( dr . order_var , : , irun ( 1 ) ) = alphatilde ( 1 : endo_nbr , : ) + ...
2016-06-03 17:09:09 +02:00
constant_part ;
2010-02-10 18:51:07 +01:00
else
2016-06-03 17:09:09 +02:00
constant_part = repmat ( SteadyState ( dr . order_var ) , 1 , gend ) ;
2010-02-10 18:51:07 +01:00
stock_smooth ( dr . order_var , : , irun ( 1 ) ) = alphahat ( 1 : endo_nbr , : ) + ...
2016-06-03 17:09:09 +02:00
constant_part ;
2010-02-10 18:51:07 +01:00
stock_update ( dr . order_var , : , irun ( 1 ) ) = alphatilde ( 1 : endo_nbr , : ) + ...
2016-06-03 17:09:09 +02:00
constant_part ;
2017-05-16 15:10:20 +02:00
end
2016-06-03 17:09:09 +02:00
stock_smoothed_constant ( dr . order_var , : , irun ( 10 ) ) = constant_part ;
2015-03-11 18:15:28 +01:00
%% Compute constant for observables
if options_ . prefilter == 1 %as mean is taken after log transformation, no distinction is needed here
constant_part = repmat ( mean_varobs ' , 1 , gend ) ;
2019-03-19 14:26:16 +01:00
elseif options_ . prefilter == 0 && options_ . loglinear %logged steady state must be used
2015-03-11 18:15:28 +01:00
constant_part = repmat ( log ( SteadyState ( IdObs ) ) , 1 , gend ) ;
2019-03-19 14:26:16 +01:00
elseif options_ . prefilter == 0 && ~ options_ . loglinear %unlogged steady state must be used
2015-03-11 18:15:28 +01:00
constant_part = repmat ( SteadyState ( IdObs ) , 1 , gend ) ;
end
2015-03-01 15:40:47 +01:00
%add trend to observables
2015-03-11 18:15:28 +01:00
if options_ . prefilter
%do correction for prefiltering for observed variables
if options_ . loglinear
mean_correction = - repmat ( log ( SteadyState ( IdObs ) ) , 1 , gend ) + constant_part ;
else
mean_correction = - repmat ( SteadyState ( IdObs ) , 1 , gend ) + constant_part ;
end
2016-06-03 17:09:09 +02:00
stock_smoothed_constant ( IdObs , : , irun ( 10 ) ) = stock_smoothed_constant ( IdObs , : , irun ( 10 ) ) + mean_correction ;
2015-03-11 18:15:28 +01:00
%smoothed variables are E_T(y_t) so no trend shift is required
stock_smooth ( IdObs , : , irun ( 1 ) ) = stock_smooth ( IdObs , : , irun ( 1 ) ) + trend_addition + mean_correction ;
%updated variables are E_t(y_t) so no trend shift is required
2017-05-16 15:10:20 +02:00
stock_update ( IdObs , : , irun ( 1 ) ) = stock_update ( IdObs , : , irun ( 1 ) ) + trend_addition + mean_correction ;
2015-03-11 18:15:28 +01:00
else
stock_smooth ( IdObs , : , irun ( 1 ) ) = stock_smooth ( IdObs , : , irun ( 1 ) ) + trend_addition ;
2017-05-16 15:10:20 +02:00
stock_update ( IdObs , : , irun ( 1 ) ) = stock_update ( IdObs , : , irun ( 1 ) ) + trend_addition ;
2015-03-11 18:15:28 +01:00
end
2010-02-10 18:51:07 +01:00
stock_innov ( : , : , irun ( 2 ) ) = etahat ;
if nvn
stock_error ( : , : , irun ( 3 ) ) = epsilonhat ;
end
if naK
2015-03-11 18:15:28 +01:00
%filtered variable E_t(y_t+k) requires to shift trend by k periods
2016-06-07 11:57:11 +02:00
%write variables into declaration order
if options_ . loglinear %reads values from smoother results, which are in dr-order and put them into declaration order
2016-09-29 12:31:59 +02:00
constant_part = repmat ( log ( SteadyState ( dr . order_var ) ) ' , [ length ( options_ . filter_step_ahead ) , 1 , gend + max ( options_ . filter_step_ahead ) ] ) ;
2016-06-07 11:57:11 +02:00
else
constant_part = repmat ( SteadyState ( dr . order_var ) ' , [ length ( options_ . filter_step_ahead ) , 1 , gend + max ( options_ . filter_step_ahead ) ] ) ;
2017-05-16 15:10:20 +02:00
end
2017-07-23 23:21:36 +02:00
stock_filter_step_ahead ( : , dr . order_var , : , irun ( 4 ) ) = aK ( options_ . filter_step_ahead , 1 : endo_nbr , : ) + constant_part ;
2016-06-07 11:57:11 +02:00
%now add trend to observables
2015-03-11 18:15:28 +01:00
for ii = 1 : length ( options_ . filter_step_ahead )
2016-06-07 11:57:11 +02:00
if options_ . prefilter
2017-07-23 23:21:36 +02:00
zdim = size ( stock_filter_step_ahead ( ii , IdObs , : , irun ( 4 ) ) ) ;
squeezed = reshape ( stock_filter_step_ahead ( ii , IdObs , : , irun ( 4 ) ) , [ zdim ( 2 : end ) 1 ] ) ;
stock_filter_step_ahead ( ii , IdObs , : , irun ( 4 ) ) = squeezed ...
+ repmat ( mean_correction ( : , 1 ) , 1 , gend + max ( options_ . filter_step_ahead ) ) ... %constant correction
2016-06-07 11:57:11 +02:00
+ [ trend_addition repmat ( trend_addition ( : , end ) , 1 , max ( options_ . filter_step_ahead ) ) + trend_coeff * [ 1 : max ( options_ . filter_step_ahead ) ] ] ; %trend
else
2017-07-23 23:21:36 +02:00
zdim = size ( stock_filter_step_ahead ( ii , IdObs , : , irun ( 4 ) ) ) ;
squeezed = reshape ( stock_filter_step_ahead ( ii , IdObs , : , irun ( 4 ) ) , [ zdim ( 2 : end ) 1 ] ) ;
stock_filter_step_ahead ( ii , IdObs , : , irun ( 4 ) ) = squeezed ...
2016-06-07 11:57:11 +02:00
+ [ trend_addition repmat ( trend_addition ( : , end ) , 1 , max ( options_ . filter_step_ahead ) ) + trend_coeff * [ 1 : max ( options_ . filter_step_ahead ) ] ] ; %trend
end
2015-03-11 18:15:28 +01:00
end
2010-02-10 18:51:07 +01:00
end
if horizon
yyyy = alphahat ( iendo , i_last_obs ) ;
yf = forcst2a ( yyyy , dr , zeros ( horizon , exo_nbr ) ) ;
2017-05-16 15:10:20 +02:00
if options_ . prefilter
2015-03-01 15:12:34 +01:00
% add mean
2015-03-02 10:45:12 +01:00
yf ( : , IdObs ) = yf ( : , IdObs ) + repmat ( mean_varobs , ...
2011-02-04 17:17:48 +01:00
horizon + maxlag , 1 ) ;
2015-03-01 15:40:47 +01:00
% add trend, taking into account that last point of sample is still included in forecasts and only cut off later
yf ( : , IdObs ) = yf ( : , IdObs ) + ( ( options_ . first_obs - 1 ) + gend + [ 1 - maxlag : horizon ] ' ) * trend_coeff ' - ...
2017-05-16 15:10:20 +02:00
repmat ( mean ( trend_coeff * [ options_ . first_obs : options_ . first_obs + gend - 1 ] , 2 ) ' , length ( 1 - maxlag : horizon ) , 1 ) ; %center trend
2015-03-01 15:12:34 +01:00
else
2015-03-01 15:40:47 +01:00
% add trend, taking into account that last point of sample is still included in forecasts and only cut off later
2017-05-16 15:10:20 +02:00
yf ( : , IdObs ) = yf ( : , IdObs ) + ( ( options_ . first_obs - 1 ) + gend + [ 1 - maxlag : horizon ] ' ) * trend_coeff ' ;
2010-02-10 18:51:07 +01:00
end
2014-01-30 13:00:25 +01:00
if options_ . loglinear
2010-02-10 18:51:07 +01:00
yf = yf + repmat ( log ( SteadyState ' ) , horizon + maxlag , 1 ) ;
else
yf = yf + repmat ( SteadyState ' , horizon + maxlag , 1 ) ;
end
yf1 = forcst2 ( yyyy , horizon , dr , 1 ) ;
if options_ . prefilter == 1
2015-03-01 15:12:34 +01:00
% add mean
2010-02-10 18:51:07 +01:00
yf1 ( : , IdObs , : ) = yf1 ( : , IdObs , : ) + ...
2015-03-02 10:45:12 +01:00
repmat ( mean_varobs , [ horizon + maxlag , 1 , 1 ] ) ;
2015-03-01 15:40:47 +01:00
% add trend, taking into account that last point of sample is still included in forecasts and only cut off later
yf1 ( : , IdObs ) = yf1 ( : , IdObs ) + ( ( options_ . first_obs - 1 ) + gend + [ 1 - maxlag : horizon ] ' ) * trend_coeff ' - ...
2017-05-16 15:10:20 +02:00
repmat ( mean ( trend_coeff * [ options_ . first_obs : options_ . first_obs + gend - 1 ] , 2 ) ' , length ( 1 - maxlag : horizon ) , 1 ) ; %center trend
2015-03-01 15:12:34 +01:00
else
2017-05-16 15:10:20 +02:00
% add trend, taking into account that last point of sample is still included in forecasts and only cut off later
yf1 ( : , IdObs , : ) = yf1 ( : , IdObs , : ) + repmat ( ( ( options_ . first_obs - 1 ) + gend + [ 1 - maxlag : horizon ] ' ) * ...
2015-03-01 15:12:34 +01:00
trend_coeff ' , [ 1 , 1 , 1 ] ) ;
2010-02-10 18:51:07 +01:00
end
2014-01-30 13:00:25 +01:00
if options_ . loglinear
2010-02-10 18:51:07 +01:00
yf1 = yf1 + repmat ( log ( SteadyState ' ) , [ horizon + maxlag , 1 , 1 ] ) ;
else
yf1 = yf1 + repmat ( SteadyState ' , [ horizon + maxlag , 1 , 1 ] ) ;
end
2011-02-02 14:13:11 +01:00
2014-10-30 20:12:45 +01:00
stock_forcst_mean ( : , : , irun ( 6 ) ) = yf ( maxlag + 1 : end , : ) ' ;
stock_forcst_point ( : , : , irun ( 7 ) ) = yf1 ( maxlag + 1 : end , : ) ' ;
2016-06-18 14:33:11 +02:00
if ~ isequal ( M_ . H , 0 )
ME_shocks = zeros ( length ( varobs ) , horizon ) ;
i_exo_var = setdiff ( [ 1 : length ( varobs ) ] , find ( diag ( M_ . H ) == 0 ) ) ;
nxs = length ( i_exo_var ) ;
chol_H = chol ( M_ . H ( i_exo_var , i_exo_var ) ) ;
if ~ isempty ( M_ . H )
ME_shocks ( i_exo_var , : ) = chol_H * randn ( nxs , horizon ) ;
end
stock_forcst_point_ME ( : , : , irun ( 12 ) ) = yf1 ( maxlag + 1 : end , IdObs ) ' + ME_shocks ;
end
2010-02-10 18:51:07 +01:00
end
2016-06-02 18:28:36 +02:00
if filter_covariance
stock_filter_covariance ( dr . order_var , dr . order_var , : , irun ( 8 ) ) = P ;
end
2016-11-02 15:55:23 +01:00
if smoothed_state_uncertainty
stock_smoothed_uncert ( dr . order_var , dr . order_var , : , irun ( 13 ) ) = state_uncertainty ;
end
2013-02-08 14:30:40 +01:00
else
2021-06-22 14:28:46 +02:00
[ T , R , SteadyState , info , M_ , oo_ ] = dynare_resolve ( M_ , options_ , oo_ ) ;
2010-02-10 18:51:07 +01:00
end
stock_param ( irun ( 5 ) , : ) = deep ;
stock_logpo ( irun ( 5 ) , 1 ) = logpo ;
stock_ys ( irun ( 5 ) , : ) = SteadyState ' ;
2011-02-02 14:13:11 +01:00
2016-06-02 18:28:36 +02:00
2016-11-02 15:55:23 +01:00
irun = irun + ones ( 13 , 1 ) ;
2011-02-02 14:13:11 +01:00
2017-05-16 12:42:01 +02:00
if run_smoother && ( irun ( 1 ) > MAX_nsmoo || b == B )
2010-02-10 18:51:07 +01:00
stock = stock_smooth ( : , : , 1 : irun ( 1 ) - 1 ) ;
ifil ( 1 ) = ifil ( 1 ) + 1 ;
save ( [ DirectoryName ' /' M_ . fname ' _smooth' int2str ( ifil ( 1 ) ) ' .mat' ] , ' stock' ) ;
stock = stock_update ( : , : , 1 : irun ( 1 ) - 1 ) ;
save ( [ DirectoryName ' /' M_ . fname ' _update' int2str ( ifil ( 1 ) ) ' .mat' ] , ' stock' ) ;
2017-05-16 12:42:01 +02:00
if RemoteFlag == 1
2011-02-02 14:13:11 +01:00
OutputFileName_smooth = [ OutputFileName_smooth ; { [ DirectoryName filesep ] , [ M_ . fname ' _smooth' int2str ( ifil ( 1 ) ) ' .mat' ] } ] ;
OutputFileName_update = [ OutputFileName_update ; { [ DirectoryName filesep ] , [ M_ . fname ' _update' int2str ( ifil ( 1 ) ) ' .mat' ] } ] ;
2010-02-10 18:51:07 +01:00
end
irun ( 1 ) = 1 ;
end
2011-02-02 14:13:11 +01:00
2013-02-08 14:30:40 +01:00
if run_smoother && ( irun ( 2 ) > MAX_ninno || b == B )
2010-02-10 18:51:07 +01:00
stock = stock_innov ( : , : , 1 : irun ( 2 ) - 1 ) ;
ifil ( 2 ) = ifil ( 2 ) + 1 ;
save ( [ DirectoryName ' /' M_ . fname ' _inno' int2str ( ifil ( 2 ) ) ' .mat' ] , ' stock' ) ;
2017-05-16 12:42:01 +02:00
if RemoteFlag == 1
2011-02-02 14:13:11 +01:00
OutputFileName_inno = [ OutputFileName_inno ; { [ DirectoryName filesep ] , [ M_ . fname ' _inno' int2str ( ifil ( 2 ) ) ' .mat' ] } ] ;
2010-02-10 18:51:07 +01:00
end
irun ( 2 ) = 1 ;
end
2011-02-02 14:13:11 +01:00
2013-02-08 14:30:40 +01:00
if run_smoother && nvn && ( irun ( 3 ) > MAX_nerro || b == B )
2010-02-10 18:51:07 +01:00
stock = stock_error ( : , : , 1 : irun ( 3 ) - 1 ) ;
ifil ( 3 ) = ifil ( 3 ) + 1 ;
save ( [ DirectoryName ' /' M_ . fname ' _error' int2str ( ifil ( 3 ) ) ' .mat' ] , ' stock' ) ;
2017-05-16 12:42:01 +02:00
if RemoteFlag == 1
2011-02-02 14:13:11 +01:00
OutputFileName_error = [ OutputFileName_error ; { [ DirectoryName filesep ] , [ M_ . fname ' _error' int2str ( ifil ( 3 ) ) ' .mat' ] } ] ;
2010-02-10 18:51:07 +01:00
end
irun ( 3 ) = 1 ;
end
2011-02-02 14:13:11 +01:00
2013-02-08 14:30:40 +01:00
if run_smoother && naK && ( irun ( 4 ) > MAX_naK || b == B )
2010-02-10 18:51:07 +01:00
stock = stock_filter_step_ahead ( : , : , : , 1 : irun ( 4 ) - 1 ) ;
ifil ( 4 ) = ifil ( 4 ) + 1 ;
save ( [ DirectoryName ' /' M_ . fname ' _filter_step_ahead' int2str ( ifil ( 4 ) ) ' .mat' ] , ' stock' ) ;
2017-05-16 12:42:01 +02:00
if RemoteFlag == 1
2011-02-02 14:13:11 +01:00
OutputFileName_filter_step_ahead = [ OutputFileName_filter_step_ahead ; { [ DirectoryName filesep ] , [ M_ . fname ' _filter_step_ahead' int2str ( ifil ( 4 ) ) ' .mat' ] } ] ;
2010-02-10 18:51:07 +01:00
end
irun ( 4 ) = 1 ;
end
2011-02-02 14:13:11 +01:00
2010-02-10 18:51:07 +01:00
if irun ( 5 ) > MAX_nruns || b == B
stock = stock_param ( 1 : irun ( 5 ) - 1 , : ) ;
2015-04-14 12:10:46 +02:00
stock_logpo = stock_logpo ( 1 : irun ( 5 ) - 1 ) ;
stock_ys = stock_ys ( 1 : irun ( 5 ) - 1 , : ) ;
2010-02-10 18:51:07 +01:00
ifil ( 5 ) = ifil ( 5 ) + 1 ;
save ( [ DirectoryName ' /' M_ . fname ' _param' int2str ( ifil ( 5 ) ) ' .mat' ] , ' stock' , ' stock_logpo' , ' stock_ys' ) ;
2017-05-16 12:42:01 +02:00
if RemoteFlag == 1
2011-02-02 14:13:11 +01:00
OutputFileName_param = [ OutputFileName_param ; { [ DirectoryName filesep ] , [ M_ . fname ' _param' int2str ( ifil ( 5 ) ) ' .mat' ] } ] ;
2010-02-10 18:51:07 +01:00
end
irun ( 5 ) = 1 ;
end
2011-02-02 14:13:11 +01:00
2013-02-08 14:30:40 +01:00
if run_smoother && horizon && ( irun ( 6 ) > MAX_nforc1 || b == B )
2010-02-10 18:51:07 +01:00
stock = stock_forcst_mean ( : , : , 1 : irun ( 6 ) - 1 ) ;
ifil ( 6 ) = ifil ( 6 ) + 1 ;
save ( [ DirectoryName ' /' M_ . fname ' _forc_mean' int2str ( ifil ( 6 ) ) ' .mat' ] , ' stock' ) ;
2017-05-16 12:42:01 +02:00
if RemoteFlag == 1
2011-02-02 14:13:11 +01:00
OutputFileName_forc_mean = [ OutputFileName_forc_mean ; { [ DirectoryName filesep ] , [ M_ . fname ' _forc_mean' int2str ( ifil ( 6 ) ) ' .mat' ] } ] ;
2010-02-10 18:51:07 +01:00
end
irun ( 6 ) = 1 ;
end
2011-02-02 14:13:11 +01:00
2013-02-08 14:30:40 +01:00
if run_smoother && horizon && ( irun ( 7 ) > MAX_nforc2 || b == B )
2010-02-10 18:51:07 +01:00
stock = stock_forcst_point ( : , : , 1 : irun ( 7 ) - 1 ) ;
ifil ( 7 ) = ifil ( 7 ) + 1 ;
save ( [ DirectoryName ' /' M_ . fname ' _forc_point' int2str ( ifil ( 7 ) ) ' .mat' ] , ' stock' ) ;
2017-05-16 12:42:01 +02:00
if RemoteFlag == 1
2011-02-02 14:13:11 +01:00
OutputFileName_forc_point = [ OutputFileName_forc_point ; { [ DirectoryName filesep ] , [ M_ . fname ' _forc_point' int2str ( ifil ( 7 ) ) ' .mat' ] } ] ;
2010-02-10 18:51:07 +01:00
end
irun ( 7 ) = 1 ;
end
2017-05-16 15:10:20 +02:00
2016-06-02 18:28:36 +02:00
if run_smoother && filter_covariance && ( irun ( 8 ) > MAX_filter_covariance || b == B )
stock = stock_filter_covariance ( : , : , : , 1 : irun ( 8 ) - 1 ) ;
ifil ( 8 ) = ifil ( 8 ) + 1 ;
save ( [ DirectoryName ' /' M_ . fname ' _filter_covar' int2str ( ifil ( 8 ) ) ' .mat' ] , ' stock' ) ;
2017-05-16 12:42:01 +02:00
if RemoteFlag == 1
2016-06-02 18:28:36 +02:00
OutputFileName_filter_covar = [ OutputFileName_filter_covar ; { [ DirectoryName filesep ] , [ M_ . fname ' _filter_covar' int2str ( ifil ( 8 ) ) ' .mat' ] } ] ;
end
irun ( 8 ) = 1 ;
end
2017-05-16 15:10:20 +02:00
2016-06-03 17:09:09 +02:00
irun_index = 9 ;
if run_smoother && ( irun ( irun_index ) > MAX_n_trend_coeff || b == B )
stock = stock_trend_coeff ( : , 1 : irun ( irun_index ) - 1 ) ;
ifil ( irun_index ) = ifil ( irun_index ) + 1 ;
save ( [ DirectoryName ' /' M_ . fname ' _trend_coeff' int2str ( ifil ( irun_index ) ) ' .mat' ] , ' stock' ) ;
2017-05-16 12:42:01 +02:00
if RemoteFlag == 1
2016-06-03 17:09:09 +02:00
OutputFileName_trend_coeff = [ OutputFileName_trend_coeff ; { [ DirectoryName filesep ] , [ M_ . fname ' _trend_coeff' int2str ( ifil ( irun_index ) ) ' .mat' ] } ] ;
end
irun ( irun_index ) = 1 ;
end
2017-05-16 15:10:20 +02:00
2016-06-03 17:09:09 +02:00
irun_index = 10 ;
if run_smoother && ( irun ( irun_index ) > MAX_n_smoothed_constant || b == B )
stock = stock_smoothed_constant ( : , : , 1 : irun ( irun_index ) - 1 ) ;
ifil ( irun_index ) = ifil ( irun_index ) + 1 ;
save ( [ DirectoryName ' /' M_ . fname ' _smoothed_constant' int2str ( ifil ( irun_index ) ) ' .mat' ] , ' stock' ) ;
2017-05-16 12:42:01 +02:00
if RemoteFlag == 1
2016-06-03 17:09:09 +02:00
OutputFileName_smoothed_constant = [ OutputFileName_smoothed_constant ; { [ DirectoryName filesep ] , [ M_ . fname ' _smoothed_constant' int2str ( ifil ( irun_index ) ) ' .mat' ] } ] ;
end
irun ( irun_index ) = 1 ;
end
irun_index = 11 ;
if run_smoother && ( irun ( irun_index ) > MAX_n_smoothed_trend || b == B )
stock = stock_smoothed_trend ( : , : , 1 : irun ( irun_index ) - 1 ) ;
ifil ( irun_index ) = ifil ( irun_index ) + 1 ;
save ( [ DirectoryName ' /' M_ . fname ' _smoothed_trend' int2str ( ifil ( irun_index ) ) ' .mat' ] , ' stock' ) ;
2017-05-16 12:42:01 +02:00
if RemoteFlag == 1
2016-06-03 17:09:09 +02:00
OutputFileName_smoothed_trend = [ OutputFileName_smoothed_trend ; { [ DirectoryName filesep ] , [ M_ . fname ' _smoothed_trend' int2str ( ifil ( irun_index ) ) ' .mat' ] } ] ;
end
irun ( irun_index ) = 1 ;
end
2017-05-16 15:10:20 +02:00
2016-06-18 14:33:11 +02:00
irun_index = 12 ;
if run_smoother && horizon && ~ isequal ( M_ . H , 0 ) && ( irun ( irun_index ) > MAX_nforc_ME || b == B )
stock = stock_forcst_point_ME ( : , : , 1 : irun ( irun_index ) - 1 ) ;
ifil ( irun_index ) = ifil ( irun_index ) + 1 ;
save ( [ DirectoryName ' /' M_ . fname ' _forc_point_ME' int2str ( ifil ( irun_index ) ) ' .mat' ] , ' stock' ) ;
2017-05-16 12:42:01 +02:00
if RemoteFlag == 1
2016-06-18 14:33:11 +02:00
OutputFileName_forc_point_ME = [ OutputFileName_forc_point_ME ; { [ DirectoryName filesep ] , [ M_ . fname ' _forc_point_ME' int2str ( ifil ( irun_index ) ) ' .mat' ] } ] ;
end
irun ( irun_index ) = 1 ;
end
2017-05-16 15:10:20 +02:00
2016-11-02 15:55:23 +01:00
irun_index = 13 ;
if run_smoother && smoothed_state_uncertainty && ( irun ( irun_index ) > MAX_n_smoothed_state_uncertainty || b == B )
stock = stock_smoothed_uncert ( : , : , : , 1 : irun ( irun_index ) - 1 ) ;
ifil ( irun_index ) = ifil ( irun_index ) + 1 ;
save ( [ DirectoryName ' /' M_ . fname ' _state_uncert' int2str ( ifil ( irun_index ) ) ' .mat' ] , ' stock' ) ;
2017-05-16 12:42:01 +02:00
if RemoteFlag == 1
2016-11-02 15:55:23 +01:00
OutputFileName_state_uncert = [ OutputFileName_state_uncert ; { [ DirectoryName filesep ] , [ M_ . fname ' _state_uncert' int2str ( ifil ( irun_index ) ) ' .mat' ] } ] ;
end
irun ( irun_index ) = 1 ;
end
2011-12-13 18:32:57 +01:00
dyn_waitbar ( ( b - fpar + 1 ) / ( B - fpar + 1 ) , h ) ;
2010-02-10 18:51:07 +01:00
end
myoutput . ifil = ifil ;
2017-05-16 12:42:01 +02:00
if RemoteFlag == 1
2011-02-02 14:13:11 +01:00
myoutput . OutputFileName = [ OutputFileName_smooth ;
2011-02-04 17:17:48 +01:00
OutputFileName_update ;
OutputFileName_inno ;
OutputFileName_error ;
OutputFileName_filter_step_ahead ;
OutputFileName_param ;
OutputFileName_forc_mean ;
2016-06-07 11:57:11 +02:00
OutputFileName_forc_point ;
2016-06-18 14:33:11 +02:00
OutputFileName_forc_point_ME ;
2016-06-07 11:57:11 +02:00
OutputFileName_filter_covar ;
OutputFileName_trend_coeff ;
OutputFileName_smoothed_trend ;
2016-11-02 15:55:23 +01:00
OutputFileName_smoothed_constant ;
OutputFileName_state_uncert ] ;
2010-02-10 18:51:07 +01:00
end
2010-10-12 08:07:57 +02:00
2011-12-13 18:32:57 +01:00
dyn_waitbar_close ( h ) ;