2013-12-09 12:32:20 +01:00
|
|
|
function ts = extended_path(initial_conditions,sample_size)
|
2011-12-05 10:58:39 +01:00
|
|
|
% Stochastic simulation of a non linear DSGE model using the Extended Path method (Fair and Taylor 1983). A time
|
2012-01-27 18:24:24 +01:00
|
|
|
% series of size T is obtained by solving T perfect foresight models.
|
|
|
|
%
|
2011-12-05 10:58:39 +01:00
|
|
|
% INPUTS
|
|
|
|
% o initial_conditions [double] m*nlags array, where m is the number of endogenous variables in the model and
|
|
|
|
% nlags is the maximum number of lags.
|
|
|
|
% o sample_size [integer] scalar, size of the sample to be simulated.
|
2012-01-27 18:24:24 +01:00
|
|
|
%
|
2011-12-05 10:58:39 +01:00
|
|
|
% OUTPUTS
|
|
|
|
% o time_series [double] m*sample_size array, the simulations.
|
2012-01-27 18:24:24 +01:00
|
|
|
%
|
2011-12-05 10:58:39 +01:00
|
|
|
% ALGORITHM
|
2012-01-27 18:24:24 +01:00
|
|
|
%
|
2011-12-05 10:58:39 +01:00
|
|
|
% SPECIAL REQUIREMENTS
|
|
|
|
|
2013-06-12 16:42:09 +02:00
|
|
|
% Copyright (C) 2009-2013 Dynare Team
|
2011-12-05 10:58:39 +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.
|
|
|
|
%
|
|
|
|
% You should have received a copy of the GNU General Public License
|
|
|
|
% along with Dynare. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
global M_ options_ oo_
|
2012-01-21 14:13:31 +01:00
|
|
|
|
2011-12-26 14:26:32 +01:00
|
|
|
options_.verbosity = options_.ep.verbosity;
|
2012-01-21 14:13:31 +01:00
|
|
|
verbosity = options_.ep.verbosity+options_.ep.debug;
|
2011-12-12 12:40:46 +01:00
|
|
|
|
2014-01-17 21:49:04 +01:00
|
|
|
% Set maximum number of iterations for the deterministic solver.
|
|
|
|
options_.simul.maxit = options_.ep.maxit;
|
|
|
|
|
2012-02-03 12:39:48 +01:00
|
|
|
% Prepare a structure needed by the matlab implementation of the perfect foresight model solver
|
2014-01-17 21:49:04 +01:00
|
|
|
pfm = setup_stochastic_perfect_foresight_model_solver(M_,options_,oo_);
|
2011-12-05 10:58:39 +01:00
|
|
|
|
2012-02-07 16:30:36 +01:00
|
|
|
exo_nbr = M_.exo_nbr;
|
|
|
|
periods = options_.periods;
|
|
|
|
ep = options_.ep;
|
|
|
|
steady_state = oo_.steady_state;
|
|
|
|
dynatol = options_.dynatol;
|
|
|
|
|
2011-12-05 10:58:39 +01:00
|
|
|
% Set default initial conditions.
|
|
|
|
if isempty(initial_conditions)
|
2013-12-11 11:27:19 +01:00
|
|
|
if isempty(M_.endo_histval)
|
|
|
|
initial_conditions = oo_.steady_state;
|
|
|
|
else
|
|
|
|
initial_conditions = M_.endo_histval;
|
|
|
|
end
|
2011-12-05 10:58:39 +01:00
|
|
|
end
|
|
|
|
|
|
|
|
|
|
|
|
% Set the number of periods for the perfect foresight model
|
2012-02-07 16:30:36 +01:00
|
|
|
periods = options_.ep.periods;
|
2014-05-12 09:40:12 +02:00
|
|
|
pfm.periods = periods;
|
2012-02-03 12:39:48 +01:00
|
|
|
pfm.i_upd = pfm.ny+(1:pfm.periods*pfm.ny);
|
2011-12-05 10:58:39 +01:00
|
|
|
|
2012-03-02 11:33:29 +01:00
|
|
|
% keep a copy of pfm.i_upd
|
|
|
|
i_upd = pfm.i_upd;
|
|
|
|
|
2011-12-09 18:03:30 +01:00
|
|
|
% Set the algorithm for the perfect foresight solver
|
|
|
|
options_.stack_solve_algo = options_.ep.stack_solve_algo;
|
|
|
|
|
2012-02-04 16:19:15 +01:00
|
|
|
% Set check_stability flag
|
|
|
|
do_not_check_stability_flag = ~options_.ep.check_stability;
|
|
|
|
|
2011-12-05 10:58:39 +01:00
|
|
|
% Compute the first order reduced form if needed.
|
|
|
|
%
|
|
|
|
% REMARK. It is assumed that the user did run the same mod file with stoch_simul(order=1) and save
|
|
|
|
% all the globals in a mat file called linear_reduced_form.mat;
|
2012-01-27 18:24:24 +01:00
|
|
|
|
2012-02-07 16:30:36 +01:00
|
|
|
dr = struct();
|
2011-12-05 10:58:39 +01:00
|
|
|
if options_.ep.init
|
2012-01-27 18:24:24 +01:00
|
|
|
options_.order = 1;
|
|
|
|
[dr,Info,M_,options_,oo_] = resol(1,M_,options_,oo_);
|
2011-12-05 10:58:39 +01:00
|
|
|
end
|
|
|
|
|
|
|
|
% Do not use a minimal number of perdiods for the perfect foresight solver (with bytecode and blocks)
|
2012-01-27 18:24:24 +01:00
|
|
|
options_.minimal_solving_period = 100;%options_.ep.periods;
|
2011-12-05 10:58:39 +01:00
|
|
|
|
|
|
|
% Initialize the exogenous variables.
|
|
|
|
make_ex_;
|
|
|
|
|
|
|
|
% Initialize the endogenous variables.
|
|
|
|
make_y_;
|
|
|
|
|
|
|
|
% Initialize the output array.
|
2011-12-19 18:04:54 +01:00
|
|
|
time_series = zeros(M_.endo_nbr,sample_size);
|
2011-12-05 10:58:39 +01:00
|
|
|
|
|
|
|
% Set the covariance matrix of the structural innovations.
|
2012-01-27 18:24:24 +01:00
|
|
|
variances = diag(M_.Sigma_e);
|
|
|
|
positive_var_indx = find(variances>0);
|
2012-01-20 16:39:12 +01:00
|
|
|
effective_number_of_shocks = length(positive_var_indx);
|
|
|
|
stdd = sqrt(variances(positive_var_indx));
|
2012-01-27 18:24:24 +01:00
|
|
|
covariance_matrix = M_.Sigma_e(positive_var_indx,positive_var_indx);
|
|
|
|
covariance_matrix_upper_cholesky = chol(covariance_matrix);
|
2011-12-05 10:58:39 +01:00
|
|
|
|
2012-02-10 12:56:52 +01:00
|
|
|
% (re)Set exo_nbr
|
2012-03-03 15:32:44 +01:00
|
|
|
%exo_nbr = effective_number_of_shocks;
|
2012-02-10 12:56:52 +01:00
|
|
|
|
2011-12-12 14:46:13 +01:00
|
|
|
% Set seed.
|
|
|
|
if options_.ep.set_dynare_seed_to_default
|
|
|
|
set_dynare_seed('default');
|
|
|
|
end
|
|
|
|
|
2012-02-04 16:26:22 +01:00
|
|
|
% Set bytecode flag
|
|
|
|
bytecode_flag = options_.ep.use_bytecode;
|
|
|
|
|
2011-12-05 10:58:39 +01:00
|
|
|
% Simulate shocks.
|
|
|
|
switch options_.ep.innovation_distribution
|
|
|
|
case 'gaussian'
|
2013-12-17 14:14:51 +01:00
|
|
|
oo_.ep.shocks = transpose(transpose(covariance_matrix_upper_cholesky)*randn(effective_number_of_shocks,sample_size));
|
2011-12-05 10:58:39 +01:00
|
|
|
otherwise
|
|
|
|
error(['extended_path:: ' options_.ep.innovation_distribution ' distribution for the structural innovations is not (yet) implemented!'])
|
|
|
|
end
|
2012-01-27 18:24:24 +01:00
|
|
|
|
2011-12-05 10:58:39 +01:00
|
|
|
|
2011-12-12 14:48:26 +01:00
|
|
|
% Set waitbar (graphic or text mode)
|
2012-01-11 16:00:42 +01:00
|
|
|
hh = dyn_waitbar(0,'Please wait. Extended Path simulations...');
|
|
|
|
set(hh,'Name','EP simulations.');
|
2011-12-12 14:01:17 +01:00
|
|
|
|
2013-03-12 09:36:58 +01:00
|
|
|
% hybrid correction
|
|
|
|
pfm.hybrid_order = options_.ep.stochastic.hybrid_order;
|
|
|
|
if pfm.hybrid_order
|
|
|
|
oo_.dr = set_state_space(oo_.dr,M_,options_);
|
|
|
|
options = options_;
|
|
|
|
options.order = pfm.hybrid_order;
|
|
|
|
pfm.dr = resol(0,M_,options,oo_);
|
|
|
|
else
|
|
|
|
pfm.dr = [];
|
|
|
|
end
|
|
|
|
|
2014-05-12 09:40:12 +02:00
|
|
|
% number of nonzero derivatives
|
|
|
|
pfm.nnzA = M_.NNZDerivatives(1);
|
|
|
|
|
|
|
|
% setting up integration nodes if order > 0
|
|
|
|
if options_.ep.stochastic.order > 0
|
|
|
|
[nodes,weights,nnodes] = setup_integration_nodes(options_.ep,pfm);
|
|
|
|
pfm.nodes = nodes;
|
2014-06-16 14:25:23 +02:00
|
|
|
pfm.weights = weights;
|
2014-05-12 09:40:12 +02:00
|
|
|
pfm.nnodes = nnodes;
|
2014-06-16 14:25:23 +02:00
|
|
|
|
|
|
|
% compute number of blocks
|
|
|
|
[block_nbr,pfm.world_nbr] = get_block_world_nbr(options_.ep.stochastic.algo,nnodes,options_.ep.ut.k,options_.ep.periods);
|
|
|
|
else
|
|
|
|
block_nbr = options_.ep.periods
|
2014-05-12 09:40:12 +02:00
|
|
|
end
|
|
|
|
|
2014-06-16 14:25:23 +02:00
|
|
|
|
2014-05-12 17:16:09 +02:00
|
|
|
% set boundaries if mcp
|
|
|
|
[lb,ub,pfm.eq_index] = get_complementarity_conditions(M_);
|
|
|
|
options_.lmmcp.lb = repmat(lb,block_nbr,1);
|
|
|
|
options_.lmmcp.ub = repmat(ub,block_nbr,1);
|
|
|
|
pfm.block_nbr = block_nbr;
|
|
|
|
|
2014-06-16 14:25:23 +02:00
|
|
|
% storage for failed draws
|
|
|
|
oo_.ep.failures.periods = [];
|
|
|
|
oo_.ep.failures.previous_period = cell(0);
|
|
|
|
oo_.ep.failures.shocks = cell(0);
|
|
|
|
|
|
|
|
% Initializes some variables.
|
|
|
|
t = 0;
|
|
|
|
tsimul = 1;
|
2011-12-05 10:58:39 +01:00
|
|
|
% Main loop.
|
|
|
|
while (t<sample_size)
|
2011-12-12 14:01:17 +01:00
|
|
|
if ~mod(t,10)
|
2012-01-11 17:10:24 +01:00
|
|
|
dyn_waitbar(t/sample_size,hh,'Please wait. Extended Path simulations...');
|
2011-12-12 14:01:17 +01:00
|
|
|
end
|
2011-12-05 10:58:39 +01:00
|
|
|
% Set period index.
|
|
|
|
t = t+1;
|
|
|
|
shocks = oo_.ep.shocks(t,:);
|
|
|
|
% Put it in oo_.exo_simul (second line).
|
|
|
|
oo_.exo_simul(2,positive_var_indx) = shocks;
|
2012-02-10 12:56:52 +01:00
|
|
|
periods1 = periods;
|
|
|
|
exo_simul_1 = zeros(periods1+2,exo_nbr);
|
|
|
|
exo_simul_1(2,:) = oo_.exo_simul(2,:);
|
|
|
|
pfm1 = pfm;
|
|
|
|
info_convergence = 0;
|
|
|
|
if ep.init% Compute first order solution (Perturbation)...
|
|
|
|
ex = zeros(size(endo_simul_1,2),size(exo_simul_1,2));
|
|
|
|
ex(1:size(exo_simul_1,1),:) = exo_simul_1;
|
|
|
|
exo_simul_1 = ex;
|
|
|
|
initial_path = simult_(initial_conditions,dr,exo_simul_1(2:end,:),1);
|
|
|
|
endo_simul_1(:,1:end-1) = initial_path(:,1:end-1)*ep.init+endo_simul_1(:,1:end-1)*(1-ep.init);
|
|
|
|
else
|
2012-03-02 11:33:29 +01:00
|
|
|
if t==1
|
|
|
|
endo_simul_1 = repmat(steady_state,1,periods1+2);
|
|
|
|
end
|
2012-02-10 12:56:52 +01:00
|
|
|
end
|
|
|
|
% Solve a perfect foresight model.
|
|
|
|
increase_periods = 0;
|
2012-03-02 11:33:29 +01:00
|
|
|
% Keep a copy of endo_simul_1
|
2012-02-10 12:56:52 +01:00
|
|
|
endo_simul = endo_simul_1;
|
2014-03-16 18:49:07 +01:00
|
|
|
if verbosity
|
|
|
|
save ep_test_1 endo_simul_1 exo_simul_1
|
|
|
|
end
|
2012-02-10 12:56:52 +01:00
|
|
|
while 1
|
|
|
|
if ~increase_periods
|
2012-03-12 10:34:28 +01:00
|
|
|
if bytecode_flag && ~options_.ep.stochastic.order
|
|
|
|
[flag,tmp] = bytecode('dynamic',endo_simul_1,exo_simul_1);
|
2012-02-10 12:56:52 +01:00
|
|
|
else
|
|
|
|
flag = 1;
|
2012-01-20 16:39:12 +01:00
|
|
|
end
|
2012-02-10 12:56:52 +01:00
|
|
|
if flag
|
2012-03-08 12:49:17 +01:00
|
|
|
if options_.ep.stochastic.order == 0
|
|
|
|
[flag,tmp,err] = solve_perfect_foresight_model(endo_simul_1,exo_simul_1,pfm1);
|
2012-03-05 12:35:55 +01:00
|
|
|
else
|
2013-03-12 09:36:58 +01:00
|
|
|
switch(options_.ep.stochastic.algo)
|
|
|
|
case 0
|
|
|
|
[flag,tmp] = ...
|
|
|
|
solve_stochastic_perfect_foresight_model(endo_simul_1,exo_simul_1,pfm1,options_.ep.stochastic.quadrature.nodes,options_.ep.stochastic.order);
|
|
|
|
case 1
|
|
|
|
[flag,tmp] = ...
|
2014-05-12 09:40:12 +02:00
|
|
|
solve_stochastic_perfect_foresight_model_1(endo_simul_1,exo_simul_1,options_,pfm1,options_.ep.stochastic.order);
|
2013-03-12 09:36:58 +01:00
|
|
|
end
|
2012-03-05 12:35:55 +01:00
|
|
|
end
|
2012-02-10 12:56:52 +01:00
|
|
|
end
|
|
|
|
info_convergence = ~flag;
|
2011-12-05 10:58:39 +01:00
|
|
|
end
|
2012-02-10 12:56:52 +01:00
|
|
|
if verbosity
|
|
|
|
if info_convergence
|
|
|
|
if t<10
|
|
|
|
disp(['Time: ' int2str(t) '. Convergence of the perfect foresight model solver!'])
|
|
|
|
elseif t<100
|
|
|
|
disp(['Time: ' int2str(t) '. Convergence of the perfect foresight model solver!'])
|
|
|
|
elseif t<1000
|
|
|
|
disp(['Time: ' int2str(t) '. Convergence of the perfect foresight model solver!'])
|
2012-02-04 16:26:22 +01:00
|
|
|
else
|
2012-02-10 12:56:52 +01:00
|
|
|
disp(['Time: ' int2str(t) '. Convergence of the perfect foresight model solver!'])
|
2012-02-04 16:26:22 +01:00
|
|
|
end
|
2012-02-10 12:56:52 +01:00
|
|
|
else
|
|
|
|
if t<10
|
|
|
|
disp(['Time: ' int2str(t) '. No convergence of the perfect foresight model solver!'])
|
|
|
|
elseif t<100
|
|
|
|
disp(['Time: ' int2str(t) '. No convergence of the perfect foresight model solver!'])
|
|
|
|
elseif t<1000
|
|
|
|
disp(['Time: ' int2str(t) '. No convergence of the perfect foresight model solver!'])
|
|
|
|
else
|
|
|
|
disp(['Time: ' int2str(t) '. No convergence of the perfect foresight model solver!'])
|
2012-02-03 11:24:14 +01:00
|
|
|
end
|
2011-12-19 18:04:54 +01:00
|
|
|
end
|
2012-02-10 12:56:52 +01:00
|
|
|
end
|
|
|
|
if do_not_check_stability_flag
|
|
|
|
% Exit from the while loop.
|
|
|
|
endo_simul_1 = tmp;
|
|
|
|
break
|
|
|
|
else
|
|
|
|
% Test if periods is big enough.
|
|
|
|
% Increase the number of periods.
|
|
|
|
periods1 = periods1 + ep.step;
|
|
|
|
pfm1.periods = periods1;
|
|
|
|
pfm1.i_upd = pfm1.ny+(1:pfm1.periods*pfm1.ny);
|
|
|
|
% Increment the counter.
|
|
|
|
increase_periods = increase_periods + 1;
|
2011-12-12 12:40:46 +01:00
|
|
|
if verbosity
|
2012-02-10 12:56:52 +01:00
|
|
|
if t<10
|
|
|
|
disp(['Time: ' int2str(t) '. I increase the number of periods to ' int2str(periods1) '.'])
|
|
|
|
elseif t<100
|
|
|
|
disp(['Time: ' int2str(t) '. I increase the number of periods to ' int2str(periods1) '.'])
|
|
|
|
elseif t<1000
|
|
|
|
disp(['Time: ' int2str(t) '. I increase the number of periods to ' int2str(periods1) '.'])
|
2011-12-05 10:58:39 +01:00
|
|
|
else
|
2012-02-10 12:56:52 +01:00
|
|
|
disp(['Time: ' int2str(t) '. I increase the number of periods to ' int2str(periods1) '.'])
|
2011-12-05 10:58:39 +01:00
|
|
|
end
|
|
|
|
end
|
2012-02-10 12:56:52 +01:00
|
|
|
if info_convergence
|
|
|
|
% If the previous call to the perfect foresight model solver exited
|
|
|
|
% announcing that the routine converged, adapt the size of endo_simul_1
|
|
|
|
% and exo_simul_1.
|
2012-03-08 12:49:17 +01:00
|
|
|
endo_simul_1 = [ tmp , repmat(steady_state,1,ep.step) ];
|
|
|
|
exo_simul_1 = [ exo_simul_1 ; zeros(ep.step,exo_nbr)];
|
2012-02-10 12:56:52 +01:00
|
|
|
tmp_old = tmp;
|
2012-02-04 16:19:15 +01:00
|
|
|
else
|
2012-02-10 12:56:52 +01:00
|
|
|
% If the previous call to the perfect foresight model solver exited
|
|
|
|
% announcing that the routine did not converge, then tmp=1... Maybe
|
|
|
|
% should change that, because in some circonstances it may usefull
|
|
|
|
% to know where the routine did stop, even if convergence was not
|
|
|
|
% achieved.
|
|
|
|
endo_simul_1 = [ endo_simul_1 , repmat(steady_state,1,ep.step) ];
|
2012-03-08 12:49:17 +01:00
|
|
|
exo_simul_1 = [ exo_simul_1 ; zeros(ep.step,exo_nbr)];
|
2012-02-10 12:56:52 +01:00
|
|
|
end
|
|
|
|
% Solve the perfect foresight model with an increased number of periods.
|
2012-03-12 10:34:28 +01:00
|
|
|
if bytecode_flag && ~options_.ep.stochastic.order
|
|
|
|
[flag,tmp] = bytecode('dynamic',endo_simul_1,exo_simul_1);
|
2012-02-10 12:56:52 +01:00
|
|
|
else
|
|
|
|
flag = 1;
|
|
|
|
end
|
|
|
|
if flag
|
2012-03-08 12:49:17 +01:00
|
|
|
if options_.ep.stochastic.order == 0
|
|
|
|
[flag,tmp,err] = solve_perfect_foresight_model(endo_simul_1,exo_simul_1,pfm1);
|
2012-03-05 12:35:55 +01:00
|
|
|
else
|
2014-01-17 21:49:04 +01:00
|
|
|
switch(options_.ep.stochastic.algo)
|
|
|
|
case 0
|
|
|
|
[flag,tmp] = ...
|
|
|
|
solve_stochastic_perfect_foresight_model(endo_simul_1,exo_simul_1,pfm1,options_.ep.stochastic.quadrature.nodes,options_.ep.stochastic.order);
|
|
|
|
case 1
|
|
|
|
[flag,tmp] = ...
|
2014-05-12 09:40:12 +02:00
|
|
|
solve_stochastic_perfect_foresight_model_1(endo_simul_1,exo_simul_1,options_,pfm1,options_.ep.stochastic.order);
|
2014-01-17 21:49:04 +01:00
|
|
|
end
|
2012-03-05 12:35:55 +01:00
|
|
|
end
|
2012-02-10 12:56:52 +01:00
|
|
|
end
|
|
|
|
info_convergence = ~flag;
|
|
|
|
if info_convergence
|
|
|
|
% If the solver achieved convergence, check that simulated paths did not
|
|
|
|
% change during the first periods.
|
|
|
|
% Compute the maximum deviation between old path and new path over the
|
|
|
|
% first periods
|
2012-03-08 12:49:17 +01:00
|
|
|
delta = max(max(abs(tmp(:,2)-tmp_old(:,2))));
|
2012-02-10 12:56:52 +01:00
|
|
|
if delta < dynatol.x
|
|
|
|
% If the maximum deviation is close enough to zero, reset the number
|
|
|
|
% of periods to ep.periods
|
|
|
|
periods1 = ep.periods;
|
|
|
|
pfm1.periods = periods1;
|
|
|
|
pfm1.i_upd = pfm1.ny+(1:pfm1.periods*pfm1.ny);
|
|
|
|
% Cut exo_simul_1 and endo_simul_1 consistently with the resetted
|
|
|
|
% number of periods and exit from the while loop.
|
|
|
|
exo_simul_1 = exo_simul_1(1:(periods1+2),:);
|
|
|
|
endo_simul_1 = endo_simul_1(:,1:(periods1+2));
|
|
|
|
break
|
2012-02-03 11:24:14 +01:00
|
|
|
end
|
2012-02-10 12:56:52 +01:00
|
|
|
else
|
|
|
|
% The solver did not converge... Try to solve the model again with a bigger
|
|
|
|
% number of periods, except if the number of periods has been increased more
|
|
|
|
% than 10 times.
|
|
|
|
if increase_periods==10;
|
|
|
|
if verbosity
|
|
|
|
if t<10
|
|
|
|
disp(['Time: ' int2str(t) '. Even with ' int2str(periods1) ', I am not able to solve the perfect foresight model. Use homotopy instead...'])
|
|
|
|
elseif t<100
|
|
|
|
disp(['Time: ' int2str(t) '. Even with ' int2str(periods1) ', I am not able to solve the perfect foresight model. Use homotopy instead...'])
|
|
|
|
elseif t<1000
|
|
|
|
disp(['Time: ' int2str(t) '. Even with ' int2str(periods1) ', I am not able to solve the perfect foresight model. Use homotopy instead...'])
|
|
|
|
else
|
|
|
|
disp(['Time: ' int2str(t) '. Even with ' int2str(periods1) ', I am not able to solve the perfect foresight model. Use homotopy instead...'])
|
2011-12-05 10:58:39 +01:00
|
|
|
end
|
2012-01-27 18:24:24 +01:00
|
|
|
end
|
2012-02-10 12:56:52 +01:00
|
|
|
% Exit from the while loop.
|
|
|
|
break
|
2012-01-27 18:24:24 +01:00
|
|
|
end
|
2012-02-10 12:56:52 +01:00
|
|
|
end% if info_convergence
|
|
|
|
end
|
|
|
|
end% while
|
2014-06-16 14:25:23 +02:00
|
|
|
if ~info_convergence && ep.homotopic_steps % If exited from the while loop without achieving
|
|
|
|
% convergence use an homotopic approach
|
2012-03-02 11:33:29 +01:00
|
|
|
if ~do_not_check_stability_flag
|
|
|
|
periods1 = ep.periods;
|
|
|
|
pfm1.periods = periods1;
|
|
|
|
pfm1.i_upd = i_upd;
|
|
|
|
exo_simul_1 = exo_simul_1(1:(periods1+2),:);
|
|
|
|
endo_simul_1 = endo_simul_1(:,1:(periods1+2));
|
|
|
|
end
|
2012-03-08 12:49:17 +01:00
|
|
|
[INFO,tmp] = homotopic_steps(endo_simul,exo_simul_1,.5,.01,pfm1);
|
2012-03-02 11:33:29 +01:00
|
|
|
if isstruct(INFO)
|
|
|
|
info_convergence = INFO.convergence;
|
|
|
|
else
|
|
|
|
info_convergence = 0;
|
|
|
|
end
|
|
|
|
if ~info_convergence
|
2012-03-08 12:49:17 +01:00
|
|
|
[INFO,tmp] = homotopic_steps(endo_simul,exo_simul_1,0,.01,pfm1);
|
2012-03-02 11:33:29 +01:00
|
|
|
if isstruct(INFO)
|
|
|
|
info_convergence = INFO.convergence;
|
|
|
|
else
|
|
|
|
info_convergence = 0;
|
|
|
|
end
|
2012-02-10 12:56:52 +01:00
|
|
|
if ~info_convergence
|
|
|
|
disp('Homotopy:: No convergence of the perfect foresight model solver!')
|
|
|
|
error('I am not able to simulate this model!');
|
2011-12-19 18:04:54 +01:00
|
|
|
else
|
2012-02-07 16:30:36 +01:00
|
|
|
endo_simul_1 = tmp;
|
|
|
|
if verbosity && info_convergence
|
2011-12-19 18:04:54 +01:00
|
|
|
disp('Homotopy:: Convergence of the perfect foresight model solver!')
|
|
|
|
end
|
|
|
|
end
|
2012-02-10 12:56:52 +01:00
|
|
|
else
|
|
|
|
info_convergence = 1;
|
|
|
|
endo_simul_1 = tmp;
|
|
|
|
if verbosity && info_convergence
|
|
|
|
disp('Homotopy:: Convergence of the perfect foresight model solver!')
|
|
|
|
end
|
2011-12-05 10:58:39 +01:00
|
|
|
end
|
|
|
|
end
|
2014-06-16 14:25:23 +02:00
|
|
|
if info_convergence
|
|
|
|
% Save results of the perfect foresight model solver.
|
|
|
|
time_series(:,tsimul) = endo_simul_1(:,2);
|
|
|
|
endo_simul_1(:,1:end-1) = endo_simul_1(:,2:end);
|
|
|
|
endo_simul_1(:,1) = time_series(:,tsimul);
|
|
|
|
endo_simul_1(:,end) = oo_.steady_state;
|
|
|
|
tsimul = tsimul+1;
|
|
|
|
else
|
|
|
|
oo_.ep.failures.periods = [oo_.ep.failures.periods t];
|
|
|
|
oo_.ep.failures.previous_period = [oo_.ep.failures.previous_period endo_simul_1(:,1)];
|
|
|
|
oo_.ep.failures.shocks = [oo_.ep.failures.shocks shocks];
|
|
|
|
end
|
2012-02-03 11:24:14 +01:00
|
|
|
end% (while) loop over t
|
2011-12-12 14:01:17 +01:00
|
|
|
|
2012-01-11 16:00:42 +01:00
|
|
|
dyn_waitbar_close(hh);
|
2011-12-19 18:04:54 +01:00
|
|
|
|
2013-12-09 12:32:20 +01:00
|
|
|
if nargout
|
|
|
|
ts = dseries(transpose([initial_conditions, time_series]),options_.initial_period,cellstr(M_.endo_names));
|
|
|
|
else
|
2013-12-09 11:07:44 +01:00
|
|
|
oo_.endo_simul = [initial_conditions, time_series];
|
2013-12-09 12:32:20 +01:00
|
|
|
ts = dseries(transpose(oo_.endo_simul),options_.initial_period,cellstr(M_.endo_names));
|
2013-12-09 11:07:44 +01:00
|
|
|
dyn2vec;
|
2013-12-09 12:32:20 +01:00
|
|
|
end
|
|
|
|
|
|
|
|
assignin('base', 'Simulated_time_series', ts);
|