diff --git a/matlab/dyn_ramsey_static.m b/matlab/dyn_ramsey_static.m
index c2b619b3f8b04735b1fb88d32d834acc816be270..dc19a6817cab3a78f80cdef0232c348820b9f587 100644
--- a/matlab/dyn_ramsey_static.m
+++ b/matlab/dyn_ramsey_static.m
@@ -143,8 +143,8 @@ end
 % the auxiliary variables before the Lagrange multipliers are treated
 % as ordinary endogenous variables
 aux_eq = [1:orig_endo_aux_nbr, orig_endo_aux_nbr+orig_eq_nbr+1:size(fJ,1)];
-A = fJ(aux_eq,orig_endo_aux_nbr+1:end);
-y = res(aux_eq);
+A = fJ(1:orig_endo_aux_nbr,orig_endo_nbr+find(aux_vars_type==6));
+y = res(1:orig_endo_aux_nbr);
 mult = -A\y;
 
 resids1 = y+A*mult;
diff --git a/matlab/evaluate_steady_state.m b/matlab/evaluate_steady_state.m
index aec3a1b3f3fc9cc46cf95377ce7fd3b662ab46ea..6f240898bab219e6f4433504a313cfc77b25ea16 100644
--- a/matlab/evaluate_steady_state.m
+++ b/matlab/evaluate_steady_state.m
@@ -125,7 +125,7 @@ function [ys,params,info] = evaluate_steady_state(ys_init,M,options,oo,steadysta
         %check whether steady state really solves the model
         resids = evaluate_static_model(ys,exo_ss,params,M,options);
 
-        n_multipliers=M.endo_nbr-M.orig_endo_nbr;
+        n_multipliers=M.orig_eq_nbr;
         nan_indices_multiplier=find(isnan(resids(1:n_multipliers)));
         nan_indices=find(isnan(resids(n_multipliers+1:end)));
 
@@ -196,9 +196,12 @@ function [ys,params,info] = evaluate_steady_state(ys_init,M,options,oo,steadysta
     elseif (options.bytecode == 0 && options.block == 0)
         if options.linear == 0
             % non linear model
-            [ys,check] = dynare_solve([M.fname '_static'],...
+            static_model = str2func([M.fname '_static']);
+            [ys,check] = dynare_solve(@static_problem,...
                                       ys_init,...
-                                      options, exo_ss, params);
+                                      options, exo_ss, params,...
+                                      M.orig_endo_nbr,...
+                                      static_model);
         else
             % linear model
             fh_static = str2func([M.fname '_static']);
@@ -294,3 +297,8 @@ function [ys,params,info] = evaluate_steady_state(ys_init,M,options,oo,steadysta
         info(2) = NaN;
         return
     end
+
+function [resids,jac] = static_problem(y,x,params,nvar,fh_static_model)
+    [r,j] = fh_static_model(y,x,params);
+    resids = r(1:nvar);
+    jac = j(1:nvar,1:nvar);
diff --git a/matlab/getH.m b/matlab/getH.m
index fdf84c4e9ea1d4b99609749c171651c7f6277cbc..902a93c5d192c82770ee78a5c91c4285b45284a9 100644
--- a/matlab/getH.m
+++ b/matlab/getH.m
@@ -130,26 +130,19 @@ else
 dyssdtheta=zeros(length(oo_.dr.ys),M_.param_nbr);
 d2yssdtheta=zeros(length(oo_.dr.ys),M_.param_nbr,M_.param_nbr);
 [residual, gg1] = feval([M_.fname,'_static'],oo_.dr.ys, oo_.exo_steady_state', M_.params);
-df = feval([M_.fname,'_params_derivs'],yy0, repmat(oo_.exo_steady_state',[M_.maximum_exo_lag+M_.maximum_exo_lead+1]), ...
-    M_.params, oo_.dr.ys, 1, dyssdtheta, d2yssdtheta);
+df = feval([M_.fname,'_static_params_derivs'],oo_.dr.ys, repmat(oo_.exo_steady_state',[M_.maximum_exo_lag+M_.maximum_exo_lead+1]), ...
+    M_.params);
 dyssdtheta = -gg1\df;
 if nargout>5,
     [residual, gg1, gg2] = feval([M_.fname,'_static'],oo_.dr.ys, oo_.exo_steady_state', M_.params);
     [residual, g1, g2, g3] = feval([M_.fname,'_dynamic'],yy0, oo_.exo_steady_state', ...
         M_.params, oo_.dr.ys, 1);
-    [nr, nc]=size(g2);
+    [nr, nc]=size(gg2);
 
-    [df, gp, d2f] = feval([M_.fname,'_params_derivs'],yy0, oo_.exo_steady_state', ...
-        M_.params, oo_.dr.ys, 1, dyssdtheta*0, d2yssdtheta);
+    [df, gpx, d2f] = feval([M_.fname,'_static_params_derivs'],oo_.dr.ys, oo_.exo_steady_state', ...
+        M_.params);%, oo_.dr.ys, 1, dyssdtheta*0, d2yssdtheta);
     d2f = get_all_resid_2nd_derivs(d2f,length(oo_.dr.ys),M_.param_nbr);
-    gpx = zeros(nr,nr,M_.param_nbr);
-    for j=1:nr,
-        for i=1:nr,
-            inx = I == i;
-            gpx(j,i,:)=sum(gp(j,inx,:),2);
-        end
-    end
-%     d2f = d2f(:,indx,indx);
+
     if isempty(find(gg2)),
         for j=1:M_.param_nbr,
         d2yssdtheta(:,:,j) = -gg1\d2f(:,:,j);
@@ -196,9 +189,9 @@ else
         M_.params, oo_.dr.ys, 1, dyssdtheta,d2yssdtheta);
     [residual, g1, g2 ] = feval([M_.fname,'_dynamic'],yy0, repmat(oo_.exo_steady_state',[M_.maximum_exo_lag+M_.maximum_exo_lead+1,1]), ...
         M_.params, oo_.dr.ys, 1);
-    [nr, nc]=size(g2);
 end
 
+[nr, nc]=size(g2);
 nc = sqrt(nc);
 Hss = dyssdtheta(oo_.dr.order_var,indx);
 dyssdtheta = dyssdtheta(I,:);
diff --git a/preprocessor/ExprNode.cc b/preprocessor/ExprNode.cc
index 057b780ccf1d9eb5d2d5ee4141ffe4bb51ee3991..3a15342009fe3e41c17b5ebf1f3686b4bd25f626 100644
--- a/preprocessor/ExprNode.cc
+++ b/preprocessor/ExprNode.cc
@@ -215,7 +215,7 @@ ExprNode::createEndoLeadAuxiliaryVarForMyself(subst_table_t &subst_table, vector
       it = subst_table.find(orig_expr);
       if (it == subst_table.end())
         {
-          int symb_id = datatree.symbol_table.addEndoLeadAuxiliaryVar(orig_expr->idx);
+          int symb_id = datatree.symbol_table.addEndoLeadAuxiliaryVar(orig_expr->idx, substexpr);
           neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(datatree.AddVariable(symb_id, 0), substexpr)));
           substexpr = datatree.AddVariable(symb_id, +1);
           assert(dynamic_cast<VariableNode *>(substexpr) != NULL);
@@ -251,7 +251,7 @@ ExprNode::createExoLeadAuxiliaryVarForMyself(subst_table_t &subst_table, vector<
       it = subst_table.find(orig_expr);
       if (it == subst_table.end())
         {
-          int symb_id = datatree.symbol_table.addExoLeadAuxiliaryVar(orig_expr->idx);
+          int symb_id = datatree.symbol_table.addExoLeadAuxiliaryVar(orig_expr->idx, substexpr);
           neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(datatree.AddVariable(symb_id, 0), substexpr)));
           substexpr = datatree.AddVariable(symb_id, +1);
           assert(dynamic_cast<VariableNode *>(substexpr) != NULL);
@@ -502,6 +502,13 @@ NumConstNode::isInStaticForm() const
   return true;
 }
 
+expr_t
+NumConstNode::substituteStaticAuxiliaryVariable() const
+{
+  return const_cast<NumConstNode *>(this);
+}
+
+
 VariableNode::VariableNode(DataTree &datatree_arg, int symb_id_arg, int lag_arg) :
   ExprNode(datatree_arg),
   symb_id(symb_id_arg),
@@ -836,6 +843,22 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
     }
 }
 
+expr_t
+VariableNode::substituteStaticAuxiliaryVariable() const
+{
+  if (type == eEndogenous)
+    {
+      try
+        {
+          return datatree.symbol_table.getAuxiliaryVarsExprNode(symb_id)->substituteStaticAuxiliaryVariable();
+        }
+      catch (SymbolTable::SearchFailedException &e)
+        {
+        }
+    }
+  return const_cast<VariableNode *>(this);
+}
+  
 double
 VariableNode::eval(const eval_context_t &eval_context) const throw (EvalException, EvalExternalFunctionException)
 {
@@ -1217,7 +1240,7 @@ VariableNode::substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector
           it = subst_table.find(orig_expr);
           if (it == subst_table.end())
             {
-              int aux_symb_id = datatree.symbol_table.addEndoLagAuxiliaryVar(symb_id, cur_lag+1);
+              int aux_symb_id = datatree.symbol_table.addEndoLagAuxiliaryVar(symb_id, cur_lag+1, substexpr);
               neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(datatree.AddVariable(aux_symb_id, 0), substexpr)));
               substexpr = datatree.AddVariable(aux_symb_id, -1);
               subst_table[orig_expr] = substexpr;
@@ -1290,7 +1313,7 @@ VariableNode::substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *
           it = subst_table.find(orig_expr);
           if (it == subst_table.end())
             {
-              int aux_symb_id = datatree.symbol_table.addExoLagAuxiliaryVar(symb_id, cur_lag+1);
+              int aux_symb_id = datatree.symbol_table.addExoLagAuxiliaryVar(symb_id, cur_lag+1, substexpr);
               neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(datatree.AddVariable(aux_symb_id, 0), substexpr)));
               substexpr = datatree.AddVariable(aux_symb_id, -1);
               subst_table[orig_expr] = substexpr;
@@ -1339,7 +1362,8 @@ VariableNode::differentiateForwardVars(const vector<string> &subset, subst_table
             diffvar = const_cast<VariableNode *>(it->second);
           else
             {
-              int aux_symb_id = datatree.symbol_table.addDiffForwardAuxiliaryVar(symb_id);
+              int aux_symb_id = datatree.symbol_table.addDiffForwardAuxiliaryVar(symb_id, datatree.AddMinus(datatree.AddVariable(symb_id, 0),
+                                                                                                            datatree.AddVariable(symb_id, -1)));
               neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(datatree.AddVariable(aux_symb_id, 0), datatree.AddMinus(datatree.AddVariable(symb_id, 0),
                                                                                                                                       datatree.AddVariable(symb_id, -1)))));
               diffvar = datatree.AddVariable(aux_symb_id, 1);
@@ -2533,6 +2557,16 @@ UnaryOpNode::isInStaticForm() const
     return arg->isInStaticForm();
 }
 
+expr_t
+UnaryOpNode::substituteStaticAuxiliaryVariable() const
+{
+  expr_t argsubst = arg->substituteStaticAuxiliaryVariable();
+  if (op_code == oExpectation)
+    return argsubst;
+  else
+    return buildSimilarUnaryOpNode(argsubst, datatree);
+}
+
 BinaryOpNode::BinaryOpNode(DataTree &datatree_arg, const expr_t arg1_arg,
                            BinaryOpcode op_code_arg, const expr_t arg2_arg) :
   ExprNode(datatree_arg),
@@ -3820,6 +3854,21 @@ BinaryOpNode::isInStaticForm() const
   return arg1->isInStaticForm() && arg2->isInStaticForm();
 }
 
+expr_t
+BinaryOpNode::substituteStaticAuxiliaryVariable() const
+{
+  expr_t arg1subst = arg1->substituteStaticAuxiliaryVariable();
+  expr_t arg2subst = arg2->substituteStaticAuxiliaryVariable();
+  return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
+}
+
+expr_t
+BinaryOpNode::substituteStaticAuxiliaryDefinition() const
+{
+  expr_t arg2subst = arg2->substituteStaticAuxiliaryVariable();
+  return buildSimilarBinaryOpNode(arg1, arg2subst, datatree);
+}
+
 TrinaryOpNode::TrinaryOpNode(DataTree &datatree_arg, const expr_t arg1_arg,
                              TrinaryOpcode op_code_arg, const expr_t arg2_arg, const expr_t arg3_arg) :
   ExprNode(datatree_arg),
@@ -4471,6 +4520,15 @@ TrinaryOpNode::isInStaticForm() const
   return arg1->isInStaticForm() && arg2->isInStaticForm() && arg3->isInStaticForm();
 }
 
+expr_t
+TrinaryOpNode::substituteStaticAuxiliaryVariable() const
+{
+  expr_t arg1subst = arg1->substituteStaticAuxiliaryVariable();
+  expr_t arg2subst = arg2->substituteStaticAuxiliaryVariable();
+  expr_t arg3subst = arg3->substituteStaticAuxiliaryVariable();
+  return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
+}
+
 AbstractExternalFunctionNode::AbstractExternalFunctionNode(DataTree &datatree_arg,
                                                            int symb_id_arg,
                                                            const vector<expr_t> &arguments_arg) :
@@ -4815,6 +4873,15 @@ AbstractExternalFunctionNode::containsExternalFunction() const
   return true;
 }
 
+expr_t
+AbstractExternalFunctionNode::substituteStaticAuxiliaryVariable() const
+{
+  vector<expr_t> arguments_subst;
+  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
+    arguments_subst.push_back((*it)->substituteStaticAuxiliaryVariable());
+  return buildSimilarExternalFunctionNode(arguments_subst, datatree);
+}
+
 ExternalFunctionNode::ExternalFunctionNode(DataTree &datatree_arg,
                                            int symb_id_arg,
                                            const vector<expr_t> &arguments_arg) :
diff --git a/preprocessor/ExprNode.hh b/preprocessor/ExprNode.hh
index e46f567e49cec1e274b2a4a20b7fcaa6a8e0ae9b..5c0914b9ee5a9d848c5168705316d48de410ccc5 100644
--- a/preprocessor/ExprNode.hh
+++ b/preprocessor/ExprNode.hh
@@ -442,6 +442,9 @@ public:
 
   //! Returns true if the expression is in static form (no lead, no lag, no expectation, no STEADY_STATE)
   virtual bool isInStaticForm() const = 0;
+
+  //! Substitute auxiliary variables by their expression in static model
+  virtual expr_t substituteStaticAuxiliaryVariable() const = 0;
 };
 
 //! Object used to compare two nodes (using their indexes)
@@ -501,6 +504,7 @@ public:
   virtual expr_t cloneDynamic(DataTree &dynamic_datatree) const;
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
   virtual bool isInStaticForm() const;
+  virtual expr_t substituteStaticAuxiliaryVariable() const;
 };
 
 //! Symbol or variable node
@@ -565,6 +569,8 @@ public:
   virtual expr_t cloneDynamic(DataTree &dynamic_datatree) const;
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
   virtual bool isInStaticForm() const;
+  //! Substitute auxiliary variables by their expression in static model
+  virtual expr_t substituteStaticAuxiliaryVariable() const;
 };
 
 //! Unary operator node
@@ -648,6 +654,8 @@ public:
   virtual expr_t cloneDynamic(DataTree &dynamic_datatree) const;
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
   virtual bool isInStaticForm() const;
+  //! Substitute auxiliary variables by their expression in static model
+  virtual expr_t substituteStaticAuxiliaryVariable() const;
 };
 
 //! Binary operator node
@@ -750,6 +758,10 @@ public:
   //! Returns the non-zero hand-side of an equation (that must have a hand side equal to zero)
   expr_t getNonZeroPartofEquation() const;
   virtual bool isInStaticForm() 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
+  virtual expr_t substituteStaticAuxiliaryDefinition() const;
 };
 
 //! Trinary operator node
@@ -820,6 +832,8 @@ public:
   virtual expr_t cloneDynamic(DataTree &dynamic_datatree) const;
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
   virtual bool isInStaticForm() const;
+  //! Substitute auxiliary variables by their expression in static model
+  virtual expr_t substituteStaticAuxiliaryVariable() const;
 };
 
 //! External function node
@@ -899,6 +913,8 @@ public:
   virtual expr_t cloneDynamic(DataTree &dynamic_datatree) const = 0;
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
   virtual bool isInStaticForm() const;
+  //! Substitute auxiliary variables by their expression in static model
+  virtual expr_t substituteStaticAuxiliaryVariable() const;
 };
 
 class ExternalFunctionNode : public AbstractExternalFunctionNode
diff --git a/preprocessor/ModFile.cc b/preprocessor/ModFile.cc
index 00eb9c37564558fd51ce78124e3150e7bee730bf..1776fd9b295341e7c9096b8c452598c07d37e08a 100644
--- a/preprocessor/ModFile.cc
+++ b/preprocessor/ModFile.cc
@@ -39,8 +39,7 @@ ModFile::ModFile(WarningConsolidation &warnings_arg)
     static_model(symbol_table, num_constants, external_functions_table),
     steady_state_model(symbol_table, num_constants, external_functions_table, static_model),
     linear(false), block(false), byte_code(false), use_dll(false), no_static(false), 
-    differentiate_forward_vars(false),
-    nonstationary_variables(false), orig_eqn_nbr(0), ramsey_eqn_nbr(0),
+    differentiate_forward_vars(false), nonstationary_variables(false),
     param_used_with_lead_lag(false), warnings(warnings_arg)
 {
 }
@@ -345,7 +344,7 @@ ModFile::transformPass(bool nostrict)
       dynamic_model.removeTrendVariableFromEquations();
     }
 
-  orig_eqn_nbr = dynamic_model.equation_number();
+  mod_file_struct.orig_eq_nbr = dynamic_model.equation_number();
   if (mod_file_struct.ramsey_model_present)
     {
       StaticModel *planner_objective = NULL;
@@ -364,7 +363,7 @@ ModFile::transformPass(bool nostrict)
       dynamic_model.cloneDynamic(ramsey_FOC_equations_dynamic_model);
       ramsey_FOC_equations_dynamic_model.computeRamseyPolicyFOCs(*planner_objective);
       ramsey_FOC_equations_dynamic_model.replaceMyEquations(dynamic_model);
-      ramsey_eqn_nbr = dynamic_model.equation_number() - orig_eqn_nbr;
+      mod_file_struct.ramsey_eq_nbr = dynamic_model.equation_number() - mod_file_struct.orig_eq_nbr;
     }
 
   if (mod_file_struct.stoch_simul_present
@@ -446,7 +445,7 @@ ModFile::transformPass(bool nostrict)
     cout << "Found " << dynamic_model.equation_number() << " equation(s)." << endl;
   else
     {
-      cout << "Found " << orig_eqn_nbr  << " equation(s)." << endl;
+      cout << "Found " << mod_file_struct.orig_eq_nbr  << " equation(s)." << endl;
       cout << "Found " << dynamic_model.equation_number() << " FOC equation(s) for Ramsey Problem." << endl;
     }
 
@@ -748,9 +747,9 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo
   if (block && !byte_code)
     mOutputFile << "addpath " << basename << ";" << endl;
 
-  mOutputFile << "M_.orig_eq_nbr = " << orig_eqn_nbr << ";" << endl
+  mOutputFile << "M_.orig_eq_nbr = " << mod_file_struct.orig_eq_nbr << ";" << endl
               << "M_.eq_nbr = " << dynamic_model.equation_number() << ";" << endl
-              << "M_.ramsey_eq_nbr = " << ramsey_eqn_nbr << ";" << endl;
+              << "M_.ramsey_eq_nbr = " << mod_file_struct.ramsey_eq_nbr << ";" << endl;
 
   if (dynamic_model.equation_number() > 0)
     {
@@ -1133,9 +1132,9 @@ ModFile::writeExternalFilesJulia(const string &basename, FileOutputType output)
                << symbol_table.exo_nbr() << ")" << endl
                << "model.correlation_matrix = ones(Float64, " << symbol_table.exo_nbr() << ", "
                << symbol_table.exo_nbr() << ")" << endl
-               << "model.orig_eq_nbr = " << orig_eqn_nbr << endl
+               << "model.orig_eq_nbr = " << mod_file_struct.orig_eq_nbr << endl
                << "model.eq_nbr = " << dynamic_model.equation_number() << endl
-               << "model.ramsey_eq_nbr = " << ramsey_eqn_nbr << endl;
+               << "model.ramsey_eq_nbr = " << mod_file_struct.ramsey_eq_nbr << endl;
 
   if (mod_file_struct.calibrated_measurement_errors)
     jlOutputFile << "model.h = zeros(Float64,"
diff --git a/preprocessor/ModFile.hh b/preprocessor/ModFile.hh
index 66c9fb6b995064125eb6899001e590f6692018f5..c720b45e9c56c1aae96d047c35cac8b4110f6b7f 100644
--- a/preprocessor/ModFile.hh
+++ b/preprocessor/ModFile.hh
@@ -101,12 +101,6 @@ public:
   /*! Filled using initval blocks and parameters initializations */
   eval_context_t global_eval_context;
 
-  //! Stores the original number of equations in the model_block
-  int orig_eqn_nbr;
-
-  //! Stores the number of equations added to the Ramsey model
-  int ramsey_eqn_nbr;
-
   //! Parameter used with lead/lag
   bool param_used_with_lead_lag;
 
diff --git a/preprocessor/Statement.hh b/preprocessor/Statement.hh
index 2ac8edbaf531d99b45247a6f40215fee47dae909..050d182a77c9146c5ccd02bcd081f91645bd7e37 100644
--- a/preprocessor/Statement.hh
+++ b/preprocessor/Statement.hh
@@ -91,7 +91,7 @@ public:
   bool bayesian_irf_present;
   //! Whether there is a data statement present
   bool estimation_data_statement_present;
-  //! Last chain number for Markov Switching statement
+  //! Last chain number for Markov Switching statement2
   int last_markov_switching_chain;
   //! Whether a calib_smoother statement is present
   bool calib_smoother_present;
@@ -118,6 +118,11 @@ public:
   bool ms_dsge_present;
   //! Whether occbin is present
   bool occbin_option;
+  //! Stores the original number of equations in the model_block
+  int orig_eq_nbr;
+   //! Stores the number of equations added to the Ramsey model
+  int ramsey_eq_nbr;
+
 };
 
 class Statement
diff --git a/preprocessor/StaticModel.cc b/preprocessor/StaticModel.cc
index 92aa43f64460acab701f03cbb3db4a52fac7f0a4..d6bef86d0da3b2457dc331f0f67b0642c94a4844 100644
--- a/preprocessor/StaticModel.cc
+++ b/preprocessor/StaticModel.cc
@@ -1051,12 +1051,31 @@ StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms
 {
   initializeVariablesAndEquations();
 
+  vector<BinaryOpNode *> neweqs;
+  for (unsigned int eq = 0; eq < equations.size() - aux_equations.size(); eq++)
+    {
+      expr_t eq_tmp = equations[eq]->substituteStaticAuxiliaryVariable();
+      neweqs.push_back(dynamic_cast<BinaryOpNode *>(eq_tmp->toStatic(*this)));
+    }
+
+  for (unsigned int eq = 0; eq < aux_equations.size();  eq++)
+    {
+      expr_t eq_tmp = aux_equations[eq]->substituteStaticAuxiliaryDefinition();
+      neweqs.push_back(dynamic_cast<BinaryOpNode *>(eq_tmp->toStatic(*this)));
+    }
+      
+  equations.clear();
+  copy(neweqs.begin(),neweqs.end(),back_inserter(equations));
   // Compute derivatives w.r. to all endogenous, and possibly exogenous and exogenous deterministic
   set<int> vars;
 
   for (int i = 0; i < symbol_table.endo_nbr(); i++)
-    vars.insert(getDerivID(symbol_table.getID(eEndogenous, i), 0));
-
+    {
+      int id = symbol_table.getID(eEndogenous, i);
+      //      if (!symbol_table.isAuxiliaryVariableButNotMultiplier(id))
+      vars.insert(getDerivID(id, 0));
+    }        
+ 
   // Launch computations
   cout << "Computing static model derivatives:" << endl
        << " - order 1" << endl;
@@ -1076,7 +1095,7 @@ StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms
       computeThirdDerivatives(vars);
     }
 
-if (paramsDerivatives)
+  if (paramsDerivatives)
     {
       cout << " - derivatives of Jacobian/Hessian w.r. to parameters" << endl;
       computeParamsDerivatives();
@@ -1720,7 +1739,7 @@ StaticModel::writeStaticFile(const string &basename, bool block, bool bytecode,
     writeStaticJuliaFile(basename);
   else
     writeStaticMFile(basename);
-  writeAuxVarRecursiveDefinitions(basename, julia);
+  writeSetAuxiliaryVariables(basename, julia);
 }
 
 void
@@ -2080,8 +2099,9 @@ StaticModel::writeAuxVarInitval(ostream &output, ExprNodeOutputType output_type)
     }
 }
 
-void StaticModel::writeAuxVarRecursiveDefinitions(const string &basename, const bool julia) const
+void StaticModel::writeSetAuxiliaryVariables(const string &basename, const bool julia) const
 {
+  
   string func_name = basename + "_set_auxiliary_variables";
   string filename = julia ? func_name + ".jl" : func_name + ".m";
   string comment = julia ? "#" : "%";
@@ -2108,10 +2128,17 @@ void StaticModel::writeAuxVarRecursiveDefinitions(const string &basename, const
     if (dynamic_cast<ExprNode *>(aux_equations[i])->containsExternalFunction())
       dynamic_cast<ExprNode *>(aux_equations[i])->writeExternalFunctionOutput(output, oMatlabStaticModel,
                                                                               temporary_terms, tef_terms);
+  writeAuxVarRecursiveDefinitions(output, oMatlabStaticModel);
+}
 
+void
+StaticModel::writeAuxVarRecursiveDefinitions(ostream &output, ExprNodeOutputType output_type) const
+{
+  deriv_node_temp_terms_t tef_terms;
+  temporary_terms_t temporary_terms;
   for (int i = 0; i < (int) aux_equations.size(); i++)
     {
-      dynamic_cast<ExprNode *>(aux_equations[i])->writeOutput(output, oMatlabStaticModel, temporary_terms, tef_terms);
+      dynamic_cast<ExprNode *>(aux_equations[i]->substituteStaticAuxiliaryDefinition())->writeOutput(output, output_type, temporary_terms, tef_terms);
       output << ";" << endl;
     }
 }
diff --git a/preprocessor/StaticModel.hh b/preprocessor/StaticModel.hh
index 90a3aebb7962c27ec1b8952d8cddf8f17103e1f4..5dfc5dda7247354ed048030e58a397b3e82f92fa 100644
--- a/preprocessor/StaticModel.hh
+++ b/preprocessor/StaticModel.hh
@@ -183,7 +183,8 @@ public:
   void writeAuxVarInitval(ostream &output, ExprNodeOutputType output_type) const;
 
   //! Writes definition of the auxiliary variables in a .m or .jl file
-  void writeAuxVarRecursiveDefinitions(const string &basename, const bool julia) const;
+  void writeSetAuxiliaryVariables(const string &basename, const bool julia) const;
+  void writeAuxVarRecursiveDefinitions(ostream &output, ExprNodeOutputType output_type) const;
 
   virtual int getDerivID(int symb_id, int lag) const throw (UnknownDerivIDException);
   virtual void addAllParamDerivId(set<int> &deriv_id_set);
diff --git a/preprocessor/SteadyStateModel.cc b/preprocessor/SteadyStateModel.cc
index 4533c41f27a2e232d2a8c77c6cd0a1d703a222db..61220a80509363cb638954aaa53a584a283074e8 100644
--- a/preprocessor/SteadyStateModel.cc
+++ b/preprocessor/SteadyStateModel.cc
@@ -161,7 +161,7 @@ SteadyStateModel::writeSteadyStateFile(const string &basename, bool ramsey_model
     output << "    % Auxiliary equations" << endl;
   else
     output << "    # Auxiliary equations" << endl;
-  static_model.writeAuxVarInitval(output, output_type);
+  static_model.writeAuxVarRecursiveDefinitions(output, output_type);
 
   if (!julia)
     output << "    check_=0;" << endl;
diff --git a/preprocessor/SymbolTable.cc b/preprocessor/SymbolTable.cc
index 73433ae9e3c7d078f81f7dfb13af88080ba285ce..0a15fc070e1bdd60c21495d8bcf2477bc9585ecc 100644
--- a/preprocessor/SymbolTable.cc
+++ b/preprocessor/SymbolTable.cc
@@ -26,14 +26,14 @@
 
 AuxVarInfo::AuxVarInfo(int symb_id_arg, aux_var_t type_arg, int orig_symb_id_arg, int orig_lead_lag_arg,
                        int equation_number_for_multiplier_arg, int information_set_arg,
-                       expr_t expectation_expr_node_arg) :
+                       expr_t expr_node_arg) :
   symb_id(symb_id_arg),
   type(type_arg),
   orig_symb_id(orig_symb_id_arg),
   orig_lead_lag(orig_lead_lag_arg),
   equation_number_for_multiplier(equation_number_for_multiplier_arg),
   information_set(information_set_arg),
-  expectation_expr_node(expectation_expr_node_arg)
+  expr_node(expr_node_arg)
 {
 }
 
@@ -333,7 +333,7 @@ SymbolTable::writeOutput(ostream &output) const throw (NotYetFrozenException)
             output << "M_.aux_vars(" << i+1 << ").orig_expr = '\\mathbb{E}_{t"
                    << (aux_vars[i].get_information_set() < 0 ? "" : "+")
                    << aux_vars[i].get_information_set() << "}(";
-            aux_vars[i].get_expectation_expr_node()->writeOutput(output, oLatexDynamicModel);
+            aux_vars[i].get_expr_node()->writeOutput(output, oLatexDynamicModel);
             output << ")';" << endl;
             break;
           }
@@ -520,7 +520,7 @@ SymbolTable::writeCCOutput(ostream &output) const throw (NotYetFrozenException)
 }
 
 int
-SymbolTable::addLeadAuxiliaryVarInternal(bool endo, int index) throw (FrozenException)
+SymbolTable::addLeadAuxiliaryVarInternal(bool endo, int index, expr_t expr_arg) throw (FrozenException)
 {
   ostringstream varname;
   if (endo)
@@ -539,13 +539,13 @@ SymbolTable::addLeadAuxiliaryVarInternal(bool endo, int index) throw (FrozenExce
       exit(EXIT_FAILURE);
     }
 
-  aux_vars.push_back(AuxVarInfo(symb_id, (endo ? avEndoLead : avExoLead), 0, 0, 0, 0, NULL));
+  aux_vars.push_back(AuxVarInfo(symb_id, (endo ? avEndoLead : avExoLead), 0, 0, 0, 0, expr_arg));
 
   return symb_id;
 }
 
 int
-SymbolTable::addLagAuxiliaryVarInternal(bool endo, int orig_symb_id, int orig_lead_lag) throw (FrozenException)
+SymbolTable::addLagAuxiliaryVarInternal(bool endo, int orig_symb_id, int orig_lead_lag, expr_t expr_arg) throw (FrozenException)
 {
   ostringstream varname;
   if (endo)
@@ -565,37 +565,37 @@ SymbolTable::addLagAuxiliaryVarInternal(bool endo, int orig_symb_id, int orig_le
       exit(EXIT_FAILURE);
     }
 
-  aux_vars.push_back(AuxVarInfo(symb_id, (endo ? avEndoLag : avExoLag), orig_symb_id, orig_lead_lag, 0, 0, NULL));
+  aux_vars.push_back(AuxVarInfo(symb_id, (endo ? avEndoLag : avExoLag), orig_symb_id, orig_lead_lag, 0, 0, expr_arg));
 
   return symb_id;
 }
 
 int
-SymbolTable::addEndoLeadAuxiliaryVar(int index) throw (FrozenException)
+SymbolTable::addEndoLeadAuxiliaryVar(int index, expr_t expr_arg) throw (FrozenException)
 {
-  return addLeadAuxiliaryVarInternal(true, index);
+  return addLeadAuxiliaryVarInternal(true, index, expr_arg);
 }
 
 int
-SymbolTable::addEndoLagAuxiliaryVar(int orig_symb_id, int orig_lead_lag) throw (FrozenException)
+SymbolTable::addEndoLagAuxiliaryVar(int orig_symb_id, int orig_lead_lag, expr_t expr_arg) throw (FrozenException)
 {
-  return addLagAuxiliaryVarInternal(true, orig_symb_id, orig_lead_lag);
+  return addLagAuxiliaryVarInternal(true, orig_symb_id, orig_lead_lag, expr_arg);
 }
 
 int
-SymbolTable::addExoLeadAuxiliaryVar(int index) throw (FrozenException)
+SymbolTable::addExoLeadAuxiliaryVar(int index, expr_t expr_arg) throw (FrozenException)
 {
-  return addLeadAuxiliaryVarInternal(false, index);
+  return addLeadAuxiliaryVarInternal(false, index, expr_arg);
 }
 
 int
-SymbolTable::addExoLagAuxiliaryVar(int orig_symb_id, int orig_lead_lag) throw (FrozenException)
+SymbolTable::addExoLagAuxiliaryVar(int orig_symb_id, int orig_lead_lag, expr_t expr_arg) throw (FrozenException)
 {
-  return addLagAuxiliaryVarInternal(false, orig_symb_id, orig_lead_lag);
+  return addLagAuxiliaryVarInternal(false, orig_symb_id, orig_lead_lag, expr_arg);
 }
 
 int
-SymbolTable::addExpectationAuxiliaryVar(int information_set, int index, expr_t exp_arg) throw (FrozenException)
+SymbolTable::addExpectationAuxiliaryVar(int information_set, int index, expr_t expr_arg) throw (FrozenException)
 {
   ostringstream varname;
   int symb_id;
@@ -613,7 +613,7 @@ SymbolTable::addExpectationAuxiliaryVar(int information_set, int index, expr_t e
       exit(EXIT_FAILURE);
     }
 
-  aux_vars.push_back(AuxVarInfo(symb_id, avExpectation, 0, 0, 0, information_set, exp_arg));
+  aux_vars.push_back(AuxVarInfo(symb_id, avExpectation, 0, 0, 0, information_set, expr_arg));
 
   return symb_id;
 }
@@ -640,7 +640,7 @@ SymbolTable::addMultiplierAuxiliaryVar(int index) throw (FrozenException)
 }
 
 int
-SymbolTable::addDiffForwardAuxiliaryVar(int orig_symb_id) throw (FrozenException)
+SymbolTable::addDiffForwardAuxiliaryVar(int orig_symb_id, expr_t expr_arg) throw (FrozenException)
 {
   ostringstream varname;
   int symb_id;
@@ -656,7 +656,7 @@ SymbolTable::addDiffForwardAuxiliaryVar(int orig_symb_id) throw (FrozenException
       exit(EXIT_FAILURE);
     }
 
-  aux_vars.push_back(AuxVarInfo(symb_id, avDiffForward, orig_symb_id, 0, 0, 0, NULL));
+  aux_vars.push_back(AuxVarInfo(symb_id, avDiffForward, orig_symb_id, 0, 0, 0, expr_arg));
   return symb_id;
 }
 
@@ -670,6 +670,22 @@ SymbolTable::searchAuxiliaryVars(int orig_symb_id, int orig_lead_lag) const thro
   throw SearchFailedException(orig_symb_id, orig_lead_lag);
 }
 
+expr_t
+SymbolTable::getAuxiliaryVarsExprNode(int symb_id) const throw (SearchFailedException)
+// throw exception if it is a Lagrange multiplier
+{
+  for (size_t i = 0; i < aux_vars.size(); i++)
+    if (aux_vars[i].get_symb_id() == symb_id)
+      {
+        expr_t expr_node = aux_vars[i].get_expr_node();
+        if (expr_node != NULL)
+          return expr_node;
+        else
+          throw SearchFailedException(symb_id);
+      }
+  throw SearchFailedException(symb_id);
+}
+
 void
 SymbolTable::markPredetermined(int symb_id) throw (UnknownSymbolIDException, FrozenException)
 {
@@ -770,6 +786,15 @@ SymbolTable::isAuxiliaryVariable(int symb_id) const
   return false;
 }
 
+bool
+SymbolTable::isAuxiliaryVariableButNotMultiplier(int symb_id) const
+{
+  for (int i = 0; i < (int) aux_vars.size(); i++)
+    if (aux_vars[i].get_symb_id() == symb_id && aux_vars[i].get_type() != avMultiplier)
+      return true;
+  return false;
+}
+
 set<int>
 SymbolTable::getOrigEndogenous() const
 {
@@ -862,7 +887,7 @@ SymbolTable::writeJuliaOutput(ostream &output) const throw (NotYetFrozenExceptio
               output << "NaN, NaN, NaN, \"\\mathbb{E}_{t"
                      << (aux_vars[i].get_information_set() < 0 ? "" : "+")
                      << aux_vars[i].get_information_set() << "}(";
-              aux_vars[i].get_expectation_expr_node()->writeOutput(output, oLatexDynamicModel);
+              aux_vars[i].get_expr_node()->writeOutput(output, oLatexDynamicModel);
               output << ")\"";
               break;
             }
diff --git a/preprocessor/SymbolTable.hh b/preprocessor/SymbolTable.hh
index f212953f4de392bb68f2237c27d1f14b144918fe..4f623dfdd6416bb912f48c676426a0f265724e69 100644
--- a/preprocessor/SymbolTable.hh
+++ b/preprocessor/SymbolTable.hh
@@ -55,16 +55,16 @@ private:
   int orig_lead_lag; //!< Lead/lag of the endo of the original model represented by this aux var. Only used for avEndoLag and avExoLag.
   int equation_number_for_multiplier; //!< Stores the original constraint equation number associated with this aux var. Only used for avMultiplier.
   int information_set; //! Argument of expectation operator. Only used for avExpectation.
-  expr_t expectation_expr_node; //! Argument of expectation operator. Only used for avExpectation.
+  expr_t expr_node; //! Auxiliary variable definition
 public:
-  AuxVarInfo(int symb_id_arg, aux_var_t type_arg, int orig_symb_id, int orig_lead_lag, int equation_number_for_multiplier_arg, int information_set_arg, expr_t expectation_expr_node_arg);
+  AuxVarInfo(int symb_id_arg, aux_var_t type_arg, int orig_symb_id, int orig_lead_lag, int equation_number_for_multiplier_arg, int information_set_arg, expr_t expr_node_arg);
   int get_symb_id() const { return symb_id; };
   aux_var_t get_type() const { return type; };
   int get_orig_symb_id() const { return orig_symb_id; };
   int get_orig_lead_lag() const { return orig_lead_lag; };
   int get_equation_number_for_multiplier() const { return equation_number_for_multiplier; };
   int get_information_set() const { return information_set; };
-  expr_t get_expectation_expr_node() const { return expectation_expr_node; } ;
+  expr_t get_expr_node() const { return expr_node; } ;
 };
 
 //! Stores the symbol table
@@ -177,18 +177,21 @@ public:
   class SearchFailedException
   {
   public:
-    int orig_symb_id, orig_lead_lag;
+    int orig_symb_id, orig_lead_lag, symb_id;
     SearchFailedException(int orig_symb_id_arg, int orig_lead_lag_arg) : orig_symb_id(orig_symb_id_arg),
                                                                          orig_lead_lag(orig_lead_lag_arg)
     {
     }
+    SearchFailedException(int symb_id_arg) : symb_id(symb_id_arg)
+    {
+    }
   };
 
 private:
   //! Factorized code for adding aux lag variables
-  int addLagAuxiliaryVarInternal(bool endo, int orig_symb_id, int orig_lead_lag) throw (FrozenException);
+  int addLagAuxiliaryVarInternal(bool endo, int orig_symb_id, int orig_lead_lag, expr_t arg) throw (FrozenException);
   //! Factorized code for adding aux lead variables
-  int addLeadAuxiliaryVarInternal(bool endo, int index) throw (FrozenException);
+  int addLeadAuxiliaryVarInternal(bool endo, int index, expr_t arg) throw (FrozenException);
 
 public:
   //! Add a symbol
@@ -201,24 +204,24 @@ public:
   /*!
     \param[in] index Used to construct the variable name
     \return the symbol ID of the new symbol */
-  int addEndoLeadAuxiliaryVar(int index) throw (FrozenException);
+  int addEndoLeadAuxiliaryVar(int index, expr_t arg) throw (FrozenException);
   //! Adds an auxiliary variable for endogenous with lag >= 2
   /*!
     \param[in] orig_symb_id symbol ID of the endogenous declared by the user that this new variable will represent
     \param[in] orig_lead_lag lag value such that this new variable will be equivalent to orig_symb_id(orig_lead_lag)
     \return the symbol ID of the new symbol */
-  int addEndoLagAuxiliaryVar(int orig_symb_id, int orig_lead_lag) throw (FrozenException);
+  int addEndoLagAuxiliaryVar(int orig_symb_id, int orig_lead_lag, expr_t arg) throw (FrozenException);
   //! Adds an auxiliary variable for endogenous with lead >= 1
   /*!
     \param[in] index Used to construct the variable name
     \return the symbol ID of the new symbol */
-  int addExoLeadAuxiliaryVar(int index) throw (FrozenException);
+  int addExoLeadAuxiliaryVar(int index, expr_t arg) throw (FrozenException);
   //! Adds an auxiliary variable for exogenous with lag >= 1
   /*!
     \param[in] orig_symb_id symbol ID of the exogenous declared by the user that this new variable will represent
     \param[in] orig_lead_lag lag value such that this new variable will be equivalent to orig_symb_id(orig_lead_lag)
     \return the symbol ID of the new symbol */
-  int addExoLagAuxiliaryVar(int orig_symb_id, int orig_lead_lag) throw (FrozenException);
+  int addExoLagAuxiliaryVar(int orig_symb_id, int orig_lead_lag, expr_t arg) throw (FrozenException);
   //! Adds an auxiliary variable for the expectation operator
   /*!
     \param[in] information_set information set (possibly negative) of the expectation operator
@@ -237,7 +240,7 @@ public:
     \param[in] orig_symb_id The symb_id of the forward variable
     \return the symbol ID of the new symbol
   */
-  int addDiffForwardAuxiliaryVar(int orig_symb_id) throw (FrozenException);
+  int addDiffForwardAuxiliaryVar(int orig_symb_id, expr_t arg) throw (FrozenException);
   //! Searches auxiliary variables which are substitutes for a given symbol_id and lead/lag
   /*!
     The search is only performed among auxiliary variables of endo/exo lag.
@@ -247,6 +250,8 @@ public:
   int searchAuxiliaryVars(int orig_symb_id, int orig_lead_lag) const throw (SearchFailedException);
   //! Returns the number of auxiliary variables
   int AuxVarsSize() const { return aux_vars.size(); };
+  //! Retruns expr_node for an auxiliary variable
+  expr_t getAuxiliaryVarsExprNode(int symb_id) const throw (SearchFailedException);
   //! Tests if symbol already exists
   inline bool exists(const string &name) const;
   //! Get symbol name (by ID)
@@ -320,6 +325,8 @@ public:
   set <int> getEndogenous() const;
   //! Is a given symbol an auxiliary variable
   bool isAuxiliaryVariable(int symb_id) const;
+  //! Is a given symbol an auxiliary variable but not a Lagrange multiplier
+  bool isAuxiliaryVariableButNotMultiplier(int symb_id) const;
   //! Get list of endogenous variables without aux vars
   set <int> getOrigEndogenous() const;
 };