From 23d49470c80b950f3beace3458ca45dd6b260960 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Fri, 4 Apr 2025 15:17:00 +0200
Subject: [PATCH] =?UTF-8?q?Make=20perfect=5Fforesight=5Fcontrolled=5Fpaths?=
 =?UTF-8?q?=20compatible=20with=20the=20=E2=80=9Clinear=E2=80=9D=20option?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 doc/manual/source/the-model-file.rst          |  4 +-
 .../perfect_foresight_solver_core.m           |  5 +-
 matlab/perfect-foresight-models/sim1_linear.m | 41 +++++++++++--
 meson.build                                   |  1 +
 .../linear_state_space_arma_condpaths.mod     | 60 +++++++++++++++++++
 5 files changed, 101 insertions(+), 10 deletions(-)
 create mode 100644 tests/deterministic_simulations/linear_state_space_arma_condpaths.mod

diff --git a/doc/manual/source/the-model-file.rst b/doc/manual/source/the-model-file.rst
index c47208e9f0..86c0533844 100644
--- a/doc/manual/source/the-model-file.rst
+++ b/doc/manual/source/the-model-file.rst
@@ -4672,8 +4672,8 @@ that have been left free.
     The ``perfect_foresight_controlled_paths`` block requires that the
     :opt:`stack_solve_algo <stack_solve_algo = INTEGER>` option be equal to
     either ``0``, ``1``, ``2``, ``3``, ``6`` or ``7``, and is incompatible with
-    the :opt:`block`, :opt:`linear` and :opt:`bytecode` options of the
-    :bck:`model` block and :comm:`model_options` command.
+    the :opt:`block` and :opt:`bytecode` options of the :bck:`model` block and
+    :comm:`model_options` command.
 
     *Options*
 
diff --git a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
index c1977137cc..879961b8e9 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_solver_core.m
@@ -90,9 +90,6 @@ if ~isempty(controlled_paths_by_period)
     if options_.block
         error('perfect_foresight_controlled_paths is not available with the block option')
     end
-    if options_.linear
-        error('perfect_foresight_controlled_paths is not available with the linear option')
-    end
     if M_.maximum_endo_lead == 0 || M_.maximum_endo_lag == 0
         error('perfect_foresight_controlled_paths is not available with purely backward, purely forward or static models')
     end
@@ -147,7 +144,7 @@ else
             switch options_.stack_solve_algo
               case {0 2 3}
                 if options_.linear_approximation
-                    [y, success, maxerror] = sim1_linear(y, exo_simul, steady_state, exo_steady_state, M_, options_);
+                    [y, success, maxerror, exo_simul] = sim1_linear(y, exo_simul, steady_state, exo_steady_state, controlled_paths_by_period, M_, options_);
                 else
                     [y, success, maxerror, iter, exo_simul] = sim1(y, exo_simul, steady_state, controlled_paths_by_period, M_, options_);
                 end
diff --git a/matlab/perfect-foresight-models/sim1_linear.m b/matlab/perfect-foresight-models/sim1_linear.m
index 30c11700cf..6524cf22f2 100644
--- a/matlab/perfect-foresight-models/sim1_linear.m
+++ b/matlab/perfect-foresight-models/sim1_linear.m
@@ -1,12 +1,13 @@
-function [endogenousvariables, success, ERR] = sim1_linear(endogenousvariables, exogenousvariables, steadystate_y, steadystate_x, M_, options_)
+function [endogenousvariables, success, ERR, exogenousvariables] = sim1_linear(endogenousvariables, exogenousvariables, steadystate_y, steadystate_x, controlled_paths_by_period, M_, options_)
 % [endogenousvariables, success, ERR] = sim1_linear(endogenousvariables, exogenousvariables, steadystate_y, steadystate_x, M_, options_)
 % Solves a linear approximation of a perfect foresight model using sparse matrix.
 %
 % INPUTS
-% - endogenousvariables [double] N*T array, paths for the endogenous variables (initial guess).
+% - endogenousvariables [double] N*T array, paths for the endogenous variables (initial condition + initial guess + terminal condition).
 % - exogenousvariables  [double] T*M array, paths for the exogenous variables.
 % - steadystate_y       [double] N*1 array, steady state for the endogenous variables.
 % - steadystate_x       [double] M*1 array, steady state for the exogenous variables.
+% - controlled_paths_by_period [struct] data from perfect_foresight_controlled_paths block
 % - M_                  [struct] contains a description of the model.
 % - options_            [struct] contains various options.
 %
@@ -14,6 +15,8 @@ function [endogenousvariables, success, ERR] = sim1_linear(endogenousvariables,
 % - endogenousvariables [double] N*T array, paths for the endogenous variables (solution of the perfect foresight model).
 % - success             [logical] Whether a solution was found
 % - ERR                 [double] ∞-norm of the residual
+% - exogenousvariables  [double] T*M array, paths for the exogenous variables
+%                                (may be modified if perfect_foresight_controlled_paths present)
 %
 % NOTATIONS
 % - N is the number of endogenous variables.
@@ -39,7 +42,7 @@ function [endogenousvariables, success, ERR] = sim1_linear(endogenousvariables,
 % to center the variables around the deterministic steady state to solve the
 % perfect foresight model.
 
-% Copyright © 2015-2024 Dynare Team
+% Copyright © 2015-2025 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -67,6 +70,19 @@ periods = get_simulation_periods(options_);
 
 params = M_.params;
 
+if ~isempty(controlled_paths_by_period)
+    for p = 1:periods
+        if isempty(controlled_paths_by_period(p).exogenize_id)
+            continue
+        end
+        endogenousvariables(controlled_paths_by_period(p).exogenize_id, p+M_.maximum_lag) = controlled_paths_by_period(p).values;
+    end
+
+    if options_.debug
+        error('Debugging not available with perfect_foresight_controlled_paths')
+    end
+end
+
 if verbose
     skipline()
     printline(80)
@@ -125,6 +141,10 @@ h2 = clock;
 
 [res, A] = linear_perfect_foresight_problem(Y, jacobian, y0, yT, exogenousvariables, params, steadystate_y, maximum_lag, periods, ny);
 
+if ~isempty(controlled_paths_by_period)
+    A = controlled_paths_substitute_stacked_jacobian(A, Y, y0, yT, exogenousvariables, steadystate_y, controlled_paths_by_period, M_);
+end
+
 % Evaluation of the maximum residual at the initial guess (steady state for the endogenous variables).
 err = max(abs(res));
 
@@ -141,7 +161,8 @@ end
 
 % Try to update the vector of endogenous variables.
 try
-    Y =  Y - A\res;
+    dY = -A\res;
+    Y = Y + dY;
 catch
     % Normally, because the model is linear, the solution of the perfect foresight model should
     % be obtained in one Newton step. This is not the case if the model is singular.
@@ -154,6 +175,18 @@ catch
     return
 end
 
+if ~isempty(controlled_paths_by_period)
+    for p = 1:periods
+        endogenize_id = controlled_paths_by_period(p).endogenize_id;
+        exogenize_id = controlled_paths_by_period(p).exogenize_id;
+        if isempty(endogenize_id)
+            continue
+        end
+        Y(exogenize_id+(p-1)*M_.endo_nbr) = controlled_paths_by_period(p).values;
+        exogenousvariables(p+M_.maximum_lag,endogenize_id) = exogenousvariables(p+M_.maximum_lag,endogenize_id) + dY(exogenize_id+(p-1)*M_.endo_nbr)';
+    end
+end
+
 YY = reshape([y0; Y; yT], ny, periods+2);
 
 i_rows = (1:ny)';
diff --git a/meson.build b/meson.build
index f52e708e87..351922ff61 100644
--- a/meson.build
+++ b/meson.build
@@ -1385,6 +1385,7 @@ mod_and_m_tests = [
   { 'test' : [ 'deterministic_simulations/Irreversible_investment.mod',
                'deterministic_simulations/Irreversible_investment_leadlag.mod' ] },
   { 'test' : [ 'deterministic_simulations/linear_state_space_arma.mod' ] },
+  { 'test' : [ 'deterministic_simulations/linear_state_space_arma_condpaths.mod' ] },
   { 'test' : [ 'deterministic_simulations/purely_forward/ar1.mod' ] },
   { 'test' : [ 'deterministic_simulations/purely_forward/nk.mod' ] },
   { 'test' : [ 'deterministic_simulations/purely_backward/ar1.mod' ] },
diff --git a/tests/deterministic_simulations/linear_state_space_arma_condpaths.mod b/tests/deterministic_simulations/linear_state_space_arma_condpaths.mod
new file mode 100644
index 0000000000..eb408209f6
--- /dev/null
+++ b/tests/deterministic_simulations/linear_state_space_arma_condpaths.mod
@@ -0,0 +1,60 @@
+% Tests perfect simulation of linear model (sim1_linear.m) with controlled paths
+
+var x
+    y
+    z;
+
+varexo u
+       v;
+
+parameters a1 a2 a3 a4
+	   b1 b2 b3
+	   c1;
+
+a1 =  .50;
+a2 =  .00;
+a3 =  .70;
+a4 =  .40;
+b1 =  .90;
+b2 =  .00;
+b3 =  .80;
+c1 =  .95;
+
+
+model(linear);
+  y = a1*x(-1) + a2*x(+1) + a3*z + a4*y(-1);
+  z = b1*z(-1) + b2*z(+1) + b3*x + u;
+  x = c1*x(-1) + v +v(-1)+v(+1);
+end;
+
+initval;
+  y=-1;
+  x=-1;
+  z=-1;
+end;
+
+endval;
+  y=0;
+  x=0;
+  z=0;
+end;
+
+perfect_foresight_controlled_paths;
+  exogenize z;
+  periods 50;
+  values 0.5;
+  endogenize u;
+end;
+
+steady;
+
+perfect_foresight_setup(periods=100);
+perfect_foresight_solver;
+
+if ~oo_.deterministic_simulation.status
+   error('Perfect foresight simulation failed')
+end
+
+if ~(oo_.endo_simul(3, 51) == 0.5 && isequal(find(oo_.exo_simul ~= 0), 51))
+   error('Linear perfect foresight simulation with controlled paths failed')
+end
-- 
GitLab