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

@@ -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)

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()

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)!!!')

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);

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')

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_);

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