stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Removed globals from extended_path routine.

time-shift
Stéphane Adjemian (Hermes) 2016-03-11 13:54:15 +01:00 committed by Stéphane Adjemian (Charybdis)
parent 635d5b704b
commit b60bd7b36b
7 changed files with 70 additions and 71 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

111
matlab/ep/extended_path.m
View File

@ -1,4 +1,4 @@
function [ts,results] = extended_path(initial_conditions,sample_size)
function [ts,results] = extended_path(initial_conditions,sample_size, DynareOptions, DynareModel, DynareResults)
% 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.
%
@ -30,36 +30,35 @@ function [ts,results] = extended_path(initial_conditions,sample_size)
%
% 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_
ep = options_.ep;
options_.verbosity = ep.verbosity;
ep = DynareOptions.ep;
DynareOptions.verbosity = ep.verbosity;
verbosity = ep.verbosity+ep.debug;
% Set maximum number of iterations for the deterministic solver.
options_.simul.maxit = ep.maxit;
DynareOptions.simul.maxit = ep.maxit;
% Prepare a structure needed by the matlab implementation of the perfect foresight model solver
pfm = setup_stochastic_perfect_foresight_model_solver(M_,options_,oo_);
pfm = setup_stochastic_perfect_foresight_model_solver(DynareModel,DynareOptions,DynareResults);
if M_.exo_det_nbr~=0
if DynareModel.exo_det_nbr~=0
error('ep: Extended path does not support varexo_det.')
end
endo_nbr = M_.endo_nbr;
exo_nbr = M_.exo_nbr;
maximum_lag = M_.maximum_lag;
maximum_lead = M_.maximum_lead;
endo_nbr = DynareModel.endo_nbr;
exo_nbr = DynareModel.exo_nbr;
maximum_lag = DynareModel.maximum_lag;
maximum_lead = DynareModel.maximum_lead;
epreplic_nbr = ep.replic_nbr;
steady_state = oo_.steady_state;
dynatol = options_.dynatol;
steady_state = DynareResults.steady_state;
dynatol = DynareOptions.dynatol;
% Set default initial conditions.
if isempty(initial_conditions)
if isempty(M_.endo_histval)
if isempty(DynareModel.endo_histval)
initial_conditions = steady_state;
else
initial_conditions = M_.endo_histval;
initial_conditions = DynareModel.endo_histval;
end
end
@ -68,13 +67,13 @@ end
periods = ep.periods;
pfm.periods = periods;
pfm.i_upd = pfm.ny+(1:pfm.periods*pfm.ny);
pfm.block = options_.block;
pfm.block = DynareOptions.block;
% keep a copy of pfm.i_upd
i_upd = pfm.i_upd;
% Set the algorithm for the perfect foresight solver
options_.stack_solve_algo = ep.stack_solve_algo;
DynareOptions.stack_solve_algo = ep.stack_solve_algo;
% Set check_stability flag
do_not_check_stability_flag = ~ep.check_stability;
@ -86,23 +85,23 @@ do_not_check_stability_flag = ~ep.check_stability;
dr = struct();
if ep.init
options_.order = 1;
oo_.dr=set_state_space(dr,M_,options_);
[dr,Info,M_,options_,oo_] = resol(0,M_,options_,oo_);
DynareOptions.order = 1;
DynareResults.dr=set_state_space(dr,DynareModel,DynareOptions);
[dr,Info,DynareModel,DynareOptions,DynareResults] = resol(0,DynareModel,DynareOptions,DynareResults);
end
% Do not use a minimal number of perdiods for the perfect foresight solver (with bytecode and blocks)
options_.minimal_solving_period = 100;%options_.ep.periods;
DynareOptions.minimal_solving_period = 100;%DynareOptions.ep.periods;
% Initialize the output array.
time_series = zeros(M_.endo_nbr,sample_size);
time_series = zeros(DynareModel.endo_nbr,sample_size);
% Set the covariance matrix of the structural innovations.
variances = diag(M_.Sigma_e);
variances = diag(DynareModel.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 = DynareModel.Sigma_e(positive_var_indx,positive_var_indx);
covariance_matrix_upper_cholesky = chol(covariance_matrix);
% (re)Set exo_nbr
@ -127,12 +126,12 @@ switch ep.innovation_distribution
shocks(:,positive_var_indx) = shocks;
case 'calibrated'
replic_nbr = 1;
shocks = zeros(sample_size,M_.exo_nbr);
for i = 1:length(M_.unanticipated_det_shocks)
k = M_.unanticipated_det_shocks(i).periods;
ivar = M_.unanticipated_det_shocks(i).exo_id;
v = M_.unanticipated_det_shocks(i).value;
if ~M_.unanticipated_det_shocks(i).multiplicative
shocks = zeros(sample_size,DynareModel.exo_nbr);
for i = 1:length(DynareModel.unanticipated_det_shocks)
k = DynareModel.unanticipated_det_shocks(i).periods;
ivar = DynareModel.unanticipated_det_shocks(i).exo_id;
v = DynareModel.unanticipated_det_shocks(i).value;
if ~DynareModel.unanticipated_det_shocks(i).multiplicative
shocks(k,ivar) = v;
else
socks(k,ivar) = shocks(k,ivar) * v;
@ -149,20 +148,20 @@ set(hh,'Name','EP simulations.');
% hybrid correction
pfm.hybrid_order = ep.stochastic.hybrid_order;
if pfm.hybrid_order
oo_.dr = set_state_space(oo_.dr,M_,options_);
options = options_;
DynareResults.dr = set_state_space(DynareResults.dr,DynareModel,DynareOptions);
options = DynareOptions;
options.order = pfm.hybrid_order;
pfm.dr = resol(0,M_,options,oo_);
pfm.dr = resol(0,DynareModel,options,DynareResults);
else
pfm.dr = [];
end
% number of nonzero derivatives
pfm.nnzA = M_.NNZDerivatives(1);
pfm.nnzA = DynareModel.NNZDerivatives(1);
% setting up integration nodes if order > 0
if ep.stochastic.order > 0
[nodes,weights,nnodes] = setup_integration_nodes(options_.ep,pfm);
[nodes,weights,nnodes] = setup_integration_nodes(DynareOptions.ep,pfm);
pfm.nodes = nodes;
pfm.weights = weights;
pfm.nnodes = nnodes;
@ -175,17 +174,17 @@ end
% set boundaries if mcp
[lb,ub,pfm.eq_index] = get_complementarity_conditions(M_, options_.ramsey_policy);
options_.lmmcp.lb = repmat(lb,block_nbr,1);
options_.lmmcp.ub = repmat(ub,block_nbr,1);
[lb,ub,pfm.eq_index] = get_complementarity_conditions(DynareModel, DynareOptions.ramsey_policy);
DynareOptions.lmmcp.lb = repmat(lb,block_nbr,1);
DynareOptions.lmmcp.ub = repmat(ub,block_nbr,1);
pfm.block_nbr = block_nbr;
% storage for failed draws
oo_.ep.failures.periods = [];
oo_.ep.failures.previous_period = cell(0);
oo_.ep.failures.shocks = cell(0);
DynareResults.ep.failures.periods = [];
DynareResults.ep.failures.previous_period = cell(0);
DynareResults.ep.failures.shocks = cell(0);
oo_.exo_simul = shocks;
DynareResults.exo_simul = shocks;
% Initializes some variables.
t = 1;
@ -196,7 +195,7 @@ for k = 1:replic_nbr
end
%make_ex_;
exo_simul_ = zeros(maximum_lag+sample_size+maximum_lead,exo_nbr);
exo_simul_(1:size(oo_.exo_simul,1),1:size(oo_.exo_simul,2)) = oo_.exo_simul;
exo_simul_(1:size(DynareResults.exo_simul,1),1:size(DynareResults.exo_simul,2)) = DynareResults.exo_simul;
% Main loop.
while (t <= sample_size)
if ~mod(t,10)
@ -207,21 +206,21 @@ while (t <= sample_size)
if replic_nbr > 1 && ep.parallel_1
parfor k = 1:replic_nbr
exo_simul = repmat(oo_.exo_steady_state',periods+2,1);
exo_simul = repmat(DynareResults.exo_steady_state',periods+2,1);
% exo_simul(1:sample_size+3-t,:) = exo_simul_(t:end,:);
exo_simul(2,:) = exo_simul_(M_.maximum_lag+t,:) + ...
exo_simul(2,:) = exo_simul_(DynareModel.maximum_lag+t,:) + ...
shocks((t-2)*replic_nbr+k,:);
initial_conditions = results{k}(:,t-1);
[results{k}(:,t), info_convergence] = extended_path_core(ep.periods,endo_nbr,exo_nbr,positive_var_indx, ...
exo_simul,ep.init,initial_conditions,...
maximum_lag,maximum_lead,steady_state, ...
ep.verbosity,bytecode_flag,ep.stochastic.order,...
M_.params,pfm,ep.stochastic.algo,ep.solve_algo,ep.stack_solve_algo,...
options_.lmmcp,options_,oo_);
DynareModel.params,pfm,ep.stochastic.algo,ep.solve_algo,ep.stack_solve_algo,...
DynareOptions.lmmcp,DynareOptions,DynareResults);
end
else
for k = 1:replic_nbr
exo_simul = repmat(oo_.exo_steady_state',periods+maximum_lag+ ...
exo_simul = repmat(DynareResults.exo_steady_state',periods+maximum_lag+ ...
maximum_lead,1);
% exo_simul(1:sample_size+maximum_lag+maximum_lead-t+1,:) = ...
% exo_simul_(t:end,:);
@ -232,8 +231,8 @@ while (t <= sample_size)
exo_simul,ep.init,initial_conditions,...
maximum_lag,maximum_lead,steady_state, ...
ep.verbosity,bytecode_flag,ep.stochastic.order,...
M_,pfm,ep.stochastic.algo,ep.solve_algo,ep.stack_solve_algo,...
options_.lmmcp,options_,oo_);
DynareModel,pfm,ep.stochastic.algo,ep.solve_algo,ep.stack_solve_algo,...
DynareOptions.lmmcp,DynareOptions,DynareResults);
end
end
if verbosity
@ -247,25 +246,25 @@ end% (while) loop over t
dyn_waitbar_close(hh);
if isnan(options_.initial_period)
if isnan(DynareOptions.initial_period)
initial_period = dates(1,1);
else
initial_period = options_.initial_period;
initial_period = DynareOptions.initial_period;
end
if nargout
if ~isnan(results{1})
ts = dseries(transpose([results{1}]), ...
initial_period,cellstr(M_.endo_names));
initial_period,cellstr(DynareModel.endo_names));
else
ts = NaN;
end
else
if ~isnan(results{1})
oo_.endo_simul = results{1};
DynareResults.endo_simul = results{1};
ts = dseries(transpose(results{1}),initial_period, ...
cellstr(M_.endo_names));
cellstr(DynareModel.endo_names));
else
oo_.endo_simul = NaN;
DynareResults.endo_simul = NaN;
ts = NaN;
end
end

2
preprocessor/ComputingTasks.cc
View File

@ -3200,7 +3200,7 @@ ExtendedPathStatement::writeOutput(ostream &output, const string &basename, bool
output << "options_." << it->first << " = " << it->second << ";" << endl;
output << "extended_path([], " << options_list.num_options.find("periods")->second
<< ");" << endl;
<< ", options_, M_, oo_);" << endl;
}
ModelDiagnosticsStatement::ModelDiagnosticsStatement()

4
tests/ep/ar.mod
View File

@ -36,7 +36,7 @@ options_.ep.stochastic.order = 0;
options_.ep.stochastic.nodes = 0;
options_.console_mode = 0;
ts = extended_path([],10);
ts = extended_path([], 10, options_, M_, oo_);
options_.ep.verbosity = 0;
options_.ep.stochastic.order = 1;
@ -44,7 +44,7 @@ options_.ep.IntegrationAlgorithm='Tensor-Gaussian-Quadrature';
options_.ep.stochastic.nodes = 3;
options_.console_mode = 0;
sts = extended_path([],10);
sts = extended_path([], 10, options_, M_, oo_);
if max(max(abs(ts-sts)))>pi*options_.dynatol.x
disp('Stochastic Extended Path:: Something is wrong here (potential bug in extended_path.m)!!!')

6
tests/ep/burnside.mod
View File

@ -50,15 +50,15 @@ options_.ep.stochastic.nodes = 2;
options_.console_mode = 0;
set_dynare_seed('default');
ts = extended_path([],5000);
ts = extended_path([], 5000, options_, M_, oo_);
options_.ep.stochastic.order = 2;
options_.ep.IntegrationAlgorithm='Tensor-Gaussian-Quadrature';
set_dynare_seed('default');
ts1_4 = extended_path([],5000);
ts1_4 = extended_path([], 5000, options_, M_, oo_);
set_dynare_seed('default');
ytrue=exact_solution(M_,oo_,800);
ytrue=exact_solution(M_,oo_, 800);
disp('True mean and standard deviation')
disp(mean(ytrue(101:end)))

4
tests/ep/linearmodel.mod
View File

@ -33,13 +33,13 @@ options_.ep.order = 0;
options_.ep.nnodes = 0;
options_.console_mode = 0;
ts = extended_path([],10);
ts = extended_path([], 10, options_, M_, oo_);
options_.ep.stochastic.status = 1;
options_.ep.IntegrationAlgorithm='Tensor-Gaussian-Quadrature';
options_.ep.order = 1;
options_.ep.nnodes = 3;
sts = extended_path([],10);
sts = extended_path([], 10, options_, M_, oo_);
if max(max(abs(ts.data-sts.data))) > 1e-12
error('extended path algorithm fails in ./tests/ep/linearmodel.mod')

10
tests/ep/rbc.mod
View File

@ -75,19 +75,19 @@ steady(nocheck);
options_.ep.verbosity = 0;
options_.ep.stochastic.order = 0;
ts0 = extended_path([],10);
ts0 = extended_path([], 10, options_, M_, oo_);
options_.ep.stochastic.order = 1;
options_.ep.stochastic.nodes = 3;
options_.ep.IntegrationAlgorithm='Tensor-Gaussian-Quadrature';
ts1_3 = extended_path([],10);
ts1_3 = extended_path([], 10, options_, M_, oo_);
options_.ep.stochastic.nodes = 5;
ts1_5 = extended_path([],10);
ts1_5 = extended_path([], 10, options_, M_, oo_);
options_.ep.stochastic.order = 2;
options_.ep.stochastic.nodes = 3;
ts2_3 = extended_path([],10);
ts2_3 = extended_path([], 10, options_, M_, oo_);
options_.ep.stochastic.nodes = 5;
ts2_5 = extended_path([],10);
ts2_5 = extended_path([], 10, options_, M_, oo_);

4
tests/ep/rbcii.mod
View File

@ -72,12 +72,12 @@ copyfile('rbcii_steady_state.m','rbcii_steadystate2.m');
options_.ep.stochastic.nodes = 2;
options_.console_mode = 0;
ts = extended_path([],20);
ts = extended_path([], 20, options_, M_, oo_);
options_.ep.stochastic.order = 1;
options_.ep.IntegrationAlgorithm='Tensor-Gaussian-Quadrature';
// profile on
ts1_4 = extended_path([],20);
ts1_4 = extended_path([], 20, options_, M_, oo_);
// profile off
// profile viewer
@#else