diff --git a/ModelTree.cc b/ModelTree.cc
index 3f2bcbcc6260f1c8995fee42619a0494ea5bad41..52830bba3893e3483987509c424e397e3f64ac04 100644
--- a/ModelTree.cc
+++ b/ModelTree.cc
@@ -240,6 +240,30 @@ ModelTree::computeXrefs()
       (*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
@@ -272,6 +296,25 @@ ModelTree::writeXrefs(ostream &output) const
         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
diff --git a/ModelTree.hh b/ModelTree.hh
index bb6ffb3fe3adf3cd1a0ed42e5b3745d62b0dab40..b68e3b3a08d907a67852e1ba020e87267a1b20ec 100644
--- a/ModelTree.hh
+++ b/ModelTree.hh
@@ -99,7 +99,12 @@ 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.
@@ -224,8 +229,12 @@ protected:
   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