function time_series = extended_path(initial_conditions,sample_size)
% Stochastic simulation of a non linear DSGE model using the Extended Path method (Fair and Taylor 1983). A time
% series of size T is obtained by solving T perfect foresight models.
%
% 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.
%
% OUTPUTS
% o time_series [double] m*sample_size array, the simulations.
%
% ALGORITHM
%
% SPECIAL REQUIREMENTS
% Copyright (C) 2009, 2010, 2011 Dynare Team
%
% 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 .
global M_ options_ oo_
debug = 0;
options_.verbosity = options_.ep.verbosity;
verbosity = options_.ep.verbosity+debug;
% Test if bytecode and block options are used (these options are mandatory)
if ~( options_.bytecode && options_.block )
error('extended_path:: Options bytecode and block are mandatory!')
end
% Set default initial conditions.
if isempty(initial_conditions)
initial_conditions = oo_.steady_state;
end
% Set maximum number of iterations for the deterministic solver.
options_.maxit_ = options_.ep.maxit;
% Set the number of periods for the perfect foresight model
options_.periods = options_.ep.periods;
% Set the algorithm for the perfect foresight solver
options_.stack_solve_algo = options_.ep.stack_solve_algo;
% 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;
if options_.ep.init
lrf = load('linear_reduced_form','oo_');
oo_.dr = lrf.oo_.dr; clear('lrf');
if options_.ep.init==2
lambda = .8;
end
end
% Do not use a minimal number of perdiods for the perfect foresight solver (with bytecode and blocks)
options_.minimal_solving_period = options_.ep.periods;
% Get indices of variables with non zero steady state
idx = find(abs(oo_.steady_state)>0);
% Initialize the exogenous variables.
make_ex_;
% Initialize the endogenous variables.
make_y_;
% Initialize the output array.
time_series = zeros(M_.endo_nbr,sample_size);
% Set the covariance matrix of the structural innovations.
variances = diag(M_.Sigma_e);
positive_var_indx = find(variances>0);
effective_number_of_shocks = length(positive_var_indx);
stdd = sqrt(variances(positive_var_indx));
covariance_matrix = M_.Sigma_e(positive_var_indx,positive_var_indx);
covariance_matrix_upper_cholesky = chol(covariance_matrix);
% Set seed.
if options_.ep.set_dynare_seed_to_default
set_dynare_seed('default');
end
% Simulate shocks.
switch options_.ep.innovation_distribution
case 'gaussian'
oo_.ep.shocks = randn(sample_size,effective_number_of_shocks)*covariance_matrix_upper_cholesky;
otherwise
error(['extended_path:: ' options_.ep.innovation_distribution ' distribution for the structural innovations is not (yet) implemented!'])
end
% Set future shocks (Stochastic Extended Path approach)
if options_.ep.stochastic.status
switch options_.ep.stochastic.method
case 'tensor'
switch options_.ep.stochastic.ortpol
case 'hermite'
[r,w] = gauss_hermite_weights_and_nodes(options_.ep.stochastic.nodes);
otherwise
error('extended_path:: Unknown orthogonal polynomial option!')
end
if options_.ep.stochastic.order*M_.exo_nbr>1
for i=1:options_.ep.stochastic.order*M_.exo_nbr
rr(k) = {r};
ww(k) = {w};
end
rrr = cartesian_product_of_sets(rr{:});
www = cartesian_product_of_sets(ww{:});
else
rrr = r;
www = w;
end
www = prod(www,2);
number_of_nodes = length(www);
relative_weights = www/max(www);
switch options_.ep.stochastic.pruned.status
case 1
jdx = find(relative_weights>options_.ep.stochastic.pruned.relative);
www = www(jdx);
www = www/sum(www);
rrr = rrr(jdx,:);
case 2
jdx = find(weights>options_.ep.stochastic.pruned.level);
www = www(jdx);
www = www/sum(www);
rrr = rrr(jdx,:);
otherwise
% Nothing to be done!
end
nnn = length(www);
otherwise
error('extended_path:: Unknown stochastic_method option!')
end
else
rrr = zeros(1,effective_number_of_shocks);
www = 1;
nnn = 1;
end
% Initializes some variables.
t = 0;
% Set waitbar (graphic or text mode)
hh = dyn_waitbar(0,'Please wait. Extended Path simulations...');
set(hh,'Name','EP simulations.');
% Main loop.
while (t