Allow lags on exogenous variables for backward models forecasts and IRFs.

TODO: Check if this is working with linear models.
(cherry picked from commit 44b9531bb7)

matlab/backward/backward_model_forecast.m

@@ -79,11 +79,27 @@ if withuncertainty
     sigma = transpose(chol(Sigma));
 end
 
+% Set initial condition.
+if isdates(initialcondition)
+    if isempty(M_.endo_histval)
+        error('backward_model_irf: histval block for setting initial condition is missing!')
+    end
+    initialcondition = dseries(transpose(M_.endo_histval), initialcondition, endo_names, cellstr(M_.endo_names_tex));
+end
+
 % Put initial conditions in a vector of doubles
 initialconditions = transpose(initialcondition{endo_names{:}}.data);
 
 % Compute forecast without shock 
-innovations = zeros(periods+1, M_.exo_nbr);
+innovations = zeros(periods+max(M_.maximum_exo_lag, 1), M_.exo_nbr);
+if M_.maximum_exo_lag
+    if isempty(M_.exo_histval)
+        error('You need to set the past values for the exogenous variables!')
+    else
+        innovations(1:M_.maximum_exo_lag, :) = M_.exo_histval;
+    end
+end
+
 oo__0 = simul_backward_model(initialconditions, periods, options_, M_, oo_, innovations);
 forecasts.pointforecast = dseries(transpose(oo__0.endo_simul(idy,:)), initialcondition.init, listofvariables);
 
@@ -91,7 +107,7 @@ if withuncertainty
     % Preallocate an array gathering the simulations.
     ArrayOfForectasts = zeros(n, periods+1, B);
     for i=1:B
-        innovations(2:end,:) = transpose(sigma*randn(M_.exo_nbr, periods));
+        innovations(max(M_.maximum_exo_lag, 1)+1:end,:) = transpose(sigma*randn(M_.exo_nbr, periods));
         oo__ = simul_backward_model(initialconditions, periods, options_, M_, oo_, innovations);     
         ArrayOfForecasts(:,:,i) = oo__.endo_simul(idy,:); 
     end

matlab/backward/backward_model_irf.m

@@ -1,9 +1,9 @@
-function irfs = backward_model_irf(initialcondition, listofshocks, listofvariables, periods, transform)
+function irfs = backward_model_irf(initialcondition, listofshocks, listofvariables, varargin)
 
 % Returns impulse response functions.
 %
 % INPUTS 
-% - initialcondition    [dseries]             Initial conditions for the endogenous variables.
+% - initialcondition    [dseries,dates]       Initial conditions for the endogenous variables, or period 0.
 % - listofshocks        [cell of strings]     The innovations for which the IRFs need to be computed.  
 % - listofvariables     [cell of strings]     The endogenous variables which will be returned.  
 % - periods             [integer]             scalar, the number of periods.
@@ -42,19 +42,12 @@ if M_.maximum_lead
     error(['simul_model_irf:: The specified model is not backward looking!'])
 end
 
-% Set the number of innovations and variables for which we want to compute the IRFs.
-nx = length(listofshocks);
-ny = length(listofvariables);
-
-% Initialization of the returned argument. Each will be a dseries
-% object containing the IRFS for the endogenous variables listed in the
-% third input argument.
-irfs = struct();
-
 % Set default value for the fourth input argument.
 if nargin<4
     periods = 40;
     notransform = true;
+else
+    periods = varargin{1};
 end
 
 % Set default value for the last input argument (no transformation).
@@ -62,6 +55,7 @@ if nargin<5
     notransform = true;
 else
     notransform = false;
+    transform = varargin{2};
 end
 
 % Get list of all exogenous variables in a cell of strings
@@ -85,20 +79,33 @@ sigma = transpose(chol(Sigma));
 % Put initial conditions in a vector of doubles
 initialconditions = transpose(initialcondition{endo_names{:}}.data);
 
+% Initialization of the returned argument. Each will be a dseries object containing the IRFS for the endogenous variables listed in the third input argument.
+irfs = struct();
+
+% Get the covariance matrix of the shocks.
+Sigma = M_.Sigma_e + 1e-14*eye(M_.exo_nbr);
+sigma = transpose(chol(Sigma));
+
+% Put initial conditions in a vector of doubles
+initialconditions = transpose(initialcondition{endo_names{:}}.data);
+
 % Compute the IRFs (loop over innovations).
 for i=1:length(listofshocks)
     % Get transition paths induced by the initial condition.
-    innovations = zeros(periods+1, M_.exo_nbr);
+    innovations = zeros(periods+M_.maximum_exo_lag, M_.exo_nbr);
+    if ~isempty(M_.exo_histval)
+        innovations(1:M_.maximum_exo_lag,:) = M_.exo_histval;
+    end
     oo__0 = simul_backward_model(initialconditions, periods, options_, M_, oo_, innovations);
     % Add the shock.
     j = strmatch(listofshocks{i}, exo_names);
     if isempty(j)
         error('backward_model_irf: Exogenous variable %s is unknown!', listofshocks{i})
     end
-    innovations(2,:) = transpose(sigma(:,j));
+    innovations(1+M_.maximum_exo_lag,:) = transpose(sigma(:,j));
     oo__1 = simul_backward_model(initialconditions, periods, options_, M_, oo_, innovations);
     % Transform the endogenous variables
-        if notransform
+    if notransform
         endo_simul__0 = oo__0.endo_simul;
         endo_simul__1 = oo__1.endo_simul;
     else

matlab/backward/simul_backward_model.m

@@ -76,10 +76,14 @@ if nargin<6
       otherwise
         error(['simul_backward_nonlinear_model:: ' DynareOption.bnlms.innovation_distribution ' distribution for the structural innovations is not (yet) implemented!'])
     end
-
     % Put the simulated innovations in DynareOutput.exo_simul.
-    DynareOutput.exo_simul = zeros(sample_size+1,number_of_shocks);
+    DynareOutput.exo_simul = zeros(sample_size, number_of_shocks);
     DynareOutput.exo_simul(2:end,positive_var_indx) = DynareOutput.bnlms.shocks;
+    if isfield(DynareModel,'exo_histval') && ~isempty(DynareModel.exo_histval)
+        DynareOutput.exo_simul = [M_.exo_histval; DynareOutput.exo_simul];
+    else
+        DynareOutput.exo_simul = [zeros(1,number_of_shocks); DynareOutput.exo_simul];
+    end
 else
     number_of_shocks = size(innovations,2);
     DynareOutput.exo_simul = innovations;

matlab/backward/simul_backward_nonlinear_model.m

@@ -64,7 +64,11 @@ if nargin<6
     % Put the simulated innovations in DynareOutput.exo_simul.
     DynareOutput.exo_simul = zeros(sample_size,number_of_shocks);
     DynareOutput.exo_simul(:,positive_var_indx) = DynareOutput.bnlms.shocks;
-    DynareOutput.exo_simul = [zeros(1,number_of_shocks); DynareOutput.exo_simul];
+    if isfield(DynareModel,'exo_histval') && ~ isempty(DynareModel.exo_histval)
+        DynareOutput.exo_simul = [M_.exo_histval; DynareOutput.exo_simul];
+    else
+        DynareOutput.exo_simul = [zeros(1,number_of_shocks); DynareOutput.exo_simul];
+    end
 else
     DynareOutput.exo_simul = innovations;
 end
@@ -86,7 +90,7 @@ y = NaN(length(idx)+ny1,1);
 % initialization of the returned simulations.
 DynareOutput.endo_simul = NaN(DynareModel.endo_nbr,sample_size+1);
 if isempty(initial_conditions)
-    if isfield(DynareModel,'endo_histval')
+    if isfield(DynareModel,'endo_histval') && ~isempty(DynareModel.endo_histval)
         DynareOutput.endo_simul(:,1:DynareModel.maximum_lag) = DynareModel.endo_histval;
     else
         warning('simul_backward_nonlinear_model:: Initial condition is zero for all variables! If the model is nonlinear, the model simulation may fail with the default initialization')
@@ -97,11 +101,12 @@ else
 end
 Y = DynareOutput.endo_simul;
 
+
 % Simulations (call a Newton-like algorithm for each period).
 for it = 2:sample_size+1
     ylag = Y(iy1,it-1);                   % Set lagged variables.
     y = Y(:,it-1);                        % A good guess for the initial conditions is the previous values for the endogenous variables.
-    Y(:,it) = dynare_solve(model_dynamic_s, y, DynareOptions, model_dynamic, ylag, DynareOutput.exo_simul, DynareModel.params, DynareOutput.steady_state, it);
+    Y(:,it) = dynare_solve(model_dynamic_s, y, DynareOptions, model_dynamic, ylag, DynareOutput.exo_simul, DynareModel.params, DynareOutput.steady_state, it+(DynareModel.maximum_exo_lag-1));
 end
 
 DynareOutput.endo_simul = Y;