From 13ca15a2783af77fe8c1b75ff77adfa0e26bb44e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Hermes=29?=
<stephane.adjemian@univ-lemans.fr>
Date: Mon, 14 Mar 2016 20:11:33 +0100
Subject: [PATCH] Allow (S)EP with arbitrary sequence of innovations.
The third input argument of extended_path Matlab/Octave's routine is the
sequence of shocks (T*n array, where n is the number of exogenous
variables and T is the size of the sample). If the third argument is
empty, the (stochastic) extended path is run with gaussian
innovations (this corresponds to the previous behaviour).
TODO:
- Fix the compatibility with ep.replic_nbr
- Check the 'calibrated' mode.
---
matlab/ep/extended_path.m | 63 +++++++++++++++++++---------------
preprocessor/ComputingTasks.cc | 2 +-
tests/ep/ar.mod | 2 +-
tests/ep/burnside.mod | 4 +--
tests/ep/linearmodel.mod | 4 +--
tests/ep/rbc.mod | 10 +++---
tests/ep/rbcii.mod | 4 +--
7 files changed, 49 insertions(+), 40 deletions(-)
diff --git a/matlab/ep/extended_path.m b/matlab/ep/extended_path.m
index 1c73e474d2..70fc89de28 100644
--- a/matlab/ep/extended_path.m
+++ b/matlab/ep/extended_path.m
@@ -1,4 +1,4 @@
-function [ts,results] = extended_path(initial_conditions,sample_size, DynareOptions, DynareModel, DynareResults)
+function [ts,results] = extended_path(initial_conditions,sample_size, exogenousvariables, 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.
%
@@ -62,7 +62,6 @@ if isempty(initial_conditions)
end
end
-
% Set the number of periods for the perfect foresight model
periods = ep.periods;
pfm.periods = periods;
@@ -97,12 +96,14 @@ DynareOptions.minimal_solving_period = 100;%DynareOptions.ep.periods;
time_series = zeros(DynareModel.endo_nbr,sample_size);
% Set the covariance matrix of the structural innovations.
-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 = DynareModel.Sigma_e(positive_var_indx,positive_var_indx);
-covariance_matrix_upper_cholesky = chol(covariance_matrix);
+if isempty(exogenousvariables)
+ 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 = DynareModel.Sigma_e(positive_var_indx,positive_var_indx);
+ covariance_matrix_upper_cholesky = chol(covariance_matrix);
+end
% (re)Set exo_nbr
%exo_nbr = effective_number_of_shocks;
@@ -118,27 +119,35 @@ bytecode_flag = ep.use_bytecode;
replic_nbr = ep.replic_nbr;
% Simulate shocks.
-switch ep.innovation_distribution
- case 'gaussian'
- shocks = transpose(transpose(covariance_matrix_upper_cholesky)* ...
- randn(effective_number_of_shocks,sample_size* ...
- replic_nbr));
- shocks(:,positive_var_indx) = shocks;
- case 'calibrated'
- replic_nbr = 1;
- 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;
+if isempty(exogenousvariables)
+ switch ep.innovation_distribution
+ case 'gaussian'
+ shocks = transpose(transpose(covariance_matrix_upper_cholesky)* ...
+ randn(effective_number_of_shocks,sample_size* ...
+ replic_nbr));
+ shocks(:,positive_var_indx) = shocks;
+ case 'calibrated'
+ replic_nbr = 1;
+ 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;
+ end
end
+ otherwise
+ error(['extended_path:: ' ep.innovation_distribution ' distribution for the structural innovations is not (yet) implemented!'])
end
- otherwise
- error(['extended_path:: ' ep.innovation_distribution ' distribution for the structural innovations is not (yet) implemented!'])
+else
+ shocks = exogenousvariables;
+ testnonzero = abs(shocks)>0;
+ testnonzero = sum(testnonzero);
+ positive_var_indx = find(testnonzero);
+ effective_number_of_shocks = length(positive_var_indx);
end
% Set waitbar (graphic or text mode)
diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc
index 24c3f12d3f..6618dce0df 100644
--- a/preprocessor/ComputingTasks.cc
+++ b/preprocessor/ComputingTasks.cc
@@ -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
- << ", options_, M_, oo_);" << endl;
+ << ", [], options_, M_, oo_);" << endl;
}
ModelDiagnosticsStatement::ModelDiagnosticsStatement()
diff --git a/tests/ep/ar.mod b/tests/ep/ar.mod
index 4cc34cf4ad..9582f5aada 100644
--- a/tests/ep/ar.mod
+++ b/tests/ep/ar.mod
@@ -44,7 +44,7 @@ options_.ep.IntegrationAlgorithm='Tensor-Gaussian-Quadrature';
options_.ep.stochastic.nodes = 3;
options_.console_mode = 0;
-sts = extended_path([], 10, options_, M_, oo_);
+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)!!!')
diff --git a/tests/ep/burnside.mod b/tests/ep/burnside.mod
index cd863b8215..4d0bb6dc37 100644
--- a/tests/ep/burnside.mod
+++ b/tests/ep/burnside.mod
@@ -50,12 +50,12 @@ options_.ep.stochastic.nodes = 2;
options_.console_mode = 0;
set_dynare_seed('default');
-ts = extended_path([], 5000, options_, M_, oo_);
+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, options_, M_, oo_);
+ts1_4 = extended_path([], 5000, [], options_, M_, oo_);
set_dynare_seed('default');
ytrue=exact_solution(M_,oo_, 800);
diff --git a/tests/ep/linearmodel.mod b/tests/ep/linearmodel.mod
index a36902dd6c..3cadb061ae 100644
--- a/tests/ep/linearmodel.mod
+++ b/tests/ep/linearmodel.mod
@@ -33,13 +33,13 @@ options_.ep.order = 0;
options_.ep.nnodes = 0;
options_.console_mode = 0;
-ts = extended_path([], 10, options_, M_, oo_);
+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, options_, M_, oo_);
+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')
diff --git a/tests/ep/rbc.mod b/tests/ep/rbc.mod
index 4c2e3dc475..a03d0a5f49 100644
--- a/tests/ep/rbc.mod
+++ b/tests/ep/rbc.mod
@@ -75,19 +75,19 @@ steady(nocheck);
options_.ep.verbosity = 0;
options_.ep.stochastic.order = 0;
-ts0 = extended_path([], 10, options_, M_, oo_);
+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, options_, M_, oo_);
+ts1_3 = extended_path([], 10, [], options_, M_, oo_);
options_.ep.stochastic.nodes = 5;
-ts1_5 = extended_path([], 10, options_, M_, oo_);
+ts1_5 = extended_path([], 10, [], options_, M_, oo_);
options_.ep.stochastic.order = 2;
options_.ep.stochastic.nodes = 3;
-ts2_3 = extended_path([], 10, options_, M_, oo_);
+ts2_3 = extended_path([], 10, [], options_, M_, oo_);
options_.ep.stochastic.nodes = 5;
-ts2_5 = extended_path([], 10, options_, M_, oo_);
+ts2_5 = extended_path([], 10, [], options_, M_, oo_);
diff --git a/tests/ep/rbcii.mod b/tests/ep/rbcii.mod
index 2fa42bffa1..a51c01a8f4 100644
--- a/tests/ep/rbcii.mod
+++ b/tests/ep/rbcii.mod
@@ -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, options_, M_, oo_);
+ 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, options_, M_, oo_);
+ ts1_4 = extended_path([], 20, [], options_, M_, oo_);
// profile off
// profile viewer
@#else
--
GitLab