From 647a4b8060ca164442730b55f1adf7704b68c5c6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Adjemian=20=28Argos=29?=
 <stepan@adjemian.eu>
Date: Fri, 14 Mar 2025 11:18:08 +0100
Subject: [PATCH] Allow hybrid>2 in stochastic extended path.

Hybrid stochastic extended path with the constant correction in the
reduced form of the model approximated at order k>1 (perturbation
around the deterministic steady state).
---
 matlab/ep/ep_problem_1.m                        |  2 +-
 matlab/ep/ep_problem_2.m                        |  2 +-
 matlab/ep/extended_path_initialization.m        |  6 +++++-
 ...solve_stochastic_perfect_foresight_model_0.m | 17 ++++++++---------
 ...solve_stochastic_perfect_foresight_model_1.m |  4 ++++
 5 files changed, 19 insertions(+), 12 deletions(-)

diff --git a/matlab/ep/ep_problem_1.m b/matlab/ep/ep_problem_1.m
index 51c7fb1a80..35757c83af 100644
--- a/matlab/ep/ep_problem_1.m
+++ b/matlab/ep/ep_problem_1.m
@@ -43,7 +43,7 @@ for i = 1:order+1
         end
         if i <= order
             for k=1:nnodes
-                if hybrid_order==2 && i==order
+                if hybrid_order && i==order
                     z = [Y(i_cols_p,i_w_p);
                          Y(i_cols_s,j);
                          Y(i_cols_f,(j-1)*nnodes+k)+h_correction(i_hc)];
diff --git a/matlab/ep/ep_problem_2.m b/matlab/ep/ep_problem_2.m
index 2fda1838bc..26ba78895c 100644
--- a/matlab/ep/ep_problem_2.m
+++ b/matlab/ep/ep_problem_2.m
@@ -79,7 +79,7 @@ for i = 1:order+1
                 else
                     k1 = (nnodes-1)*(i-1)+k;
                 end
-                if hybrid_order == 2 && (k > 1 || i == order)
+                if hybrid_order && (k > 1 || i == order)
                     z = [Y(i_cols_p,1);
                          Y(i_cols_s,1);
                              Y(i_cols_f,k1)+h_correction(i_hc)];
diff --git a/matlab/ep/extended_path_initialization.m b/matlab/ep/extended_path_initialization.m
index 90828dc4b6..63197f3824 100644
--- a/matlab/ep/extended_path_initialization.m
+++ b/matlab/ep/extended_path_initialization.m
@@ -100,7 +100,11 @@ end
 
 % hybrid correction
 pfm.hybrid_order = ep.stochastic.hybrid_order;
-if pfm.hybrid_order
+if pfm.hybrid_order==1
+    warning('extended_path:: hybrid=1 is equivalent to hybrid=0 (option value must be an integer greater than 1 to be effective).')
+end
+
+if pfm.hybrid_order>1
     oo_.dr = set_state_space(oo_.dr, M_);
     options = options_;
     options.order = pfm.hybrid_order;
diff --git a/matlab/ep/solve_stochastic_perfect_foresight_model_0.m b/matlab/ep/solve_stochastic_perfect_foresight_model_0.m
index efdaf62472..9196ea79bc 100644
--- a/matlab/ep/solve_stochastic_perfect_foresight_model_0.m
+++ b/matlab/ep/solve_stochastic_perfect_foresight_model_0.m
@@ -62,18 +62,18 @@ if update_pfm_struct
     pfm.nodes = nodes;
     pfm.weights = weights;
 
-    hybrid_order = pfm.hybrid_order;
-    dr = pfm.dr;
-    if hybrid_order > 0
-        if hybrid_order == 2
-            h_correction = 0.5*dr.ghs2(dr.inv_order_var);
+    if pfm.hybrid_order>0
+        if pfm.hybrid_order==2
+            pfm.h_correction = 0.5*pfm.dr.ghs2(dr.inv_order_var);
+        elseif pfm.hybrid_order>2
+            pfm.h_correction = pfm.dr.g_0(pfm.dr.inv_order_var);
+        else
+            pfm.h_correction = 0;
         end
     else
-        h_correction = 0;
+        pfm.h_correction = 0;
     end
 
-    pfm.h_correction = h_correction;
-
     z = endo_simul(lead_lag_incidence_t(:)>0);
     [~, jacobian] = dynamic_model(z, exo_simul, pfm.params, pfm.steady_state, 2);
 
@@ -134,7 +134,6 @@ if update_pfm_struct
     pfm.order = order;
     pfm.world_nbr = world_nbr;
 
-    pfm.hybrid_order = hybrid_order;
     pfm.i_cols_1 = i_cols_1;
     pfm.i_cols_h = i_cols_j;
     pfm.icA = icA;
diff --git a/matlab/ep/solve_stochastic_perfect_foresight_model_1.m b/matlab/ep/solve_stochastic_perfect_foresight_model_1.m
index 19eb067380..d9f0deda76 100644
--- a/matlab/ep/solve_stochastic_perfect_foresight_model_1.m
+++ b/matlab/ep/solve_stochastic_perfect_foresight_model_1.m
@@ -60,6 +60,10 @@ if update_pfm_struct
     if pfm.hybrid_order > 0
         if pfm.hybrid_order == 2
             pfm.h_correction = 0.5*pfm.dr.ghs2(pfm.dr.inv_order_var);
+        elseif pfm.hybrid_order>2
+            pfm.h_correction = pfm.dr.g_0(pfm.dr.inv_order_var);
+        else
+            pfm.h_correction = 0;
         end
     else
         pfm.h_correction = 0;
-- 
GitLab