diff --git a/DynamicModel.cc b/DynamicModel.cc
index ee8ffad8b0e1df00529d81941179bfaee515a4d0..da88f079f2d03d59272665efe24cad40d91a2a65 100644
--- a/DynamicModel.cc
+++ b/DynamicModel.cc
@@ -1558,14 +1558,13 @@ DynamicModel::writeDynamicMFile(const string &dynamic_basename) const
                     << "% Warning : this file is generated automatically by Dynare" << endl
                     << "%           from model file (.mod)" << endl << endl;
 
-  writeVarExpectationCalls(mDynamicModelFile);
   writeDynamicModel(mDynamicModelFile, false, false);
   mDynamicModelFile << "end" << endl; // Close *_dynamic function
   mDynamicModelFile.close();
 }
 
 void
-DynamicModel::writeVarExpectationCalls(ofstream &output) const
+DynamicModel::writeVarExpectationCalls(ostream &output) const
 {
   map<string, int> alreadyWritten;
   for (size_t i = 0; i < equations.size(); i++)
@@ -2369,8 +2368,9 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll, bool julia
       DynamicOutput << "%" << endl
                     << "% Model equations" << endl
                     << "%" << endl
-                    << endl
-                    << "residual = zeros(" << nrows << ", 1);" << endl
+                    << endl;
+      writeVarExpectationCalls(DynamicOutput);
+      DynamicOutput << "residual = zeros(" << nrows << ", 1);" << endl
                     << model_local_vars_output.str()
                     << model_output.str()
         // Writing initialization instruction for matrix g1
diff --git a/DynamicModel.hh b/DynamicModel.hh
index 5655c454c46f1e02c0da0deca8ffca54617ffbfe..64c368cef9acee5f440b309ed35c1d406f546d09 100644
--- a/DynamicModel.hh
+++ b/DynamicModel.hh
@@ -86,7 +86,7 @@ private:
   //! Writes dynamic model file (Julia version)
   void writeDynamicJuliaFile(const string &dynamic_basename) const;
   //! Write Var Expectation calls
-  void writeVarExpectationCalls(ofstream &output) const;
+  void writeVarExpectationCalls(ostream &output) const;
   //! Writes dynamic model file (C version)
   /*! \todo add third derivatives handling */
   void writeDynamicCFile(const string &dynamic_basename, const int order) const;
diff --git a/ExprNode.cc b/ExprNode.cc
index 4d7296824c2fd23d70f265bef1b2fe05746140e2..ef1df3e3f3a190c93b026088097dcf092ef66c62 100644
--- a/ExprNode.cc
+++ b/ExprNode.cc
@@ -514,7 +514,7 @@ NumConstNode::setVarExpectationIndex(map<string, SymbolList> var_model_info)
 }
 
 void
-NumConstNode::writeVarExpectationCalls(ofstream &output, map<string, int> alreadyWritten) const
+NumConstNode::writeVarExpectationCalls(ostream &output, map<string, int> &alreadyWritten) const
 {
 }
 
@@ -1518,7 +1518,7 @@ VariableNode::setVarExpectationIndex(map<string, SymbolList> var_model_info)
 }
 
 void
-VariableNode::writeVarExpectationCalls(ofstream &output, map<string, int> alreadyWritten) const
+VariableNode::writeVarExpectationCalls(ostream &output, map<string, int> &alreadyWritten) const
 {
 }
 
@@ -2601,7 +2601,7 @@ UnaryOpNode::setVarExpectationIndex(map<string, SymbolList> var_model_info)
 }
 
 void
-UnaryOpNode::writeVarExpectationCalls(ofstream &output, map<string, int> alreadyWritten) const
+UnaryOpNode::writeVarExpectationCalls(ostream &output, map<string, int> &alreadyWritten) const
 {
   arg->writeVarExpectationCalls(output, alreadyWritten);
 }
@@ -3917,7 +3917,7 @@ BinaryOpNode::setVarExpectationIndex(map<string, SymbolList> var_model_info)
 }
 
 void
-BinaryOpNode::writeVarExpectationCalls(ofstream &output, map<string, int> alreadyWritten) const
+BinaryOpNode::writeVarExpectationCalls(ostream &output, map<string, int> &alreadyWritten) const
 {
   arg1->writeVarExpectationCalls(output, alreadyWritten);
   arg2->writeVarExpectationCalls(output, alreadyWritten);
@@ -4604,7 +4604,7 @@ TrinaryOpNode::setVarExpectationIndex(map<string, SymbolList> var_model_info)
 }
 
 void
-TrinaryOpNode::writeVarExpectationCalls(ofstream &output, map<string, int> alreadyWritten) const
+TrinaryOpNode::writeVarExpectationCalls(ostream &output, map<string, int> &alreadyWritten) const
 {
   arg1->writeVarExpectationCalls(output, alreadyWritten);
   arg2->writeVarExpectationCalls(output, alreadyWritten);
@@ -4926,7 +4926,7 @@ AbstractExternalFunctionNode::setVarExpectationIndex(map<string, SymbolList> var
 }
 
 void
-AbstractExternalFunctionNode::writeVarExpectationCalls(ofstream &output, map<string, int> alreadyWritten) const
+AbstractExternalFunctionNode::writeVarExpectationCalls(ostream &output, map<string, int> &alreadyWritten) const
 {
   for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
     (*it)->writeVarExpectationCalls(output, alreadyWritten);
@@ -5929,7 +5929,7 @@ VarExpectationNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
 }
 
 void
-VarExpectationNode::writeVarExpectationCalls(ofstream &output, map<string, int> alreadyWritten) const
+VarExpectationNode::writeVarExpectationCalls(ostream &output, map<string, int> &alreadyWritten) const
 {
   map<string, int>::iterator it = alreadyWritten.find(model_name);
   if (it != alreadyWritten.end() && alreadyWritten[model_name] == forecast_horizon)
diff --git a/ExprNode.hh b/ExprNode.hh
index 5021c96cb9ae66a325d9542a3346df199af183b3..7ba9031062ff18cc9a76d06c86e4512e38cd7b50 100644
--- a/ExprNode.hh
+++ b/ExprNode.hh
@@ -456,7 +456,7 @@ public:
   virtual void setVarExpectationIndex(map<string, SymbolList> var_model_info) = 0;
 
   // Write calls to var forecast and place in temporary variable
-  virtual void writeVarExpectationCalls(ofstream &output, map<string, int>) const = 0;
+  virtual void writeVarExpectationCalls(ostream &output, map<string, int> &alreadyWritten) const = 0;
 };
 
 //! Object used to compare two nodes (using their indexes)
@@ -518,7 +518,7 @@ public:
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
   virtual bool isInStaticForm() const;
   virtual void setVarExpectationIndex(map<string, SymbolList> var_model_info);
-  virtual void writeVarExpectationCalls(ofstream &output, map<string, int>) const;
+  virtual void writeVarExpectationCalls(ostream &output, map<string, int> &alreadyWritten) const;
   virtual expr_t substituteStaticAuxiliaryVariable() const;
 };
 
@@ -586,7 +586,7 @@ public:
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
   virtual bool isInStaticForm() const;
   virtual void setVarExpectationIndex(map<string, SymbolList> var_model_info);
-  virtual void writeVarExpectationCalls(ofstream &output, map<string, int>) const;
+  virtual void writeVarExpectationCalls(ostream &output, map<string, int> &alreadyWritten) const;
   //! Substitute auxiliary variables by their expression in static model
   virtual expr_t substituteStaticAuxiliaryVariable() const;
 };
@@ -674,7 +674,7 @@ public:
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
   virtual bool isInStaticForm() const;
   virtual void setVarExpectationIndex(map<string, SymbolList> var_model_info);
-  virtual void writeVarExpectationCalls(ofstream &output, map<string, int>) const;
+  virtual void writeVarExpectationCalls(ostream &output, map<string, int> &alreadyWritten) const;
   //! Substitute auxiliary variables by their expression in static model
   virtual expr_t substituteStaticAuxiliaryVariable() const;
 };
@@ -781,7 +781,7 @@ public:
   expr_t getNonZeroPartofEquation() const;
   virtual bool isInStaticForm() const;
   virtual void setVarExpectationIndex(map<string, SymbolList> var_model_info);
-  virtual void writeVarExpectationCalls(ofstream &output, map<string, int>) const;
+  virtual void writeVarExpectationCalls(ostream &output, map<string, int> &alreadyWritten) const;
   //! Substitute auxiliary variables by their expression in static model
   virtual expr_t substituteStaticAuxiliaryVariable() const;
   //! Substitute auxiliary variables by their expression in static model auxiliary variable definition
@@ -858,7 +858,7 @@ public:
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
   virtual bool isInStaticForm() const;
   virtual void setVarExpectationIndex(map<string, SymbolList> var_model_info);
-  virtual void writeVarExpectationCalls(ofstream &output, map<string, int>) const;
+  virtual void writeVarExpectationCalls(ostream &output, map<string, int> &alreadyWritten) const;
   //! Substitute auxiliary variables by their expression in static model
   virtual expr_t substituteStaticAuxiliaryVariable() const;
 };
@@ -942,7 +942,7 @@ public:
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
   virtual bool isInStaticForm() const;
   virtual void setVarExpectationIndex(map<string, SymbolList> var_model_info);
-  virtual void writeVarExpectationCalls(ofstream &output, map<string, int>) const;
+  virtual void writeVarExpectationCalls(ostream &output, map<string, int> &alreadyWritten) const;
   //! Substitute auxiliary variables by their expression in static model
   virtual expr_t substituteStaticAuxiliaryVariable() const;
 };
@@ -1110,7 +1110,7 @@ public:
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
   virtual bool isInStaticForm() const;
   virtual void setVarExpectationIndex(map<string, SymbolList> var_model_info);
-  virtual void writeVarExpectationCalls(ofstream &output, map<string, int>) const;
+  virtual void writeVarExpectationCalls(ostream &output, map<string, int> &alreadyWritten) const;
   virtual expr_t substituteStaticAuxiliaryVariable() const;
 };