diff --git a/DynamicModel.cc b/DynamicModel.cc
index 10db3a8b9bece10aef9bf2192ac20feb7277b0b4..1bd65d7299c6002ea1574ca5f9cc2444d39f300a 100644
--- a/DynamicModel.cc
+++ b/DynamicModel.cc
@@ -3164,7 +3164,7 @@ DynamicModel::runTrendTest(const eval_context_t &eval_context)
 void
 DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivatives, int paramsDerivsOrder,
                             const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll,
-                            bool bytecode, bool compute_xrefs)
+                            bool bytecode)
 {
   assert(jacobianExo || !(hessian || thirdDerivatives || paramsDerivsOrder));
 
@@ -3272,9 +3272,113 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
         if (bytecode)
           computeTemporaryTermsMapping();
       }
+}
 
-  if (compute_xrefs)
-    computeXrefs();
+void
+DynamicModel::computeXrefs()
+{
+  int i = 0;
+  for (vector<BinaryOpNode *>::iterator it = equations.begin();
+       it != equations.end(); it++)
+    {
+      ExprNode::EquationInfo ei;
+      (*it)->computeXrefs(ei);
+      xrefs[i++] = ei;
+    }
+
+  i = 0;
+  for (map<int, ExprNode::EquationInfo>::const_iterator it = xrefs.begin();
+       it != xrefs.end(); it++, i++)
+    {
+      computeRevXref(xref_param, it->second.param, i);
+      computeRevXref(xref_endo, it->second.endo, i);
+      computeRevXref(xref_exo, it->second.exo, i);
+      computeRevXref(xref_exo_det, it->second.exo_det, i);
+    }
+}
+
+void
+DynamicModel::computeRevXref(map<pair<int, int>, set<int> > &xrefset, const set<pair<int, int> > &eiref, int eqn)
+{
+  for (set<pair<int, int> >::const_iterator it = eiref.begin();
+       it != eiref.end(); it++)
+    {
+      set<int> eq;
+      if (xrefset.find(*it) != xrefset.end())
+        eq = xrefset[*it];
+      eq.insert(eqn);
+      xrefset[*it] = eq;
+    }
+}
+
+void
+DynamicModel::writeXrefs(ostream &output) const
+{
+  output << "M_.xref1.param = cell(1, M_.eq_nbr);" << endl
+         << "M_.xref1.endo = cell(1, M_.eq_nbr);" << endl
+         << "M_.xref1.exo = cell(1, M_.eq_nbr);" << endl
+         << "M_.xref1.exo_det = cell(1, M_.eq_nbr);" << endl;
+  int i = 1;
+  for (map<int, ExprNode::EquationInfo>::const_iterator it = xrefs.begin();
+       it != xrefs.end(); it++, i++)
+    {
+      output << "M_.xref1.param{" << i << "} = [ ";
+      for (set<pair<int, int> >::const_iterator it1 = it->second.param.begin();
+           it1 != it->second.param.end(); it1++)
+        output << symbol_table.getTypeSpecificID(it1->first) + 1 << " ";
+      output << "];" << endl;
+
+      output << "M_.xref1.endo{" << i << "} = [ ";
+      for (set<pair<int, int> >::const_iterator it1 = it->second.endo.begin();
+           it1 != it->second.endo.end(); it1++)
+        output << "struct('id', " << symbol_table.getTypeSpecificID(it1->first) + 1 << ", 'shift', " << it1->second << ");";
+      output << "];" << endl;
+
+      output << "M_.xref1.exo{" << i << "} = [ ";
+      for (set<pair<int, int> >::const_iterator it1 = it->second.exo.begin();
+           it1 != it->second.exo.end(); it1++)
+        output << "struct('id', " << symbol_table.getTypeSpecificID(it1->first) + 1 << ", 'shift', " << it1->second << ");";
+      output << "];" << endl;
+
+      output << "M_.xref1.exo_det{" << i << "} = [ ";
+      for (set<pair<int, int> >::const_iterator it1 = it->second.exo_det.begin();
+           it1 != it->second.exo_det.end(); it1++)
+        output << "struct('id', " << symbol_table.getTypeSpecificID(it1->first) + 1 << ", 'shift', " << it1->second << ");";
+      output << "];" << endl;
+    }
+
+  output << "M_.xref2.param = cell(1, M_.param_nbr);" << endl
+         << "M_.xref2.endo = cell(1, M_.endo_nbr);" << endl
+         << "M_.xref2.exo = cell(1, M_.exo_nbr);" << endl
+         << "M_.xref2.exo_det = cell(1, M_.exo_det_nbr);" << endl;
+  writeRevXrefs(output, xref_param, "param");
+  writeRevXrefs(output, xref_endo, "endo");
+  writeRevXrefs(output, xref_exo, "exo");
+  writeRevXrefs(output, xref_exo_det, "exo_det");
+}
+
+void
+DynamicModel::writeRevXrefs(ostream &output, const map<pair<int, int>, set<int> > &xrefmap, const string &type) const
+{
+  int last_tsid = -1;
+  for (map<pair<int, int>, set<int> >::const_iterator it = xrefmap.begin();
+       it != xrefmap.end(); it++)
+    {
+      int tsid = symbol_table.getTypeSpecificID(it->first.first) + 1;
+      output << "M_.xref2." << type << "{" << tsid << "} = [ ";
+      if (last_tsid == tsid)
+        output << "M_.xref2." << type << "{" << tsid << "}; ";
+      else
+        last_tsid = tsid;
+
+      for (set<int>::const_iterator it1 = it->second.begin();
+           it1 != it->second.end(); it1++)
+        if (type == "param")
+          output << *it1 + 1 << " ";
+        else
+          output << "struct('shift', " << it->first.second << ", 'eq', " << *it1+1 << ");";
+      output << "];" << endl;
+    }
 }
 
 map<pair<pair<int, pair<int, int> >, pair<int, int> >, int>
diff --git a/DynamicModel.hh b/DynamicModel.hh
index b0828a2b23e6de49aefcd80c2850d5b112b55de3..6001c017a4ba79c269c28551227f35eeb95820a2 100644
--- a/DynamicModel.hh
+++ b/DynamicModel.hh
@@ -62,6 +62,13 @@ private:
   /*! Set by computeDerivIDs() */
   int max_exo_det_lag, max_exo_det_lead;
 
+  //! Cross reference information
+  map<int, ExprNode::EquationInfo> xrefs;
+  map<pair<int, int>, set<int> > xref_param;
+  map<pair<int, int>, set<int> > xref_endo;
+  map<pair<int, int>, set<int> > xref_exo;
+  map<pair<int, int>, set<int> > xref_exo_det;
+
   //! Number of columns of dynamic jacobian
   /*! Set by computeDerivID()s and computeDynJacobianCols() */
   int dynJacobianColsNbr;
@@ -189,6 +196,12 @@ private:
   /*! pair< pair<static, forward>, pair<backward,mixed> > */
   vector<pair< pair<int, int>, pair<int, int> > > block_col_type;
 
+  //! Help computeXrefs to compute the reverse references (i.e. param->eqs, endo->eqs, etc)
+  void computeRevXref(map<pair<int, int>, set<int> > &xrefset, const set<pair<int, int> > &eiref, int eqn);
+
+  //! Write reverse cross references
+  void writeRevXrefs(ostream &output, const map<pair<int, int>, set<int> > &xrefmap, const string &type) const;
+
   //! List for each variable its block number and its maximum lag and lead inside the block
   vector<pair<int, pair<int, int> > > variable_block_lead_lag;
   //! List for each equation its block number
@@ -202,7 +215,13 @@ public:
   //! Adds a variable node
   /*! This implementation allows for non-zero lag */
   virtual VariableNode *AddVariable(int symb_id, int lag = 0);
-  
+
+  //! Compute cross references
+  void computeXrefs();
+
+  //! Write cross references
+  void writeXrefs(ostream &output) const;
+
   //! Execute computations (variable sorting + derivation)
   /*!
     \param jacobianExo whether derivatives w.r. to exo and exo_det should be in the Jacobian (derivatives w.r. to endo are always computed)
@@ -213,7 +232,7 @@ public:
     \param no_tmp_terms if true, no temporary terms will be computed in the dynamic files
   */
   void computingPass(bool jacobianExo, bool hessian, bool thirdDerivatives, int paramsDerivsOrder,
-                     const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll, bool bytecode, bool compute_xrefs);
+                     const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll, bool bytecode);
   //! Writes model initialization and lead/lag incidence matrix to output
   void writeOutput(ostream &output, const string &basename, bool block, bool byte_code, bool use_dll, int order, bool estimation_present, bool compute_xrefs, bool julia) const;
 
diff --git a/DynareMain2.cc b/DynareMain2.cc
index c3a2493652691dd269e52b07f4e06ed7638fbaf2..5a870a62e11a25e4a3c17219f443e798f5b797ce 100644
--- a/DynareMain2.cc
+++ b/DynareMain2.cc
@@ -50,7 +50,7 @@ main2(stringstream &in, string &basename, bool debug, bool clear_all, bool clear
     mod_file->writeJsonOutput(basename, json, json_output_mode, onlyjson);
 
   // Perform transformations on the model (creation of auxiliary vars and equations)
-  mod_file->transformPass(nostrict);
+  mod_file->transformPass(nostrict, compute_xrefs);
   if (json == transformpass)
     mod_file->writeJsonOutput(basename, json, json_output_mode, onlyjson);
 
@@ -58,7 +58,7 @@ main2(stringstream &in, string &basename, bool debug, bool clear_all, bool clear
   mod_file->evalAllExpressions(warn_uninit);
 
   // Do computations
-  mod_file->computingPass(no_tmp_terms, output_mode, compute_xrefs, params_derivs_order);
+  mod_file->computingPass(no_tmp_terms, output_mode, params_derivs_order);
   if (json == computingpass)
     mod_file->writeJsonOutput(basename, json, json_output_mode, onlyjson, jsonprintderivdetail);
 
diff --git a/ExprNode.cc b/ExprNode.cc
index 5cc0298bdca44223c3b8b021eee6761b34f09ef7..1c03bc03ad3b79483338a169935bb1e3d10fed76 100644
--- a/ExprNode.cc
+++ b/ExprNode.cc
@@ -1109,16 +1109,16 @@ VariableNode::computeXrefs(EquationInfo &ei) const
   switch (type)
     {
     case eEndogenous:
-      ei.endo.insert(symb_id);
+      ei.endo.insert(make_pair(symb_id, lag));
       break;
     case eExogenous:
-      ei.exo.insert(symb_id);
+      ei.exo.insert(make_pair(symb_id, lag));
       break;
     case eExogenousDet:
-      ei.exo_det.insert(symb_id);
+      ei.exo_det.insert(make_pair(symb_id, lag));
       break;
     case eParameter:
-      ei.param.insert(symb_id);
+      ei.param.insert(make_pair(symb_id, 0));
       break;
     case eTrend:
     case eLogTrend:
diff --git a/ExprNode.hh b/ExprNode.hh
index e4e5d9bcd0fd7d9d4dc1fbfae070116ac15e1af5..7b5e246331ac4473fa591d8f9d5aca33daaf66ad 100644
--- a/ExprNode.hh
+++ b/ExprNode.hh
@@ -165,10 +165,10 @@ protected:
   //! For creating equation cross references
   struct EquationInfo
   {
-    set<int> param;
-    set<int> endo;
-    set<int> exo;
-    set<int> exo_det;
+    set<pair<int, int> > param;
+    set<pair<int, int> > endo;
+    set<pair<int, int> > exo;
+    set<pair<int, int> > exo_det;
   };
 
 public:
diff --git a/ModFile.cc b/ModFile.cc
index 7ee204ad42bed09eac3225e775a939c01242eca2..3c6ece94867dd7c35bc4f32fdad062beed6dfeb3 100644
--- a/ModFile.cc
+++ b/ModFile.cc
@@ -319,7 +319,7 @@ ModFile::checkPass(bool nostrict)
 }
 
 void
-ModFile::transformPass(bool nostrict)
+ModFile::transformPass(bool nostrict, bool compute_xrefs)
 {
   // Save the original model (must be done before any model transformations by preprocessor)
   dynamic_model.cloneDynamic(original_model);
@@ -414,6 +414,9 @@ ModFile::transformPass(bool nostrict)
   // Freeze the symbol table
   symbol_table.freeze();
 
+  if (compute_xrefs)
+    dynamic_model.computeXrefs();
+
   /*
     Enforce the same number of equations and endogenous, except in three cases:
     - ramsey_model, ramsey_policy or discretionary_policy is used
@@ -483,7 +486,7 @@ ModFile::transformPass(bool nostrict)
 }
 
 void
-ModFile::computingPass(bool no_tmp_terms, FileOutputType output, bool compute_xrefs, int params_derivs_order)
+ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_derivs_order)
 {
   // Mod file may have no equation (for example in a standalone BVAR estimation)
   if (dynamic_model.equation_number() > 0)
@@ -517,7 +520,7 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, bool compute_xr
 	  || mod_file_struct.calib_smoother_present)
 	{
 	  if (mod_file_struct.perfect_foresight_solver_present)
-	    dynamic_model.computingPass(true, false, false, none, global_eval_context, no_tmp_terms, block, use_dll, byte_code, compute_xrefs);
+	    dynamic_model.computingPass(true, false, false, none, global_eval_context, no_tmp_terms, block, use_dll, byte_code);
 	      else
 		{
 		  if (mod_file_struct.stoch_simul_present
@@ -542,13 +545,13 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, bool compute_xr
                   int paramsDerivsOrder = 0;
                   if (mod_file_struct.identification_present || mod_file_struct.estimation_analytic_derivation)
                     paramsDerivsOrder = params_derivs_order;
-		  dynamic_model.computingPass(true, hessian, thirdDerivatives, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code, compute_xrefs);
+		  dynamic_model.computingPass(true, hessian, thirdDerivatives, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code);
                   if (linear && mod_file_struct.ramsey_model_present)
-                    orig_ramsey_dynamic_model.computingPass(true, true, false, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code, compute_xrefs);
+                    orig_ramsey_dynamic_model.computingPass(true, true, false, paramsDerivsOrder, global_eval_context, no_tmp_terms, block, use_dll, byte_code);
 		}
 	    }
 	  else // No computing task requested, compute derivatives up to 2nd order by default
-	    dynamic_model.computingPass(true, true, false, none, global_eval_context, no_tmp_terms, block, use_dll, byte_code, compute_xrefs);
+	    dynamic_model.computingPass(true, true, false, none, global_eval_context, no_tmp_terms, block, use_dll, byte_code);
 
       if ((linear && !mod_file_struct.ramsey_model_present && !dynamic_model.checkHessianZero()) ||
           (linear && mod_file_struct.ramsey_model_present && !orig_ramsey_dynamic_model.checkHessianZero()))
diff --git a/ModFile.hh b/ModFile.hh
index c04f4de5e6cec6d5587d909dcf0bb2b683421831..935c5231283a21387eab8eb4903df9d7e67726a4 100644
--- a/ModFile.hh
+++ b/ModFile.hh
@@ -133,12 +133,12 @@ public:
   /*! \todo add check for number of equations and endogenous if ramsey_policy is present */
   void checkPass(bool nostrict);
   //! Perform some transformations on the model (creation of auxiliary vars and equations)
-  void transformPass(bool nostrict);
+  /*! \param compute_xrefs if true, equation cross references will be computed */
+  void transformPass(bool nostrict, bool compute_xrefs);
   //! Execute computations
   /*! \param no_tmp_terms if true, no temporary terms will be computed in the static and dynamic files */
-  /*! \param compute_xrefs if true, equation cross references will be computed */
   /*! \param params_derivs_order compute this order of derivs wrt parameters */
-  void computingPass(bool no_tmp_terms, FileOutputType output, bool compute_xrefs, int params_derivs_order);
+  void computingPass(bool no_tmp_terms, FileOutputType output, int params_derivs_order);
   //! Writes Matlab/Octave output files
   /*!
     \param basename The base name used for writing output files. Should be the name of the mod file without its extension
diff --git a/ModelTree.cc b/ModelTree.cc
index d1a270208f347d9365b02716aaecdfda7cbaada7..5b262870edcf0400ebeab69da04cc2cb7066ae67 100644
--- a/ModelTree.cc
+++ b/ModelTree.cc
@@ -229,103 +229,6 @@ ModelTree::computeNonSingularNormalization(jacob_map_t &contemporaneous_jacobian
     }
 }
 
-void
-ModelTree::computeXrefs()
-{
-  int i = 0;
-  for (vector<BinaryOpNode *>::iterator it = equations.begin();
-       it != equations.end(); it++)
-    {
-      ExprNode::EquationInfo ei;
-      (*it)->computeXrefs(ei);
-      xrefs[i++] = ei;
-    }
-
-  i = 0;
-  for (map<int, ExprNode::EquationInfo>::const_iterator it = xrefs.begin();
-       it != xrefs.end(); it++, i++)
-    {
-      computeRevXref(xref_param, it->second.param, i);
-      computeRevXref(xref_endo, it->second.endo, i);
-      computeRevXref(xref_exo, it->second.exo, i);
-      computeRevXref(xref_exo_det, it->second.exo_det, i);
-    }
-}
-
-void
-ModelTree::computeRevXref(map<int, set<int> > &xrefset, const set<int> &eiref, int eqn)
-{
-  for (set<int>::const_iterator it1 = eiref.begin();
-       it1 != eiref.end(); it1++)
-    {
-      set<int> eq;
-      if (xrefset.find(symbol_table.getTypeSpecificID(*it1)) != xrefset.end())
-        eq = xrefset[symbol_table.getTypeSpecificID(*it1)];
-      eq.insert(eqn);
-      xrefset[symbol_table.getTypeSpecificID(*it1)] = eq;
-    }
-}
-
-void
-ModelTree::writeXrefs(ostream &output) const
-{
-  output << "M_.xref1.param = cell(1, M_.eq_nbr);" << endl
-         << "M_.xref1.endo = cell(1, M_.eq_nbr);" << endl
-         << "M_.xref1.exo = cell(1, M_.eq_nbr);" << endl
-         << "M_.xref1.exo_det = cell(1, M_.eq_nbr);" << endl
-         << "M_.xref2.param = cell(1, M_.eq_nbr);" << endl
-         << "M_.xref2.endo = cell(1, M_.eq_nbr);" << endl
-         << "M_.xref2.exo = cell(1, M_.eq_nbr);" << endl
-         << "M_.xref2.exo_det = cell(1, M_.eq_nbr);" << endl;
-  int i = 1;
-  for (map<int, ExprNode::EquationInfo>::const_iterator it = xrefs.begin();
-       it != xrefs.end(); it++, i++)
-    {
-      output << "M_.xref1.param{" << i << "} = [ ";
-      for (set<int>::const_iterator it1 = it->second.param.begin();
-           it1 != it->second.param.end(); it1++)
-        output << symbol_table.getTypeSpecificID(*it1) + 1 << " ";
-      output << "];" << endl;
-
-      output << "M_.xref1.endo{" << i << "} = [ ";
-      for (set<int>::const_iterator it1 = it->second.endo.begin();
-           it1 != it->second.endo.end(); it1++)
-        output << symbol_table.getTypeSpecificID(*it1) + 1 << " ";
-      output << "];" << endl;
-
-      output << "M_.xref1.exo{" << i << "} = [ ";
-      for (set<int>::const_iterator it1 = it->second.exo.begin();
-           it1 != it->second.exo.end(); it1++)
-        output << symbol_table.getTypeSpecificID(*it1) + 1 << " ";
-      output << "];" << endl;
-
-      output << "M_.xref1.exo_det{" << i << "} = [ ";
-      for (set<int>::const_iterator it1 = it->second.exo_det.begin();
-           it1 != it->second.exo_det.end(); it1++)
-        output << symbol_table.getTypeSpecificID(*it1) + 1 << " ";
-      output << "];" << endl;
-    }
-
-  writeRevXrefs(output, xref_param, "param");
-  writeRevXrefs(output, xref_endo, "endo");
-  writeRevXrefs(output, xref_exo, "exo");
-  writeRevXrefs(output, xref_exo_det, "exo_det");
-}
-
-void
-ModelTree::writeRevXrefs(ostream &output, const map<int, set<int> > &xrefmap, const string &type) const
-{
-  for (map<int, set<int> >::const_iterator it = xrefmap.begin();
-       it != xrefmap.end(); it++)
-    {
-      output << "M_.xref2." << type << "{" << it->first + 1 << "} = [ ";
-      for (set<int>::const_iterator it1 = it->second.begin();
-           it1 != it->second.end(); it1++)
-        output << *it1 + 1 << " ";
-      output << "];" << endl;
-    }
-}
-
 void
 ModelTree::computeNormalizedEquations(multimap<int, int> &endo2eqs) const
 {
diff --git a/ModelTree.hh b/ModelTree.hh
index ac6bfca8e006a7687ecc5163b16c94174051106c..057bfc8e64d5b2453bb5dd098c8195c0013dea96 100644
--- a/ModelTree.hh
+++ b/ModelTree.hh
@@ -100,13 +100,6 @@ protected:
   */
   first_derivatives_t residuals_params_derivatives;
 
-  //! Cross reference information
-  map<int, ExprNode::EquationInfo> xrefs;
-  map<int, set<int> > xref_param;
-  map<int, set<int> > xref_endo;
-  map<int, set<int> > xref_exo;
-  map<int, set<int> > xref_exo_det;
-
   //! Second derivatives of the residuals w.r. to parameters
   /*! First index is equation number, second and third indeces are parameters.
     Only non-null derivatives are stored in the map.
@@ -246,14 +239,6 @@ protected:
 
   //! Try to normalized each unnormalized equation (matched endogenous variable only on the LHS)
   void computeNormalizedEquations(multimap<int, int> &endo2eqs) const;
-  //! Compute cross references
-  void computeXrefs();
-  //! Help computeXrefs to compute the reverse references (i.e. param->eqs, endo->eqs, etc)
-  void computeRevXref(map<int, set<int> > &xrefset, const set<int> &eiref, int eqn);
-  //! Write cross references
-  void writeXrefs(ostream &output) const;
-  //! Write reverse cross references
-  void writeRevXrefs(ostream &output, const map<int, set<int> > &xrefmap, const string &type) const;
   //! Evaluate the jacobian and suppress all the elements below the cutoff
   void evaluateAndReduceJacobian(const eval_context_t &eval_context, jacob_map_t &contemporaneous_jacobian, jacob_map_t &static_jacobian, dynamic_jacob_map_t &dynamic_jacobian, double cutoff, bool verbose);
   //! Search the equations and variables belonging to the prologue and the epilogue of the model