diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 1a2def69ddba42dde17ee2b896b2a3c141fc66d3..0c50e7dbcc549f0a644d7a56ad3514630a69948a 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -1350,6 +1350,8 @@ void
 DynamicModel::writeBlockDriverOutput(ostream &output, const string &basename,
                                      const vector<int> &state_var, bool estimation_present) const
 {
+  output << "M_.block_structure.time_recursive = " << boolalpha << time_recursive_block_decomposition << ";" << endl;
+
   for (int blk = 0; blk < static_cast<int>(blocks.size()); blk++)
     {
       int block_size = blocks[blk].size;
@@ -2959,6 +2961,11 @@ DynamicModel::computingPass(int derivsOrder, int paramsDerivsOrder, const eval_c
   if (paramsDerivsOrder > 0 && !no_tmp_terms)
     computeParamsDerivativesTemporaryTerms();
 
+  if (!block && (max_endo_lag == 0 || max_endo_lead == 0))
+    {
+      time_recursive_block_decomposition = true;
+      mfs = 3; // FIXME: remove this line when mfs=3 becomes the global default
+    }
   computingPassBlock(eval_context, no_tmp_terms);
   if (block_decomposed)
     computeBlockDynJacobianCols();
@@ -3077,7 +3084,9 @@ DynamicModel::determineBlockDerivativesType(int blk)
   map<tuple<int, int, int>, BlockDerivativeType> derivType;
   int size = blocks[blk].size;
   int nb_recursive = blocks[blk].getRecursiveSize();
-  for (int lag = -blocks[blk].max_lag; lag <= blocks[blk].max_lead; lag++)
+  for (int lag {time_recursive_block_decomposition ? 0 : -blocks[blk].max_lag};
+       lag <= (time_recursive_block_decomposition ? 0 : blocks[blk].max_lead);
+       lag++)
     for (int eq = 0; eq < size; eq++)
       {
         set<pair<int, int>> endos_and_lags;
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index 372ca70692b8cfd0af3dc4037be344a0a8983713..51060f45385f0dac00ba7d12d0deb9ee6ed45d40 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -166,6 +166,7 @@ ModelTree::ModelTree(const ModelTree &m) :
   eq_idx_orig2block{m.eq_idx_orig2block},
   endo_idx_orig2block{m.endo_idx_orig2block},
   block_decomposed{m.block_decomposed},
+  time_recursive_block_decomposition{m.time_recursive_block_decomposition},
   blocks{m.blocks},
   endo2block{m.endo2block},
   eq2block{m.eq2block},
@@ -210,6 +211,7 @@ ModelTree::operator=(const ModelTree &m)
   equation_type_and_normalized_equation.clear();
   blocks_derivatives.clear();
   block_decomposed = m.block_decomposed;
+  time_recursive_block_decomposition = m.time_recursive_block_decomposition;
   blocks = m.blocks;
   endo2block = m.endo2block;
   eq2block = m.eq2block;
@@ -408,7 +410,8 @@ ModelTree::computePrologueAndEpilogue()
       set<pair<int, int>> endos_and_lags;
       equations[i]->collectEndogenous(endos_and_lags);
       for (auto [endo, lag] : endos_and_lags)
-        IM[i * n + endo2eq[endo]] = true;
+        if (!time_recursive_block_decomposition || lag == 0)
+          IM[i * n + endo2eq[endo]] = true;
     }
 
   bool something_has_been_done;
@@ -662,7 +665,7 @@ ModelTree::computeBlockDecomposition(int prologue, int epilogue)
      For detecting dependencies between variables, use the symbolic adjacency
      matrix */
   VariableDependencyGraph G(nb_simvars);
-  for (const auto &[key, value] : computeSymbolicJacobian(false))
+  for (const auto &[key, value] : computeSymbolicJacobian(time_recursive_block_decomposition))
     {
       auto [eq, endo] = key;
       if (eq_idx_orig2block[eq] >= prologue
@@ -722,10 +725,11 @@ ModelTree::computeBlockDecomposition(int prologue, int epilogue)
      feedback set */
   for (int i = 0; i < nb_simvars; i++)
     if (equation_type_and_normalized_equation[eq_idx_block2orig[i+prologue]].first == EquationType::solve
-        || variable_lag_lead[endo_idx_block2orig[i+prologue]].first > 0
-        || variable_lag_lead[endo_idx_block2orig[i+prologue]].second > 0
-        || equation_lag_lead[eq_idx_block2orig[i+prologue]].first > 0
-        || equation_lag_lead[eq_idx_block2orig[i+prologue]].second > 0
+        || (!time_recursive_block_decomposition &&
+            (variable_lag_lead[endo_idx_block2orig[i+prologue]].first > 0
+             || variable_lag_lead[endo_idx_block2orig[i+prologue]].second > 0
+             || equation_lag_lead[eq_idx_block2orig[i+prologue]].first > 0
+             || equation_lag_lead[eq_idx_block2orig[i+prologue]].second > 0))
         || mfs == 0)
       add_edge(vertex(i, G), vertex(i, G), G);
 
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index c6312782224264f8f4ac455d0732da0895f82495..de892f53a4f8a34c97ceeff698e120dbfbb2d069 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -212,6 +212,25 @@ protected:
   // Whether block decomposition has been successfully computed
   bool block_decomposed {false};
 
+  /* Whether the block decomposition to compute is time-recursive (i.e. the
+     model can be simulated as a whole period-by-period).
+
+     If true, only contemporaneous occurrences of variables are considered when
+     computing the block structure; leads and lags are essentially treated as
+     exogenous, i.e. they are ignored. Such a decomposition only makes sense
+     for models that are purely backward/forward/static. When using the
+     resulting block decomposition to simulate the model, periods must be the
+     outer loop, and blocks the inner loop.
+
+     If false, then the full lead/lag structure is taken into account when
+     computing the block structure. This is the only option if there are both
+     leads and lags. When using the
+     resulting block decomposition to simulate the model, blocks must be the
+     outer loop, and periods the inner loop.
+
+     Of course, this setting does not make any difference on StaticModel. */
+  bool time_recursive_block_decomposition {false};
+
   // Stores various informations on the blocks
   vector<BlockInfo> blocks;