From 9c6e219990d1f3a80a0301fb3c920865801ffaa3 Mon Sep 17 00:00:00 2001 From: Michel Juillard Date: Fri, 6 Feb 2015 12:36:09 +0100 Subject: [PATCH] new implementation for extended path --- matlab/ep/extended_path.m | 243 +++++++++++++++++++++------------ matlab/global_initialization.m | 6 +- 2 files changed, 162 insertions(+), 87 deletions(-) diff --git a/matlab/ep/extended_path.m b/matlab/ep/extended_path.m index b35eb478f..f6c4177d5 100644 --- a/matlab/ep/extended_path.m +++ b/matlab/ep/extended_path.m @@ -1,4 +1,4 @@ -function ts = extended_path(initial_conditions,sample_size) +function [ts,results] = 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. % @@ -32,24 +32,28 @@ function ts = extended_path(initial_conditions,sample_size) % along with Dynare. If not, see . global M_ options_ oo_ -options_.verbosity = options_.ep.verbosity; -verbosity = options_.ep.verbosity+options_.ep.debug; +ep = options_.ep; +options_.verbosity = ep.verbosity; +verbosity = ep.verbosity+ep.debug; % Set maximum number of iterations for the deterministic solver. -options_.simul.maxit = options_.ep.maxit; +options_.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_); +endo_nbr = M_.endo_nbr; exo_nbr = M_.exo_nbr; -ep = options_.ep; +maximum_lag = M_.maximum_lag; +maximum_lead = M_.maximum_lead; +epreplic_nbr = ep.replic_nbr; steady_state = oo_.steady_state; dynatol = options_.dynatol; % Set default initial conditions. if isempty(initial_conditions) if isempty(M_.endo_histval) - initial_conditions = oo_.steady_state; + initial_conditions = steady_state; else initial_conditions = M_.endo_histval; end @@ -57,18 +61,19 @@ end % Set the number of periods for the perfect foresight model -periods = options_.ep.periods; +periods = ep.periods; pfm.periods = periods; pfm.i_upd = pfm.ny+(1:pfm.periods*pfm.ny); +pfm.block = options_.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 = options_.ep.stack_solve_algo; +options_.stack_solve_algo = ep.stack_solve_algo; % Set check_stability flag -do_not_check_stability_flag = ~options_.ep.check_stability; +do_not_check_stability_flag = ~ep.check_stability; % Compute the first order reduced form if needed. % @@ -76,7 +81,7 @@ do_not_check_stability_flag = ~options_.ep.check_stability; % all the globals in a mat file called linear_reduced_form.mat; dr = struct(); -if options_.ep.init +if ep.init options_.order = 1; [dr,Info,M_,options_,oo_] = resol(1,M_,options_,oo_); end @@ -85,12 +90,10 @@ end options_.minimal_solving_period = 100;%options_.ep.periods; % Initialize the exogenous variables. -% !!!!!!!! Needs to fixed -options_.periods = periods; make_ex_; % Initialize the endogenous variables. -make_y_; +%make_y_; % Initialize the output array. time_series = zeros(M_.endo_nbr,sample_size); @@ -107,19 +110,36 @@ covariance_matrix_upper_cholesky = chol(covariance_matrix); %exo_nbr = effective_number_of_shocks; % Set seed. -if options_.ep.set_dynare_seed_to_default +if ep.set_dynare_seed_to_default set_dynare_seed('default'); end % Set bytecode flag -bytecode_flag = options_.ep.use_bytecode; +bytecode_flag = ep.use_bytecode; +% Set number of replications +replic_nbr = ep.replic_nbr; % Simulate shocks. -switch options_.ep.innovation_distribution +switch ep.innovation_distribution case 'gaussian' - oo_.ep.shocks = transpose(transpose(covariance_matrix_upper_cholesky)*randn(effective_number_of_shocks,sample_size)); + shocks = transpose(transpose(covariance_matrix_upper_cholesky)* ... + randn(effective_number_of_shocks,sample_size*replic_nbr)); + case 'calibrated' + replic_nbr = 1; + shocks = zeros(sample_size,effective_number_of_shocks); + 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(k,ivar) = v; + else + socks(k,ivar) = shocks(k,ivar) * v; + end + end + shocks = shocks(:,positive_var_indx); otherwise - error(['extended_path:: ' options_.ep.innovation_distribution ' distribution for the structural innovations is not (yet) implemented!']) + error(['extended_path:: ' ep.innovation_distribution ' distribution for the structural innovations is not (yet) implemented!']) end @@ -128,7 +148,7 @@ hh = dyn_waitbar(0,'Please wait. Extended Path simulations...'); set(hh,'Name','EP simulations.'); % hybrid correction -pfm.hybrid_order = options_.ep.stochastic.hybrid_order; +pfm.hybrid_order = ep.stochastic.hybrid_order; if pfm.hybrid_order oo_.dr = set_state_space(oo_.dr,M_,options_); options = options_; @@ -142,16 +162,16 @@ end pfm.nnzA = M_.NNZDerivatives(1); % setting up integration nodes if order > 0 -if options_.ep.stochastic.order > 0 +if ep.stochastic.order > 0 [nodes,weights,nnodes] = setup_integration_nodes(options_.ep,pfm); pfm.nodes = nodes; pfm.weights = weights; pfm.nnodes = nnodes; % compute number of blocks - [block_nbr,pfm.world_nbr] = get_block_world_nbr(options_.ep.stochastic.algo,nnodes,options_.ep.stochastic.order,options_.ep.periods); + [block_nbr,pfm.world_nbr] = get_block_world_nbr(ep.stochastic.algo,nnodes,ep.stochastic.order,ep.periods); else - block_nbr = options_.ep.periods + block_nbr = ep.periods; end @@ -167,74 +187,54 @@ oo_.ep.failures.previous_period = cell(0); oo_.ep.failures.shocks = cell(0); % Initializes some variables. -t = 0; +t = 1; tsimul = 1; +nx = length(oo_.exo_simul); +for k = 1:replic_nbr + results{k} = zeros(endo_nbr,sample_size+1); + results{k}(:,1) = initial_conditions; +end % Main loop. -while (t 1 && ep.parallel_1 + parfor k = 1:replic_nbr + exo_simul = repmat(oo_.exo_steady_state',periods+2,1); + exo_simul(1:sample_size+3-t,:) = oo_.exo_simul(t:end,:); + exo_simul(2,positive_var_indx) = exo_simul(2+1,positive_var_indx) + ... + shocks((t-2)*replic_nbr+k,:); + initial_conditions = results{k}(:,t-1); + results{k}(:,t) = extended_path_core(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.stock_solve_algo,... + options_.lmmcp); + end else - if t==1 - endo_simul_1 = repmat(steady_state,1,periods+2); + for k = 1:replic_nbr + exo_simul = repmat(oo_.exo_steady_state',periods+maximum_lag+ ... + maximum_lead,1); + exo_simul(1:sample_size+maximum_lag+maximum_lead-t+1,:) = ... + oo_.exo_simul(t:end,:); + exo_simul(maximum_lag+1,positive_var_indx) = ... + exo_simul(maximum_lag+1,positive_var_indx) + shocks((t-2)*replic_nbr+k,:); + initial_conditions = results{k}(:,t-1); + results{k}(:,t) = extended_path_core(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_,pfm,ep.stochastic.algo,ep.stack_solve_algo,... + options_.lmmcp); end end - % Solve a perfect foresight model. - % Keep a copy of endo_simul_1 - endo_simul = endo_simul_1; - if verbosity - save ep_test_1 endo_simul_1 exo_simul_1 - end - if bytecode_flag && ~options_.ep.stochastic.order - [flag,tmp] = bytecode('dynamic',endo_simul_1,exo_simul_1, M_.params, endo_simul_1, options_.ep.periods); - else - flag = 1; - end - if flag - if options_.ep.stochastic.order == 0 - [flag,tmp,err] = solve_perfect_foresight_model(endo_simul_1,exo_simul_1,pfm); - else - switch(options_.ep.stochastic.algo) - case 0 - [flag,tmp] = ... - solve_stochastic_perfect_foresight_model(endo_simul_1,exo_simul_1,pfm,options_.ep.stochastic.quadrature.nodes,options_.ep.stochastic.order); - case 1 - [flag,tmp] = ... - solve_stochastic_perfect_foresight_model_1(endo_simul_1,exo_simul_1,options_,pfm,options_.ep.stochastic.order); - end - end - end - info_convergence = ~flag; - if verbosity - if info_convergence - disp(['Time: ' int2str(t) '. Convergence of the perfect foresight model solver!']) - else - disp(['Time: ' int2str(t) '. No convergence of the perfect foresight model solver!']) - end - end - endo_simul_1 = tmp; - 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]; - endo_simul_1 = repmat(steady_state,1,periods+2); - endo_simul_1(:,1) = time_series(:,tsimul-1); - end + + end% (while) loop over t dyn_waitbar_close(hh); @@ -242,14 +242,85 @@ dyn_waitbar_close(hh); if isnan(options_.initial_period) initial_period = dates(1,1); else - initial_period = optins_.initial_period; + initial_period = options_.initial_period; end if nargout - ts = dseries(transpose([initial_conditions, time_series]),initial_period,cellstr(M_.endo_names)); + if ~isnan(results{1}) + ts = dseries(transpose([results{1}]), ... + initial_period,cellstr(M_.endo_names)); + else + ts = NaN; + end else - oo_.endo_simul = [initial_conditions, time_series]; - ts = dseries(transpose(oo_.endo_simul),initial_period,cellstr(M_.endo_names)); - dyn2vec; + if ~isnan(results{1}) + oo_.endo_simul = results{1}; + ts = dseries(transpose(results{1}),initial_period, ... + cellstr(M_.endo_names)); + else + oo_.endo_simul = NaN; + ts = NaN; + end end - assignin('base', 'Simulated_time_series', ts); \ No newline at end of file + assignin('base', 'Simulated_time_series', ts); + + +function y = extended_path_core(periods,endo_nbr,exo_nbr,positive_var_indx, ... + exo_simul,init,initial_conditions,... + maximum_lag,maximum_lead,steady_state, ... + verbosity,bytecode_flag,order,M,pfm,algo,stack_solve_algo,... + olmmcp) + +if init% Compute first order solution (Perturbation)... + endo_simul = simult_(initial_conditions,dr,exo_simul(2:end,:),1); +else + endo_simul = [initial_conditions repmat(steady_state,1,periods+1)]; +end +oo.endo_simul = endo_simul; +oo_.endo_simul = endo_simul; +% Solve a perfect foresight model. +% Keep a copy of endo_simul_1 +if verbosity + save ep_test_1 endo_simul exo_simul +end +if bytecode_flag && ~ep.stochastic.order + [flag,tmp] = bytecode('dynamic',endo_simul,exo_simul, M_.params, endo_simul, ep.periods); +else + flag = 1; +end +if flag + if order == 0 + options.periods = periods; + options.block = pfm.block; + oo.endo_simul = endo_simul; + oo.exo_simul = exo_simul; + oo.steady_state = steady_state; + options.bytecode = bytecode_flag; + options.lmmcp = olmmcp; + options.stack_solve_algo = stack_solve_algo; + [tmp,flag] = perfect_foresight_solver_core(M,oo,options); + else + switch(algo) + case 0 + [flag,tmp] = ... + solve_stochastic_perfect_foresight_model(endo_simul,exo_simul,pfm,ep.stochastic.quadrature.nodes,ep.stochastic.order); + case 1 + [flag,tmp] = ... + solve_stochastic_perfect_foresight_model_1(endo_simul,exo_simul,options_,pfm,ep.stochastic.order); + end + end +end +info_convergence = ~flag; +if verbosity + if info_convergence + disp(['Time: ' int2str(t) '. Convergence of the perfect foresight model solver!']) + else + disp(['Time: ' int2str(t) '. No convergence of the perfect foresight model solver!']) + end +end +endo_simul = tmp; +if info_convergence + y = endo_simul(:,2); +else + y = NaN(size(endo_nbr,1)); +end \ No newline at end of file diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m index ebe2a04d2..bd9e9a5ae 100644 --- a/matlab/global_initialization.m +++ b/matlab/global_initialization.m @@ -187,7 +187,11 @@ ep.innovation_distribution = 'gaussian'; % Set flag for the seed ep.set_dynare_seed_to_default = 1; % Set algorithm for the perfect foresight solver -ep.stack_solve_algo = 4; +ep.stack_solve_algo = 7; +% Number of replications +ep.replic_nbr = 1; +% Parallel execution of replications +ep.parallel_1 = false; % Stochastic extended path related options. ep.stochastic.method = ''; ep.stochastic.algo = 0;