diff --git a/src/ComputingTasks.cc b/src/ComputingTasks.cc
index f980691b..12f198ff 100644
--- a/src/ComputingTasks.cc
+++ b/src/ComputingTasks.cc
@@ -2111,10 +2111,8 @@ ModelComparisonStatement::writeJsonOutput(ostream &output) const
   output << "}";
 }
 
-PlannerObjectiveStatement::PlannerObjectiveStatement(SymbolTable &symbol_table,
-                                                     NumericalConstants &num_constants,
-                                                     ExternalFunctionsTable &external_functions_table) :
-  model_tree{symbol_table, num_constants, external_functions_table}
+PlannerObjectiveStatement::PlannerObjectiveStatement(const StaticModel &model_tree_arg) :
+  model_tree {model_tree_arg}
 {
 }
 
@@ -2132,8 +2130,8 @@ PlannerObjectiveStatement::checkPass(ModFileStructure &mod_file_struct, WarningC
   mod_file_struct.planner_objective_present = true;
 }
 
-StaticModel &
-PlannerObjectiveStatement::getPlannerObjective()
+const StaticModel &
+PlannerObjectiveStatement::getPlannerObjective() const
 {
   return model_tree;
 }
diff --git a/src/ComputingTasks.hh b/src/ComputingTasks.hh
index 595cd981..8d7982d7 100644
--- a/src/ComputingTasks.hh
+++ b/src/ComputingTasks.hh
@@ -523,9 +523,7 @@ private:
   StaticModel model_tree;
   bool computing_pass_called{false};
 public:
-  PlannerObjectiveStatement(SymbolTable &symbol_table,
-                            NumericalConstants &num_constants,
-                            ExternalFunctionsTable &external_functions_table);
+  PlannerObjectiveStatement(const StaticModel &model_tree_arg);
   /*! \todo check there are only endogenous variables at the current period in the objective
     (no exogenous, no lead/lag) */
   void checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings) override;
@@ -534,7 +532,7 @@ public:
   void writeOutput(ostream &output, const string &basename, bool minimal_workspace) const override;
   void writeJsonOutput(ostream &output) const override;
   //! Return a reference the Planner Objective model tree
-  StaticModel &getPlannerObjective();
+  const StaticModel &getPlannerObjective() const;
 };
 
 class BVARDensityStatement : public Statement
diff --git a/src/DataTree.cc b/src/DataTree.cc
index 3cf53845..f63b9ede 100644
--- a/src/DataTree.cc
+++ b/src/DataTree.cc
@@ -26,12 +26,8 @@
 
 #include "DataTree.hh"
 
-DataTree::DataTree(SymbolTable &symbol_table_arg,
-                   NumericalConstants &num_constants_arg,
-                   ExternalFunctionsTable &external_functions_table_arg) :
-  symbol_table{symbol_table_arg},
-  num_constants{num_constants_arg},
-  external_functions_table{external_functions_table_arg}
+void
+DataTree::initConstants()
 {
   Zero = AddNonNegativeConstant("0");
   One = AddNonNegativeConstant("1");
@@ -46,8 +42,76 @@ DataTree::DataTree(SymbolTable &symbol_table_arg,
   Pi = AddNonNegativeConstant("3.141592653589793");
 }
 
+DataTree::DataTree(SymbolTable &symbol_table_arg,
+                   NumericalConstants &num_constants_arg,
+                   ExternalFunctionsTable &external_functions_table_arg,
+                   bool is_dynamic_arg) :
+  symbol_table{symbol_table_arg},
+  num_constants{num_constants_arg},
+  external_functions_table{external_functions_table_arg},
+  is_dynamic {is_dynamic_arg}
+{
+  initConstants();
+}
+
 DataTree::~DataTree() = default;
 
+DataTree::DataTree(const DataTree &d) :
+  symbol_table {d.symbol_table},
+  num_constants {d.num_constants},
+  external_functions_table {d.external_functions_table},
+  is_dynamic {d.is_dynamic},
+  local_variables_vector {d.local_variables_vector}
+{
+  // Constants must be initialized first because they are used in some Add* methods
+  initConstants();
+
+  for (const auto & it : d.node_list)
+    it->cloneDynamic(*this);
+
+  assert(node_list.size() == d.node_list.size());
+
+  for (const auto & it : d.local_variables_table)
+    local_variables_table[it.first] = it.second->cloneDynamic(*this);
+}
+
+DataTree &
+DataTree::operator=(const DataTree &d)
+{
+  assert (&symbol_table == &d.symbol_table);
+  assert (&num_constants == &d.num_constants);
+  assert (&external_functions_table == &d.external_functions_table);
+  assert (is_dynamic == d.is_dynamic);
+
+  num_const_node_map.clear();
+  variable_node_map.clear();
+  unary_op_node_map.clear();
+  binary_op_node_map.clear();
+  trinary_op_node_map.clear();
+  external_function_node_map.clear();
+  var_expectation_node_map.clear();
+  pac_expectation_node_map.clear();
+  first_deriv_external_function_node_map.clear();
+  second_deriv_external_function_node_map.clear();
+
+  node_list.clear();
+
+  // Constants must be initialized first because they are used in some Add* methods
+  initConstants();
+
+  for (const auto & it : d.node_list)
+    it->cloneDynamic(*this);
+
+  assert(node_list.size() == d.node_list.size());
+
+  local_variables_vector = d.local_variables_vector;
+
+  for (const auto & it : d.local_variables_table)
+    local_variables_table[it.first] = it.second->cloneDynamic(*this);
+
+  return *this;
+}
+
 expr_t
 DataTree::AddNonNegativeConstant(const string &value)
 {
@@ -65,8 +129,14 @@ DataTree::AddNonNegativeConstant(const string &value)
 }
 
 VariableNode *
-DataTree::AddVariableInternal(int symb_id, int lag)
+DataTree::AddVariable(int symb_id, int lag)
 {
+  if (lag != 0 && !is_dynamic)
+    {
+      cerr << "Leads/lags not authorized in this DataTree" << endl;
+      exit(EXIT_FAILURE);
+    }
+
   auto it = variable_node_map.find({ symb_id, lag });
   if (it != variable_node_map.end())
     return it->second;
@@ -87,13 +157,6 @@ DataTree::ParamUsedWithLeadLagInternal() const
   return false;
 }
 
-VariableNode *
-DataTree::AddVariable(int symb_id, int lag)
-{
-  assert(lag == 0);
-  return AddVariableInternal(symb_id, lag);
-}
-
 expr_t
 DataTree::AddPlus(expr_t iArg1, expr_t iArg2)
 {
diff --git a/src/DataTree.hh b/src/DataTree.hh
index cea53af2..cbfaca81 100644
--- a/src/DataTree.hh
+++ b/src/DataTree.hh
@@ -45,6 +45,8 @@ public:
   NumericalConstants &num_constants;
   //! A reference to the external functions table
   ExternalFunctionsTable &external_functions_table;
+  //! Is it possible to use leads/lags on variable nodes?
+  const bool is_dynamic;
 
 protected:
   //! num_constant_id -> NumConstNode
@@ -92,9 +94,6 @@ protected:
   //! Stores the order of appearance of local variables in the model block. Needed following change in #563
   vector<int> local_variables_vector;
 
-  //! Internal implementation of AddVariable(), without the check on the lag
-  VariableNode *AddVariableInternal(int symb_id, int lag);
-
   //! Internal implementation of ParamUsedWithLeadLag()
   bool ParamUsedWithLeadLagInternal() const;
 
@@ -115,17 +114,21 @@ private:
   inline expr_t AddBinaryOp(expr_t arg1, BinaryOpcode op_code, expr_t arg2, int powerDerivOrder = 0);
   inline expr_t AddTrinaryOp(expr_t arg1, TrinaryOpcode op_code, expr_t arg2, expr_t arg3);
 
+  //! Initializes the predefined constants, used only from the constructors
+  void initConstants();
+
 public:
   DataTree(SymbolTable &symbol_table_arg,
            NumericalConstants &num_constants_arg,
-           ExternalFunctionsTable &external_functions_table_arg);
+           ExternalFunctionsTable &external_functions_table_arg,
+           bool is_static_args = false);
 
   virtual
   ~DataTree();
 
-  DataTree(const DataTree &) = delete;
+  DataTree(const DataTree &d);
   DataTree(DataTree &&) = delete;
-  DataTree & operator=(const DataTree &) = delete;
+  DataTree & operator=(const DataTree &d);
   DataTree & operator=(DataTree &&) = delete;
 
   //! Some predefined constants
@@ -149,8 +152,7 @@ public:
   //! Adds a non-negative numerical constant (possibly Inf or NaN)
   expr_t AddNonNegativeConstant(const string &value);
   //! Adds a variable
-  /*! The default implementation of the method refuses any lag != 0 */
-  virtual VariableNode *AddVariable(int symb_id, int lag = 0);
+  VariableNode *AddVariable(int symb_id, int lag = 0);
   //! Adds "arg1+arg2" to model tree
   expr_t AddPlus(expr_t iArg1, expr_t iArg2);
   //! Adds "arg1-arg2" to model tree
diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 2153a3bc..1cf02763 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -31,21 +31,217 @@
 
 #include "DynamicModel.hh"
 
+void
+DynamicModel::copyHelper(const DynamicModel &m)
+{
+  auto f = [this](const ExprNode *e) { return e->cloneDynamic(*this); };
+
+  for (const auto &it : m.static_only_equations)
+    static_only_equations.push_back(dynamic_cast<BinaryOpNode *>(f(it)));
+
+  auto convert_vector_tt = [this,f](vector<temporary_terms_t> vtt)
+    {
+      vector<temporary_terms_t> vtt2;
+      for (const auto &tt : vtt)
+        {
+          temporary_terms_t tt2;
+          for (const auto &it : tt)
+            tt2.insert(f(it));
+          vtt2.push_back(tt2);
+        }
+      return vtt2;
+    };
+
+  for (const auto &it : m.v_temporary_terms)
+    v_temporary_terms.push_back(convert_vector_tt(it));
+
+  for (const auto &it : m.first_chain_rule_derivatives)
+    first_chain_rule_derivatives[it.first] = f(it.second);
+
+  for (const auto &it : m.equation_type_and_normalized_equation)
+    equation_type_and_normalized_equation.push_back(make_pair(it.first, f(it.second)));
+
+  for (const auto &it : m.blocks_derivatives)
+    {
+      block_derivatives_equation_variable_laglead_nodeid_t v;
+      for (const auto &it2 : it)
+        v.push_back(make_pair(it2.first, make_pair(it2.second.first, f(it2.second.second))));
+      blocks_derivatives.push_back(v);
+    }
+
+  for (const auto &it : m.dynamic_jacobian)
+    dynamic_jacobian[it.first] = f(it.second);
+
+  auto convert_derivative_t = [this,f](derivative_t dt)
+    {
+      derivative_t dt2;
+      for (const auto &it : dt)
+        dt2[it.first] = f(it.second);
+      return dt2;
+    };
+  for (const auto &it : m.derivative_endo)
+    derivative_endo.push_back(convert_derivative_t(it));
+  for (const auto &it : m.derivative_other_endo)
+    derivative_other_endo.push_back(convert_derivative_t(it));
+  for (const auto &it : m.derivative_exo)
+    derivative_exo.push_back(convert_derivative_t(it));
+  for (const auto &it : m.derivative_exo_det)
+    derivative_exo_det.push_back(convert_derivative_t(it));
+
+  for (const auto &it : m.pac_expectation_info)
+    pac_expectation_info.insert(dynamic_cast<const PacExpectationNode *>(f(it)));
+}
+
+
 DynamicModel::DynamicModel(SymbolTable &symbol_table_arg,
                            NumericalConstants &num_constants_arg,
                            ExternalFunctionsTable &external_functions_table_arg,
                            TrendComponentModelTable &trend_component_model_table_arg,
                            VarModelTable &var_model_table_arg) :
-  ModelTree{symbol_table_arg, num_constants_arg, external_functions_table_arg},
+  ModelTree {symbol_table_arg, num_constants_arg, external_functions_table_arg, true},
   trend_component_model_table{trend_component_model_table_arg},
   var_model_table{var_model_table_arg}
 {
 }
 
-VariableNode *
-DynamicModel::AddVariable(int symb_id, int lag)
+DynamicModel::DynamicModel(const DynamicModel &m) :
+  ModelTree {m},
+  trend_component_model_table {m.trend_component_model_table},
+  var_model_table {m.var_model_table},
+  static_only_equations_lineno {m.static_only_equations_lineno},
+  static_only_equations_equation_tags {m.static_only_equations_equation_tags},
+  deriv_id_table {m.deriv_id_table},
+  inv_deriv_id_table {m.inv_deriv_id_table},
+  dyn_jacobian_cols_table {m.dyn_jacobian_cols_table},
+  max_lag {m.max_lag},
+  max_lead {m.max_lead},
+  max_endo_lag {m.max_endo_lag},
+  max_endo_lead {m.max_endo_lead},
+  max_exo_lag {m.max_exo_lag},
+  max_exo_lead {m.max_exo_lead},
+  max_exo_det_lag {m.max_exo_det_lag},
+  max_exo_det_lead {m.max_exo_det_lead},
+  max_lag_orig {m.max_lag_orig},
+  max_lead_orig {m.max_lead_orig},
+  max_endo_lag_orig {m.max_endo_lag_orig},
+  max_endo_lead_orig {m.max_endo_lead_orig},
+  max_exo_lag_orig {m.max_exo_lag_orig},
+  max_exo_lead_orig {m.max_exo_lead_orig},
+  max_exo_det_lag_orig {m.max_exo_det_lag_orig},
+  max_exo_det_lead_orig {m.max_exo_det_lead_orig},
+  xrefs {m.xrefs},
+  xref_param  {m.xref_param},
+  xref_endo {m.xref_endo},
+  xref_exo {m.xref_exo},
+  xref_exo_det {m.xref_exo_det},
+  nonzero_hessian_eqs {m.nonzero_hessian_eqs},
+  v_temporary_terms_inuse {m.v_temporary_terms_inuse},
+  map_idx {m.map_idx},
+  global_temporary_terms {m.global_temporary_terms},
+  block_type_firstequation_size_mfs {m.block_type_firstequation_size_mfs},
+  blocks_linear {m.blocks_linear},
+  other_endo_block {m.other_endo_block},
+  exo_block {m.exo_block},
+  exo_det_block {m.exo_det_block},
+  block_var_exo {m.block_var_exo},
+  block_exo_index {m.block_exo_index},
+  block_det_exo_index {m.block_det_exo_index},
+  block_other_endo_index {m.block_other_endo_index},
+  block_col_type {m.block_col_type},
+  variable_block_lead_lag {m.variable_block_lead_lag},
+  equation_block {m.equation_block},
+  var_expectation_functions_to_write {m.var_expectation_functions_to_write},
+  endo_max_leadlag_block {m.endo_max_leadlag_block},
+  other_endo_max_leadlag_block {m.other_endo_max_leadlag_block},
+  exo_max_leadlag_block {m.exo_max_leadlag_block},
+  exo_det_max_leadlag_block {m.exo_det_max_leadlag_block},
+  max_leadlag_block {m.max_leadlag_block}
 {
-  return AddVariableInternal(symb_id, lag);
+  copyHelper(m);
+}
+
+DynamicModel &
+DynamicModel::operator=(const DynamicModel &m)
+{
+  ModelTree::operator=(m);
+
+  assert(&trend_component_model_table == &m.trend_component_model_table);
+  assert(&var_model_table == &m.var_model_table);
+
+  static_only_equations_lineno = m.static_only_equations_lineno;
+  static_only_equations_equation_tags = m.static_only_equations_equation_tags;
+  deriv_id_table = m.deriv_id_table;
+  inv_deriv_id_table = m.inv_deriv_id_table;
+  dyn_jacobian_cols_table = m.dyn_jacobian_cols_table;
+  max_lag = m.max_lag;
+  max_lead = m.max_lead;
+  max_endo_lag = m.max_endo_lag;
+  max_endo_lead = m.max_endo_lead;
+  max_exo_lag = m.max_exo_lag;
+  max_exo_lead = m.max_exo_lead;
+  max_exo_det_lag = m.max_exo_det_lag;
+  max_exo_det_lead = m.max_exo_det_lead;
+  max_lag_orig = m.max_lag_orig;
+  max_lead_orig = m.max_lead_orig;
+  max_endo_lag_orig = m.max_endo_lag_orig;
+  max_endo_lead_orig = m.max_endo_lead_orig;
+  max_exo_lag_orig = m.max_exo_lag_orig;
+  max_exo_lead_orig = m.max_exo_lead_orig;
+  max_exo_det_lag_orig = m.max_exo_det_lag_orig;
+  max_exo_det_lead_orig = m.max_exo_det_lead_orig;
+  xrefs = m.xrefs;
+  xref_param  = m.xref_param;
+  xref_endo = m.xref_endo;
+  xref_exo = m.xref_exo;
+  xref_exo_det = m.xref_exo_det;
+  nonzero_hessian_eqs = m.nonzero_hessian_eqs;
+
+  v_temporary_terms.clear();
+
+  v_temporary_terms_inuse = m.v_temporary_terms_inuse;
+
+  first_chain_rule_derivatives.clear();
+
+  map_idx = m.map_idx;
+  global_temporary_terms = m.global_temporary_terms;
+
+  equation_type_and_normalized_equation.clear();
+
+  block_type_firstequation_size_mfs = m.block_type_firstequation_size_mfs;
+
+  blocks_derivatives.clear();
+  dynamic_jacobian.clear();
+
+  blocks_linear = m.blocks_linear;
+
+  derivative_endo.clear();
+  derivative_other_endo.clear();
+  derivative_exo.clear();
+  derivative_exo_det.clear();
+
+  other_endo_block = m.other_endo_block;
+  exo_block = m.exo_block;
+  exo_det_block = m.exo_det_block;
+  block_var_exo = m.block_var_exo;
+  block_exo_index = m.block_exo_index;
+  block_det_exo_index = m.block_det_exo_index;
+  block_other_endo_index = m.block_other_endo_index;
+  block_col_type = m.block_col_type;
+  variable_block_lead_lag = m.variable_block_lead_lag;
+  equation_block = m.equation_block;
+  var_expectation_functions_to_write = m.var_expectation_functions_to_write;
+
+  pac_expectation_info.clear();
+
+  endo_max_leadlag_block = m.endo_max_leadlag_block;
+  other_endo_max_leadlag_block = m.other_endo_max_leadlag_block;
+  exo_max_leadlag_block = m.exo_max_leadlag_block;
+  exo_det_max_leadlag_block = m.exo_det_max_leadlag_block;
+  max_leadlag_block = m.max_leadlag_block;
+
+  copyHelper(m);
+
+  return *this;
 }
 
 void
@@ -4827,40 +5023,6 @@ DynamicModel::writeAuxVarRecursiveDefinitions(ostream &output, ExprNodeOutputTyp
     }
 }
 
-void
-DynamicModel::cloneDynamic(DynamicModel &dynamic_model) const
-{
-  /* Ensure that we are using the same symbol table, because at many places we manipulate
-     symbol IDs rather than strings */
-  assert(&symbol_table == &dynamic_model.symbol_table);
-
-  // Convert model local variables (need to be done first)
-  for (int it : local_variables_vector)
-    dynamic_model.AddLocalVariable(it, local_variables_table.find(it)->second->cloneDynamic(dynamic_model));
-
-  // Convert equations
-  for (size_t i = 0; i < equations.size(); i++)
-    {
-      vector<pair<string, string>> eq_tags;
-      for (const auto & equation_tag : equation_tags)
-        if (equation_tag.first == (int)i)
-          eq_tags.push_back(equation_tag.second);
-      dynamic_model.addEquation(equations[i]->cloneDynamic(dynamic_model), equations_lineno[i], eq_tags);
-    }
-
-  // Convert auxiliary equations
-  for (auto aux_equation : aux_equations)
-    dynamic_model.addAuxEquation(aux_equation->cloneDynamic(dynamic_model));
-
-  // Convert static_only equations
-  for (size_t i = 0; i < static_only_equations.size(); i++)
-    dynamic_model.addStaticOnlyEquation(static_only_equations[i]->cloneDynamic(dynamic_model),
-                                        static_only_equations_lineno[i],
-                                        static_only_equations_equation_tags[i]);
-
-  dynamic_model.setLeadsLagsOrig();
-}
-
 void
 DynamicModel::replaceMyEquations(DynamicModel &dynamic_model) const
 {
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 8ff5c215..75a8096e 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -254,15 +254,22 @@ private:
 
   void findPacExpectationEquationNumbers(vector<int> &eqnumber) const;
 
+  //! Internal helper for the copy constructor and assignment operator
+  /*! Copies all the structures that contain ExprNode*, by the converting the
+      pointers into their equivalent in the new tree */
+  void copyHelper(const DynamicModel &m);
+
 public:
   DynamicModel(SymbolTable &symbol_table_arg,
                NumericalConstants &num_constants_arg,
                ExternalFunctionsTable &external_functions_table_arg,
                TrendComponentModelTable &trend_component_model_table_arg,
                VarModelTable &var_model_table_arg);
-  //! Adds a variable node
-  /*! This implementation allows for non-zero lag */
-  VariableNode *AddVariable(int symb_id, int lag = 0) override;
+
+  DynamicModel(const DynamicModel &m);
+  DynamicModel(DynamicModel &&) = delete;
+  DynamicModel & operator=(const DynamicModel &m);
+  DynamicModel & operator=(DynamicModel &&) = delete;
 
   //! Compute cross references
   void computeXrefs();
@@ -374,10 +381,6 @@ public:
   //! Set the max leads/lags of the original model
   void setLeadsLagsOrig();
 
-  //! Copies a dynamic model (only the equations)
-  /*! It assumes that the dynamic model given in argument has just been allocated */
-  void cloneDynamic(DynamicModel &dynamic_model) const;
-
   //! Replaces model equations with derivatives of Lagrangian w.r.t. endogenous
   void computeRamseyPolicyFOCs(const StaticModel &static_model, const bool nopreprocessoroutput);
   //! Replaces the model equations in dynamic_model with those in this model
diff --git a/src/ModFile.cc b/src/ModFile.cc
index 03f99f07..dc6a4841 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -367,7 +367,7 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
   // - except adl and diff which we always want expanded
   dynamic_model.substituteAdl();
   dynamic_model.setLeadsLagsOrig();
-  dynamic_model.cloneDynamic(original_model);
+  original_model = dynamic_model;
 
   if (nostrict)
     {
@@ -457,7 +457,7 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
   if (nonstationary_variables)
     {
       dynamic_model.detrendEquations();
-      dynamic_model.cloneDynamic(trend_dynamic_model);
+      trend_dynamic_model = dynamic_model;
       dynamic_model.removeTrendVariableFromEquations();
     }
 
@@ -485,8 +485,8 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
         we have to call computeDerivIDs (in computeRamseyPolicyFOCs and computingPass)
       */
       if (linear)
-        dynamic_model.cloneDynamic(orig_ramsey_dynamic_model);
-      dynamic_model.cloneDynamic(ramsey_FOC_equations_dynamic_model);
+        orig_ramsey_dynamic_model = dynamic_model;
+      ramsey_FOC_equations_dynamic_model = dynamic_model;
       ramsey_FOC_equations_dynamic_model.computeRamseyPolicyFOCs(planner_objective, nopreprocessoroutput);
       ramsey_FOC_equations_dynamic_model.replaceMyEquations(dynamic_model);
       mod_file_struct.ramsey_eq_nbr = dynamic_model.equation_number() - mod_file_struct.orig_eq_nbr;
@@ -696,7 +696,7 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_deri
       dynamic_model.toStatic(static_model);
 	  if (linear_decomposition)
         {
-          dynamic_model.cloneDynamic(non_linear_equations_dynamic_model);
+          non_linear_equations_dynamic_model = dynamic_model;
           non_linear_equations_dynamic_model.set_cutoff_to_zero();
           non_linear_equations_dynamic_model.computingPass(true, false, false, 0, global_eval_context, no_tmp_terms, block, use_dll, byte_code, nopreprocessoroutput, linear_decomposition);
         }
diff --git a/src/ModelEquationBlock.cc b/src/ModelEquationBlock.cc
index 2f0b859c..7af50a3a 100644
--- a/src/ModelEquationBlock.cc
+++ b/src/ModelEquationBlock.cc
@@ -31,6 +31,28 @@ SteadyStateModel::SteadyStateModel(SymbolTable &symbol_table_arg,
 {
 }
 
+SteadyStateModel::SteadyStateModel(const SteadyStateModel &m) :
+  DataTree {m},
+  static_model {m.static_model}
+{
+  for (const auto &it : m.def_table)
+    def_table.push_back(make_pair(it.first, it.second->cloneDynamic(*this)));
+}
+
+SteadyStateModel &
+SteadyStateModel::operator=(const SteadyStateModel &m)
+{
+  DataTree::operator=(m);
+
+  assert(&static_model == &m.static_model);
+
+  def_table.clear();
+  for (const auto &it : m.def_table)
+    def_table.push_back(make_pair(it.first, it.second->cloneDynamic(*this)));
+
+  return *this;
+}
+
 void
 SteadyStateModel::addDefinition(int symb_id, expr_t expr)
 {
@@ -277,6 +299,30 @@ Epilogue::Epilogue(SymbolTable &symbol_table_arg,
 {
 }
 
+Epilogue::Epilogue(const Epilogue &m) :
+  DynamicModel {m},
+  endogs {m.endogs},
+  exogs {m.exogs}
+{
+  for (const auto &it : m.def_table)
+    def_table.push_back(make_pair(it.first, it.second->cloneDynamic(*this)));
+}
+
+Epilogue &
+Epilogue::operator=(const Epilogue &m)
+{
+  DynamicModel::operator=(m);
+
+  endogs = m.endogs;
+  exogs = m.exogs;
+
+  def_table.clear();
+  for (const auto &it : m.def_table)
+    def_table.push_back(make_pair(it.first, it.second->cloneDynamic(*this)));
+
+  return *this;
+}
+
 void
 Epilogue::addDefinition(int symb_id, expr_t expr)
 {
diff --git a/src/ModelEquationBlock.hh b/src/ModelEquationBlock.hh
index 67a4c878..4391c8ad 100644
--- a/src/ModelEquationBlock.hh
+++ b/src/ModelEquationBlock.hh
@@ -40,6 +40,12 @@ public:
                    NumericalConstants &num_constants_arg,
                    ExternalFunctionsTable &external_functions_table_arg,
                    const StaticModel &static_model_arg);
+
+  SteadyStateModel(const SteadyStateModel &m);
+  SteadyStateModel(SteadyStateModel &&) = delete;
+  SteadyStateModel & operator=(const SteadyStateModel &m);
+  SteadyStateModel & operator=(SteadyStateModel &&) = delete;
+
   //! Add an expression of the form "var = expr;"
   void addDefinition(int symb_id, expr_t expr);
   //! Add an expression of the form "[ var1, var2, ... ] = expr;"
@@ -76,6 +82,11 @@ public:
            TrendComponentModelTable &trend_component_model_table_arg,
            VarModelTable &var_model_table_arg);
 
+  Epilogue(const Epilogue &m);
+  Epilogue(Epilogue &&) = delete;
+  Epilogue & operator=(const Epilogue &m);
+  Epilogue & operator=(Epilogue &&) = delete;
+
   //! Add an expression of the form "var = expr;"
   void addDefinition(int symb_id, expr_t expr);
 
diff --git a/src/ModelTree.cc b/src/ModelTree.cc
index d2d35b65..bb129b42 100644
--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
@@ -32,6 +32,148 @@
 
 using namespace MFS;
 
+void
+ModelTree::copyHelper(const ModelTree &m)
+{
+  auto f = [this](expr_t e) { return e->cloneDynamic(*this); };
+
+  // Equations
+  for (const auto & it : m.equations)
+    equations.push_back(dynamic_cast<BinaryOpNode *>(f(it)));
+  for (const auto & it : m.aux_equations)
+    aux_equations.push_back(dynamic_cast<BinaryOpNode *>(f(it)));
+
+  // Derivatives
+  for (const auto & it : m.first_derivatives)
+    first_derivatives[it.first] = f(it.second);
+  for (const auto & it : m.second_derivatives)
+    second_derivatives[it.first] = f(it.second);
+  for (const auto & it : m.third_derivatives)
+    third_derivatives[it.first] = f(it.second);
+  for (const auto & it : m.residuals_params_derivatives)
+    residuals_params_derivatives[it.first] = f(it.second);
+  for (const auto & it : m.residuals_params_second_derivatives)
+    residuals_params_second_derivatives[it.first] = f(it.second);
+  for (const auto & it : m.jacobian_params_derivatives)
+    jacobian_params_derivatives[it.first] = f(it.second);
+  for (const auto & it : m.jacobian_params_second_derivatives)
+    jacobian_params_second_derivatives[it.first] = f(it.second);
+  for (const auto & it : m.hessian_params_derivatives)
+    hessian_params_derivatives[it.first] = f(it.second);
+
+  // Temporary terms
+  for (const auto & it : m.temporary_terms)
+    temporary_terms.insert(f(it));
+  for (const auto & it : m.temporary_terms_mlv)
+    temporary_terms_mlv[f(it.first)] = f(it.second);
+  for (const auto & it : m.temporary_terms_res)
+    temporary_terms_res.insert(f(it));
+  for (const auto & it : m.temporary_terms_g1)
+    temporary_terms_g1.insert(f(it));
+  for (const auto & it : m.temporary_terms_g2)
+    temporary_terms_g2.insert(f(it));
+  for (const auto & it : m.temporary_terms_g3)
+    temporary_terms_g3.insert(f(it));
+  for (const auto & it : m.temporary_terms_idxs)
+    temporary_terms_idxs[f(it.first)] = it.second;
+  for (const auto & it : m.params_derivs_temporary_terms)
+    params_derivs_temporary_terms.insert(f(it));
+  for (const auto & it : m.params_derivs_temporary_terms_res)
+    params_derivs_temporary_terms_res.insert(f(it));
+  for (const auto & it : m.params_derivs_temporary_terms_g1)
+    params_derivs_temporary_terms_g1.insert(f(it));
+  for (const auto & it : m.params_derivs_temporary_terms_res2)
+    params_derivs_temporary_terms_res2.insert(f(it));
+  for (const auto & it : m.params_derivs_temporary_terms_g12)
+    params_derivs_temporary_terms_g12.insert(f(it));
+  for (const auto & it : m.params_derivs_temporary_terms_g2)
+    params_derivs_temporary_terms_g2.insert(f(it));
+  for (const auto & it : m.params_derivs_temporary_terms_idxs)
+    params_derivs_temporary_terms_idxs[f(it.first)] = it.second;
+
+  // Other stuff
+  for (const auto & it : m.trend_symbols_map)
+    trend_symbols_map[it.first] = f(it.second);
+  for (const auto & it : m.nonstationary_symbols_map)
+    nonstationary_symbols_map[it.first] = make_pair(it.second.first, f(it.second.second));
+}
+
+ModelTree::ModelTree(const ModelTree &m) :
+  DataTree {m},
+  equations_lineno {m.equations_lineno},
+  equation_tags {m.equation_tags},
+  NNZDerivatives {m.NNZDerivatives},
+  equation_reordered {m.equation_reordered},
+  variable_reordered {m.variable_reordered},
+  inv_equation_reordered {m.inv_equation_reordered},
+  inv_variable_reordered {m.inv_variable_reordered},
+  is_equation_linear {m.is_equation_linear},
+  endo2eq {m.endo2eq},
+  epilogue {m.epilogue},
+  prologue {m.prologue},
+  block_lag_lead {m.block_lag_lead},
+  cutoff {m.cutoff},
+  mfs {m.mfs}
+{
+  copyHelper(m);
+}
+
+ModelTree &
+ModelTree::operator=(const ModelTree &m)
+{
+  DataTree::operator=(m);
+
+  equations.clear();
+  equations_lineno = m.equations_lineno;
+  aux_equations.clear();
+  equation_tags = m.equation_tags;
+  NNZDerivatives = m.NNZDerivatives;
+
+  first_derivatives.clear();
+  second_derivatives.clear();
+  third_derivatives.clear();
+  residuals_params_derivatives.clear();
+  residuals_params_second_derivatives.clear();
+  jacobian_params_derivatives.clear();
+  jacobian_params_second_derivatives.clear();
+  hessian_params_derivatives.clear();
+
+  temporary_terms.clear();
+  temporary_terms_mlv.clear();
+  temporary_terms_res.clear();
+  temporary_terms_g1.clear();
+  temporary_terms_g2.clear();
+  temporary_terms_g3.clear();
+  temporary_terms_idxs.clear();
+  params_derivs_temporary_terms.clear();
+  params_derivs_temporary_terms_res.clear();
+  params_derivs_temporary_terms_g1.clear();
+  params_derivs_temporary_terms_res2.clear();
+  params_derivs_temporary_terms_g12.clear();
+  params_derivs_temporary_terms_g2.clear();
+  params_derivs_temporary_terms_idxs.clear();
+
+  trend_symbols_map.clear();
+  nonstationary_symbols_map.clear();
+
+  equation_reordered = m.equation_reordered;
+  variable_reordered = m.variable_reordered;
+  inv_equation_reordered = m.inv_equation_reordered;
+  inv_variable_reordered = m.inv_variable_reordered;
+  is_equation_linear = m.is_equation_linear;
+  endo2eq = m.endo2eq;
+  epilogue = m.epilogue;
+  prologue = m.prologue;
+  block_lag_lead = m.block_lag_lead;
+  cutoff = m.cutoff;
+  mfs = m.mfs;
+
+  copyHelper(m);
+
+  return *this;
+}
+
+
 bool
 ModelTree::computeNormalization(const jacob_map_t &contemporaneous_jacobian, bool verbose)
 {
@@ -1133,8 +1275,9 @@ ModelTree::BlockLinear(const blocks_derivatives_t &blocks_derivatives, const vec
 
 ModelTree::ModelTree(SymbolTable &symbol_table_arg,
                      NumericalConstants &num_constants_arg,
-                     ExternalFunctionsTable &external_functions_table_arg) :
-  DataTree{symbol_table_arg, num_constants_arg, external_functions_table_arg}
+                     ExternalFunctionsTable &external_functions_table_arg,
+                     bool is_dynamic_arg) :
+  DataTree {symbol_table_arg, num_constants_arg, external_functions_table_arg, is_dynamic_arg}
 {
   for (int & NNZDerivative : NNZDerivatives)
     NNZDerivative = 0;
diff --git a/src/ModelTree.hh b/src/ModelTree.hh
index bdbe1cc4..8f746a77 100644
--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
@@ -74,7 +74,7 @@ protected:
   vector<pair<int, pair<string, string>>> equation_tags;
 
   //! Number of non-zero derivatives
-  int NNZDerivatives[3];
+  array<int, 3> NNZDerivatives;
 
   using first_derivatives_t = map<pair<int, int>, expr_t>;
   //! First order derivatives
@@ -331,10 +331,24 @@ protected:
   virtual int getBlockInitialOtherEndogenousID(int block_number, int variable_number) const = 0;
   //! Initialize equation_reordered & variable_reordered
   void initializeVariablesAndEquations();
+
+private:
+  //! Internal helper for the copy constructor and assignment operator
+  /*! Copies all the structures that contain ExprNode*, by the converting the
+      pointers into their equivalent in the new tree */
+  void copyHelper(const ModelTree &m);
+
 public:
   ModelTree(SymbolTable &symbol_table_arg,
             NumericalConstants &num_constants_arg,
-            ExternalFunctionsTable &external_functions_table_arg);
+            ExternalFunctionsTable &external_functions_table_arg,
+            bool is_dynamic_arg = false);
+
+  ModelTree(const ModelTree &m);
+  ModelTree(ModelTree &&) = delete;
+  ModelTree & operator=(const ModelTree &m);
+  ModelTree & operator=(ModelTree &&) = delete;
+
   //! Absolute value under which a number is considered to be zero
   double cutoff{1e-15};
   //! Compute the minimum feedback set
diff --git a/src/ParsingDriver.cc b/src/ParsingDriver.cc
index 0cbf7e4a..4af2922c 100644
--- a/src/ParsingDriver.cc
+++ b/src/ParsingDriver.cc
@@ -2132,10 +2132,10 @@ ParsingDriver::run_model_comparison()
 void
 ParsingDriver::begin_planner_objective()
 {
-  planner_objective_statement = make_unique<PlannerObjectiveStatement>(mod_file->symbol_table,
-                                                                       mod_file->num_constants,
-                                                                       mod_file->external_functions_table);
-  set_current_data_tree(&planner_objective_statement->getPlannerObjective());
+  planner_objective = make_unique<StaticModel>(mod_file->symbol_table,
+                                               mod_file->num_constants,
+                                               mod_file->external_functions_table);
+  set_current_data_tree(planner_objective.get());
 }
 
 void
@@ -2145,7 +2145,7 @@ ParsingDriver::end_planner_objective(expr_t expr)
   expr_t eq = model_tree->AddEqual(expr, model_tree->Zero);
   model_tree->addEquation(eq, location.begin.line);
 
-  mod_file->addStatement(move(planner_objective_statement));
+  mod_file->addStatement(make_unique<PlannerObjectiveStatement>(*planner_objective));
 
   reset_data_tree();
 }
diff --git a/src/ParsingDriver.hh b/src/ParsingDriver.hh
index 597f41e4..9eef5f57 100644
--- a/src/ParsingDriver.hh
+++ b/src/ParsingDriver.hh
@@ -119,7 +119,7 @@ private:
   SymbolList symbol_list;
 
   //! Temporary store for the planner objective
-  unique_ptr<PlannerObjectiveStatement> planner_objective_statement;
+  unique_ptr<StaticModel> planner_objective;
 
   //! The data tree in which to add expressions currently parsed
   /*! The object pointed to is not owned by the parsing driver. It is essentially a
diff --git a/src/StaticModel.cc b/src/StaticModel.cc
index b164dc14..1c97aebf 100644
--- a/src/StaticModel.cc
+++ b/src/StaticModel.cc
@@ -29,6 +29,63 @@
 
 #include "StaticModel.hh"
 
+void
+StaticModel::copyHelper(const StaticModel &m)
+{
+  auto f = [this](expr_t e) { return e->cloneDynamic(*this); };
+
+  auto convert_vector_tt = [this,f](vector<temporary_terms_t> vtt)
+    {
+      vector<temporary_terms_t> vtt2;
+      for (const auto &tt : vtt)
+        {
+          temporary_terms_t tt2;
+          for (const auto &it : tt)
+            tt2.insert(f(it));
+          vtt2.push_back(tt2);
+        }
+      return vtt2;
+    };
+
+  for (const auto &it : m.v_temporary_terms)
+    v_temporary_terms.push_back(convert_vector_tt(it));
+  for (const auto &it : m.v_temporary_terms_local)
+    v_temporary_terms_local.push_back(convert_vector_tt(it));
+
+  for (const auto &it : m.first_chain_rule_derivatives)
+    first_chain_rule_derivatives[it.first] = f(it.second);
+
+  for (const auto &it : m.equation_type_and_normalized_equation)
+    equation_type_and_normalized_equation.push_back(make_pair(it.first, f(it.second)));
+
+  for (const auto &it : m.blocks_derivatives)
+    {
+      block_derivatives_equation_variable_laglead_nodeid_t v;
+      for (const auto &it2 : it)
+        v.push_back(make_pair(it2.first, make_pair(it2.second.first, f(it2.second.second))));
+      blocks_derivatives.push_back(v);
+    }
+
+  for (const auto &it : m.dynamic_jacobian)
+    dynamic_jacobian[it.first] = f(it.second);
+
+  auto convert_derivative_t = [this,f](derivative_t dt)
+    {
+      derivative_t dt2;
+      for (const auto &it : dt)
+        dt2[it.first] = f(it.second);
+      return dt2;
+    };
+  for (const auto &it : m.derivative_endo)
+    derivative_endo.push_back(convert_derivative_t(it));
+  for (const auto &it : m.derivative_other_endo)
+    derivative_other_endo.push_back(convert_derivative_t(it));
+  for (const auto &it : m.derivative_exo)
+    derivative_exo.push_back(convert_derivative_t(it));
+  for (const auto &it : m.derivative_exo_det)
+    derivative_exo_det.push_back(convert_derivative_t(it));
+}
+
 StaticModel::StaticModel(SymbolTable &symbol_table_arg,
                          NumericalConstants &num_constants_arg,
                          ExternalFunctionsTable &external_functions_table_arg) :
@@ -36,6 +93,76 @@ StaticModel::StaticModel(SymbolTable &symbol_table_arg,
 {
 }
 
+StaticModel::StaticModel(const StaticModel &m) :
+  ModelTree {m},
+  v_temporary_terms_inuse {m.v_temporary_terms_inuse},
+  map_idx {m.map_idx},
+  map_idx2 {m.map_idx2},
+  global_temporary_terms {m.global_temporary_terms},
+  block_type_firstequation_size_mfs {m.block_type_firstequation_size_mfs},
+  blocks_linear {m.blocks_linear},
+  other_endo_block {m.other_endo_block},
+  exo_block {m.exo_block},
+  exo_det_block {m.exo_det_block},
+  block_col_type {m.block_col_type},
+  variable_block_lead_lag {m.variable_block_lead_lag},
+  equation_block {m.equation_block},
+  endo_max_leadlag_block {m.endo_max_leadlag_block},
+  other_endo_max_leadlag_block {m.other_endo_max_leadlag_block},
+  exo_max_leadlag_block {m.exo_max_leadlag_block},
+  exo_det_max_leadlag_block {m.exo_det_max_leadlag_block},
+  max_leadlag_block {m.max_leadlag_block}
+{
+  copyHelper(m);
+}
+
+StaticModel &
+StaticModel::operator=(const StaticModel &m)
+{
+  ModelTree::operator=(m);
+
+  v_temporary_terms.clear();
+  v_temporary_terms_local.clear();
+
+  v_temporary_terms_inuse = m.v_temporary_terms_inuse;
+
+  first_chain_rule_derivatives.clear();
+
+  map_idx = m.map_idx;
+  map_idx2 = m.map_idx2;
+  global_temporary_terms = m.global_temporary_terms;
+
+  equation_type_and_normalized_equation.clear();
+
+  block_type_firstequation_size_mfs = m.block_type_firstequation_size_mfs;
+
+  blocks_derivatives.clear();
+  dynamic_jacobian.clear();
+
+  blocks_linear = m.blocks_linear;
+
+  derivative_endo.clear();
+  derivative_other_endo.clear();
+  derivative_exo.clear();
+  derivative_exo_det.clear();
+
+  other_endo_block = m.other_endo_block;
+  exo_block = m.exo_block;
+  exo_det_block = m.exo_det_block;
+  block_col_type = m.block_col_type;
+  variable_block_lead_lag = m.variable_block_lead_lag;
+  equation_block = m.equation_block;
+  endo_max_leadlag_block = m.endo_max_leadlag_block;
+  other_endo_max_leadlag_block = m.other_endo_max_leadlag_block;
+  exo_max_leadlag_block = m.exo_max_leadlag_block;
+  exo_det_max_leadlag_block = m.exo_det_max_leadlag_block;
+  max_leadlag_block = m.max_leadlag_block;
+
+  copyHelper(m);
+
+  return *this;
+}
+
 void
 StaticModel::compileDerivative(ofstream &code_file, unsigned int &instruction_number, int eq, int symb_id, map_idx_t &map_idx, temporary_terms_t temporary_terms) const
 {
diff --git a/src/StaticModel.hh b/src/StaticModel.hh
index ca6b0f63..300d3dad 100644
--- a/src/StaticModel.hh
+++ b/src/StaticModel.hh
@@ -104,7 +104,6 @@ private:
   //! Collecte the derivatives w.r. to endogenous of the block, to endogenous of previouys blocks and to exogenous
   void collect_block_first_order_derivatives();
 
-protected:
   //! Indicate if the temporary terms are computed for the overall model (true) or not (false). Default value true
   bool global_temporary_terms{true};
 
@@ -159,11 +158,22 @@ protected:
 
   void writeStaticModel(ostream &DynamicOutput, bool use_dll, bool julia) const;
   void writeStaticModel(const string &dynamic_basename, bool use_dll, bool julia) const;
+
+  //! Internal helper for the copy constructor and assignment operator
+  /*! Copies all the structures that contain ExprNode*, by the converting the
+      pointers into their equivalent in the new tree */
+  void copyHelper(const StaticModel &m);
+
 public:
   StaticModel(SymbolTable &symbol_table_arg,
               NumericalConstants &num_constants,
               ExternalFunctionsTable &external_functions_table_arg);
 
+  StaticModel(const StaticModel &m);
+  StaticModel(StaticModel &&) = delete;
+  StaticModel & operator=(const StaticModel &m);
+  StaticModel & operator=(StaticModel &&) = delete;
+
   //! Writes information on block decomposition when relevant
   void writeOutput(ostream &output, bool block) const;