diff --git a/CodeInterpreter.hh b/CodeInterpreter.hh
index f1af1d8a257b4cd14fac6c1268b6c1299b2f57b9..4a289bf89abef7473570e4ec8711bbb590a45508 100644
--- a/CodeInterpreter.hh
+++ b/CodeInterpreter.hh
@@ -167,7 +167,8 @@ enum BinaryOpcode
     oLessEqual,
     oGreaterEqual,
     oEqualEqual,
-    oDifferent
+    oDifferent,
+    oExpectation
   };
 
 enum TrinaryOpcode
@@ -180,8 +181,7 @@ enum TrinaryOpcode
     int Equation, Variable, Own_Derivative;
   };
 
-#pragma pack(push)
-#pragma pack(1)
+#pragma pack(push, 1)
 class TagWithoutArgument
 {
 protected:
diff --git a/DataTree.cc b/DataTree.cc
index 34c97ddb6d70b7545c434cc4d2bf555752c6a70f..1ff69f76b0872b39529bc11997acb68973847d35 100644
--- a/DataTree.cc
+++ b/DataTree.cc
@@ -244,6 +244,17 @@ DataTree::AddPower(NodeID iArg1, NodeID iArg2)
     return Zero;
 }
 
+NodeID
+DataTree::AddExpectation(int iArg1, NodeID iArg2)
+{
+  ostringstream period;
+  period << abs(iArg1);
+  if (iArg1 >= 0)
+    return AddBinaryOp(AddNumConstant(period.str()), oExpectation, iArg2);
+  else
+    return AddBinaryOp(AddUMinus(AddNumConstant(period.str())), oExpectation, iArg2);
+}
+
 NodeID
 DataTree::AddExp(NodeID iArg1)
 {
diff --git a/DataTree.hh b/DataTree.hh
index f9cfd42cdf37876f4e531fadfe7d90379423cc0b..c2dbe4789b07de3452183a7de4c088b04d84eb4f 100644
--- a/DataTree.hh
+++ b/DataTree.hh
@@ -128,6 +128,8 @@ public:
   NodeID AddDifferent(NodeID iArg1, NodeID iArg2);
   //! Adds "arg1^arg2" to model tree
   NodeID AddPower(NodeID iArg1, NodeID iArg2);
+  //! Adds "E(arg1)(arg2)" to model tree
+  NodeID AddExpectation(int iArg1, NodeID iArg2);
   //! Adds "exp(arg)" to model tree
   NodeID AddExp(NodeID iArg1);
   //! Adds "log(arg)" to model tree
diff --git a/DynamicModel.cc b/DynamicModel.cc
index 3a53492728bea784c2ac06ee556e292dfe873997..1c8e9deeb3cc62620ec87e105f5d9f27a8d6ae89 100644
--- a/DynamicModel.cc
+++ b/DynamicModel.cc
@@ -1201,7 +1201,7 @@ DynamicModel::writeDynamicMFile(const string &dynamic_basename) const
     << "%" << endl
     << "% Warning : this file is generated automatically by Dynare" << endl
     << "%           from model file (.mod)" << endl << endl;
-
+    
     if (containsSteadyStateOperator())
       mDynamicModelFile << "global oo_;" << endl << endl;
 
@@ -1215,7 +1215,7 @@ DynamicModel::writeDynamicCFile(const string &dynamic_basename) const
   {
     string filename = dynamic_basename + ".c";
     ofstream mDynamicModelFile;
-
+    
     mDynamicModelFile.open(filename.c_str(), ios::out | ios::binary);
     if (!mDynamicModelFile.is_open())
       {
@@ -1223,14 +1223,17 @@ DynamicModel::writeDynamicCFile(const string &dynamic_basename) const
         exit(EXIT_FAILURE);
       }
     mDynamicModelFile << "/*" << endl
-    << " * " << filename << " : Computes dynamic model for Dynare" << endl
-    << " *" << endl
-    << " * Warning : this file is generated automatically by Dynare" << endl
-    << " *           from model file (.mod)" << endl
-    << endl
-    << " */" << endl
-    << "#include <math.h>" << endl
-    << "#include \"mex.h\"" << endl;
+                      << " * " << filename << " : Computes dynamic model for Dynare" << endl
+                      << " *" << endl
+                      << " * Warning : this file is generated automatically by Dynare" << endl
+                      << " *           from model file (.mod)" << endl
+                      << endl
+                      << " */" << endl
+                      << "#include <math.h>" << endl
+                      << "#include \"mex.h\"" << endl
+                      << endl
+                      << "#define max(a, b) (((a) > (b)) ? (a) : (b))" << endl
+                      << "#define min(a, b) (((a) > (b)) ? (b) : (a))" << endl;
 
     // Writing the function body
     writeDynamicModel(mDynamicModelFile, true);
@@ -1789,11 +1792,9 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll) const
         int var = it->first.second;
         NodeID d1 = it->second;
 
-        ostringstream g1;
-        g1 << "  g1";
-        jacobianHelper(g1, eq, getDynJacobianCol(var), output_type);
-
-        jacobian_output << g1.str() << "=" << g1.str() << "+";
+        jacobian_output << "g1";
+        jacobianHelper(jacobian_output, eq, getDynJacobianCol(var), output_type);
+        jacobian_output << "=";
         d1->writeOutput(jacobian_output, output_type, temporary_terms);
         jacobian_output << ";" << endl;
       }
@@ -2912,7 +2913,7 @@ DynamicModel::jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOut
 {
   output << LEFT_ARRAY_SUBSCRIPT(output_type);
   if (IS_MATLAB(output_type))
-    output << eq_nb + 1 << ", " << col_nb + 1;
+    output << eq_nb + 1 << "," << col_nb + 1;
   else
     output << eq_nb + col_nb * equations.size();
   output << RIGHT_ARRAY_SUBSCRIPT(output_type);
@@ -2923,7 +2924,7 @@ DynamicModel::sparseHelper(int order, ostream &output, int row_nb, int col_nb, E
 {
   output << "v" << order << LEFT_ARRAY_SUBSCRIPT(output_type);
   if (IS_MATLAB(output_type))
-    output << row_nb + 1 << ", " << col_nb + 1;
+    output << row_nb + 1 << "," << col_nb + 1;
   else
     output << row_nb + col_nb * NNZDerivatives[order-1];
   output << RIGHT_ARRAY_SUBSCRIPT(output_type);
@@ -3035,6 +3036,39 @@ DynamicModel::substituteLeadLagInternal(aux_var_t type)
     }
 }
 
+void
+DynamicModel::substituteExpectation(bool partial_information_model)
+{
+  ExprNode::subst_table_t subst_table;
+  vector<BinaryOpNode *> neweqs;
+
+  // Substitute in model binary op node map
+  for(binary_op_node_map_type::reverse_iterator it = binary_op_node_map.rbegin();
+      it != binary_op_node_map.rend(); it++)
+    it->second->substituteExpectation(subst_table, neweqs, partial_information_model);
+
+  // Substitute in equations
+  for(int i = 0; i < (int) equations.size(); i++)
+    {
+      BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(equations[i]->substituteExpectation(subst_table, neweqs, partial_information_model));
+      assert(substeq != NULL);
+      equations[i] = substeq;
+    }
+
+  // Add new equations
+  for(int i = 0; i < (int) neweqs.size(); i++)
+    addEquation(neweqs[i]);
+
+  // Add the new set of equations at the *beginning* of aux_equations
+  copy(neweqs.rbegin(), neweqs.rend(), front_inserter(aux_equations));
+
+  if (neweqs.size() > 0)
+    if (partial_information_model)
+      cout << "Substitution of Expectation operator: added auxiliary variables and equations." << endl; //FIX to reflect correct number of equations
+    else
+      cout << "Substitution of Expectation operator: added " << neweqs.size() << " auxiliary variables and equations." << endl;
+}
+
 void
 DynamicModel::fillEvalContext(eval_context_type &eval_context) const
 {
diff --git a/DynamicModel.hh b/DynamicModel.hh
index 9ec0f5f6cfb03b7622cd73bac12a5f165f272af9..a17fab51cbf69a29ad1e75adb8dd746824abf4e5 100644
--- a/DynamicModel.hh
+++ b/DynamicModel.hh
@@ -217,6 +217,9 @@ public:
   //! Transforms the model by removing all lags on exos
   void substituteExoLag();
 
+  //! Transforms the model by removing all oExpectation
+  void substituteExpectation(bool partial_information_model);
+
   //! Fills eval context with values of model local variables and auxiliary variables
   void fillEvalContext(eval_context_type &eval_context) const;
 };
diff --git a/DynareBison.yy b/DynareBison.yy
index 026124fa5e641c40a7e6732e87d6861264f118ec..74afd2d00889839fa94e28637996a9cf60bcc873 100644
--- a/DynareBison.yy
+++ b/DynareBison.yy
@@ -98,7 +98,7 @@ class ParsingDriver;
 %token FILENAME FILTER_STEP_AHEAD FILTERED_VARS FIRST_OBS
 %token <string_val> FLOAT_NUMBER
 %token FORECAST
-%token GAMMA_PDF GRAPH
+%token GAMMA_PDF GRAPH CONDITIONAL_VARIANCE_DECOMPOSITION
 %token HISTVAL HOMOTOPY_SETUP HOMOTOPY_MODE HOMOTOPY_STEPS HP_FILTER HP_NGRID
 %token IDENTIFICATION INF_CONSTANT INITVAL INITVAL_FILE
 %token <string_val> INT_NUMBER
@@ -134,7 +134,7 @@ class ParsingDriver;
 %left UMINUS UPLUS
 %nonassoc POWER
 %token EXP LOG LN LOG10 SIN COS TAN ASIN ACOS ATAN SINH COSH TANH
-%token ASINH ACOSH ATANH SQRT NORMCDF STEADY_STATE
+%token ASINH ACOSH ATANH SQRT NORMCDF STEADY_STATE EXPECTATION
 /* GSA analysis */
 %token DYNARE_SENSITIVITY MORRIS STAB REDFORM PPRIOR PRIOR_RANGE PPOST ILPTAU GLUE MORRIS_NLIV
 %token MORRIS_NTRA NSAM LOAD_REDFORM LOAD_RMSE LOAD_STAB ALPHA2_STAB KSSTAT LOGTRANS_REDFORM THRESHOLD_REDFORM
@@ -364,6 +364,8 @@ expression : '(' expression ')'
              { $$ = driver.add_equal_equal($1, $3); }
            | expression EXCLAMATION_EQUAL expression
              { $$ = driver.add_different($1, $3); }
+           | EXPECTATION '(' signed_integer ')''(' expression ')'
+	     { $$ = driver.add_expectation($3, $6); }
            | MINUS expression %prec UMINUS
              { $$ = driver.add_uminus($2); }
            | PLUS expression %prec UPLUS
@@ -512,6 +514,8 @@ hand_side : '(' hand_side ')'
             { $$ = driver.add_different($1, $3); }
           | hand_side POWER hand_side
             { $$ = driver.add_power($1, $3); }
+          | EXPECTATION '(' signed_integer ')''(' hand_side ')'
+	    { $$ = driver.add_expectation($3, $6); }
           | MINUS hand_side %prec UMINUS
             { $$ = driver.add_uminus($2); }
           | PLUS hand_side
@@ -561,14 +565,11 @@ model_var : symbol
 
 shocks : SHOCKS ';' shock_list END { driver.end_shocks(); };
 
-mshocks : MSHOCKS ';' shock_list END { driver.end_mshocks(); };
-
 shock_list : shock_list shock_elem
            | shock_elem
            ;
 
-shock_elem : VAR symbol ';' PERIODS period_list ';' VALUES value_list ';'
-             { driver.add_det_shock($2, false); }
+shock_elem : det_shock_elem
            | VAR symbol ';' STDERR expression ';'
              { driver.add_stderr_shock($2, $5); }
            | VAR symbol EQUAL expression ';'
@@ -579,6 +580,16 @@ shock_elem : VAR symbol ';' PERIODS period_list ';' VALUES value_list ';'
              { driver.add_correl_shock($2, $4, $6); }
            ;
 
+det_shock_elem : VAR symbol ';' PERIODS period_list ';' VALUES value_list ';'
+                 { driver.add_det_shock($2, false); }
+               ;
+
+mshocks : MSHOCKS ';' mshock_list END { driver.end_mshocks(); };
+
+mshock_list : mshock_list det_shock_elem
+            | det_shock_elem
+            ;
+
 period_list : period_list COMMA INT_NUMBER
               { driver.add_period($3); }
             | period_list INT_NUMBER
@@ -718,6 +729,7 @@ stoch_simul_options : o_dr_algo
                     | o_noprint
                     | o_aim_solver
                     | o_partial_information
+                    | o_conditional_variance_decomposition
                     ;
 
 symbol_list : symbol_list symbol
@@ -1548,9 +1560,14 @@ o_dr_algo : DR_ALGO EQUAL INT_NUMBER {
                                          driver.warning("dr_algo option is now deprecated, and may be removed in a future version of Dynare");
                                        else
                                          driver.error("dr_algo=1 option is no longer supported");
-                                     }
+                                     };
 o_solve_algo : SOLVE_ALGO EQUAL INT_NUMBER { driver.option_num("solve_algo", $3); };
-o_simul_algo : SIMUL_ALGO EQUAL INT_NUMBER { driver.option_num("simul_algo", $3); };
+o_simul_algo : SIMUL_ALGO EQUAL INT_NUMBER {
+                                             if (*$3 == string("0"))
+                                               driver.warning("simul_algo option is now deprecated, and may be removed in a future version of Dynare");
+                                             else
+                                               driver.error("simul_algo=1 option is no longer supported");
+                                           };
 o_stack_solve_algo : STACK_SOLVE_ALGO EQUAL INT_NUMBER { driver.option_num("stack_solve_algo", $3); };
 o_linear : LINEAR { driver.linear(); };
 o_order : ORDER EQUAL INT_NUMBER { driver.option_num("order", $3); };
@@ -1578,6 +1595,11 @@ o_nobs : NOBS EQUAL vec_int
        | NOBS EQUAL INT_NUMBER
          { driver.option_num("nobs", $3); }
        ;
+o_conditional_variance_decomposition : CONDITIONAL_VARIANCE_DECOMPOSITION EQUAL vec_int
+                                       { driver.option_num("conditional_variance_decomposition", $3); }
+                                     | CONDITIONAL_VARIANCE_DECOMPOSITION EQUAL INT_NUMBER
+                                       { driver.option_num("conditional_variance_decomposition", $3); }
+                                     ; 
 o_first_obs : FIRST_OBS EQUAL INT_NUMBER { driver.option_num("first_obs", $3); };
 o_prefilter : PREFILTER EQUAL INT_NUMBER { driver.option_num("prefilter", $3); };
 o_presample : PRESAMPLE EQUAL INT_NUMBER { driver.option_num("presample", $3); };
diff --git a/DynareFlex.ll b/DynareFlex.ll
index c871c7987d7f4d4836604fcc1ebf7ec741a623c1..35bf2b8480f043c91853defd2ee8fff4477c1ecf 100644
--- a/DynareFlex.ll
+++ b/DynareFlex.ll
@@ -223,6 +223,7 @@ int sigma_e = 0;
 <DYNARE_STATEMENT>plot_priors   {return token::PLOT_PRIORS;}
 <DYNARE_STATEMENT>aim_solver {return token::AIM_SOLVER;}
 <DYNARE_STATEMENT>partial_information {return token::PARTIAL_INFORMATION;}
+<DYNARE_STATEMENT>conditional_variance_decomposition {return token::CONDITIONAL_VARIANCE_DECOMPOSITION;}
 
 <DYNARE_STATEMENT>freq {return token::FREQ;}
 <DYNARE_STATEMENT>initial_year {return token::INITIAL_YEAR;}
@@ -439,6 +440,7 @@ int sigma_e = 0;
 <DYNARE_STATEMENT,DYNARE_BLOCK>min {return token::MIN;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>normcdf {return token::NORMCDF;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>steady_state {return token::STEADY_STATE;}
+<DYNARE_STATEMENT,DYNARE_BLOCK>expectation {return token::EXPECTATION;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>nan {return token::NAN_CONSTANT;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>inf {return token::INF_CONSTANT;}
 
diff --git a/ExprNode.cc b/ExprNode.cc
index e2fd293dc0512f5454edec317d2028a0e49ca4bb..96adcdefa6509c2330e2d3550d1f2d11edaf337b 100644
--- a/ExprNode.cc
+++ b/ExprNode.cc
@@ -33,6 +33,7 @@ using namespace __gnu_cxx;
 #include "ExprNode.hh"
 #include "DataTree.hh"
 #include "BlockTriangular.hh"
+#include "ModFile.hh"
 
 ExprNode::ExprNode(DataTree &datatree_arg) : datatree(datatree_arg), preparedForDerivation(false)
 {
@@ -219,7 +220,6 @@ ExprNode::createExoLeadAuxiliaryVarForMyself(subst_table_t &subst_table, vector<
   return dynamic_cast<VariableNode *>(substexpr);
 }
 
-
 NumConstNode::NumConstNode(DataTree &datatree_arg, int id_arg) :
     ExprNode(datatree_arg),
     id(id_arg)
@@ -342,6 +342,12 @@ NumConstNode::substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *
   return const_cast<NumConstNode *>(this);
 }
 
+NodeID
+NumConstNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) 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),
@@ -954,6 +960,12 @@ VariableNode::substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *
     }
 }
 
+NodeID
+VariableNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const
+{
+  return const_cast<VariableNode *>(this);
+}
+
 UnaryOpNode::UnaryOpNode(DataTree &datatree_arg, UnaryOpcode op_code_arg, const NodeID arg_arg) :
     ExprNode(datatree_arg),
     arg(arg_arg),
@@ -1360,13 +1372,20 @@ UnaryOpNode::eval_opcode(UnaryOpcode op_code, double v) throw (EvalException)
       return(sinh(v));
     case oTanh:
       return(tanh(v));
-#ifndef WIN64
+#ifndef _WIN64
     case oAcosh:
       return(acosh(v));
     case oAsinh:
       return(asinh(v));
     case oAtanh:
       return(atanh(v));
+#else
+    case oAcosh:
+      throw EvalException();
+    case oAsinh:
+      throw EvalException();
+    case oAtanh:
+      throw EvalException();
 #endif
     case oSqrt:
       return(sqrt(v));
@@ -1648,6 +1667,13 @@ UnaryOpNode::substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *>
   return buildSimilarUnaryOpNode(argsubst, datatree);
 }
 
+NodeID
+UnaryOpNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const
+{
+  NodeID argsubst = arg->substituteExpectation(subst_table, neweqs, partial_information_model);
+  return buildSimilarUnaryOpNode(argsubst, datatree);
+}
+
 BinaryOpNode::BinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg,
                            BinaryOpcode op_code_arg, const NodeID arg2_arg) :
     ExprNode(datatree_arg),
@@ -1794,6 +1820,7 @@ BinaryOpNode::precedence(ExprNodeOutputType output_type, const temporary_terms_t
           return 5;
       case oMin:
       case oMax:
+      case oExpectation:
         return 100;
       }
     // Suppress GCC warning
@@ -1834,6 +1861,7 @@ BinaryOpNode::cost(const temporary_terms_type &temporary_terms, bool is_matlab)
         case oPower:
           return cost + 1160;
         case oEqual:
+        case oExpectation:
           return cost;
         }
     else
@@ -1859,6 +1887,7 @@ BinaryOpNode::cost(const temporary_terms_type &temporary_terms, bool is_matlab)
         case oPower:
           return cost + 520;
         case oEqual:
+        case oExpectation:
           return cost;
         }
     // Suppress GCC warning
@@ -1957,6 +1986,7 @@ BinaryOpNode::eval_opcode(double v1, BinaryOpcode op_code, double v2) throw (Eva
     case oDifferent:
       return (v1 != v2);
     case oEqual:
+    case oExpectation:
       throw EvalException();
     }
   // Suppress GCC warning
@@ -2494,6 +2524,8 @@ BinaryOpNode::buildSimilarBinaryOpNode(NodeID alt_arg1, NodeID alt_arg2, DataTre
       return alt_datatree.AddEqualEqual(alt_arg1, alt_arg2);
     case oDifferent:
       return alt_datatree.AddDifferent(alt_arg1, alt_arg2);
+    case oExpectation:
+      return alt_datatree.AddExpectation((int)(alt_arg1->eval(map<int, double>())), alt_arg2);
     }
   // Suppress GCC warning
   exit(EXIT_FAILURE);
@@ -2615,6 +2647,59 @@ BinaryOpNode::substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *
   return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
 }
 
+NodeID
+BinaryOpNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const
+{
+  switch(op_code)
+    {
+    case oExpectation:
+      {
+        int period = (int)(arg1->eval(map<int, double>()));
+        subst_table_t::iterator it = subst_table.find(const_cast<BinaryOpNode *>(this));
+
+        //IF should evaluate to true when substituting Exp operators out of equations in second pass
+        if (it != subst_table.end())
+          return const_cast<VariableNode *>(it->second);
+
+        //Arriving here, we need to create an auxiliary variable for this Expectation Operator:
+        int symb_id = datatree.symbol_table.addExpectationAuxiliaryVar(arg1->idx, arg2->idx); //AUXE_arg1.idx_arg2.idx
+        NodeID newAuxE = datatree.AddVariable(symb_id, 0);
+        assert(dynamic_cast<VariableNode *>(newAuxE) != NULL);
+
+        if (partial_information_model && period==0)
+          {
+            //Ensure x is a single variable as opposed to an expression
+            if (dynamic_cast<VariableNode *>(arg2) == NULL)
+              {
+                cerr << "In Partial Information models, EXPECTATION(0)(X) can only be used when X is a single variable." << endl;
+                exit(EXIT_FAILURE);
+              }
+          }
+        else
+          {
+            //take care of any nested expectation operators by calling arg2->substituteExpectation(.), then decreaseLeadsLags for this oExp operator
+            //arg2(lag-period) (holds entire subtree of arg2(lag-period)
+            NodeID substexpr = (arg2->substituteExpectation(subst_table, neweqs, partial_information_model))->decreaseLeadsLags(period);
+            assert(substexpr != NULL);
+
+            neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(newAuxE, substexpr))); //AUXE_arg1.idx_arg2.idx = arg2(lag-period)
+
+            newAuxE = newAuxE->decreaseLeadsLags(-1*period);
+            assert(dynamic_cast<VariableNode *>(newAuxE) != NULL);
+          }
+
+        subst_table[this] = dynamic_cast<VariableNode *>(newAuxE);
+        return newAuxE;
+      }
+    default:
+      {
+        NodeID arg1subst = arg1->substituteExpectation(subst_table, neweqs, partial_information_model);
+        NodeID arg2subst = arg2->substituteExpectation(subst_table, neweqs, partial_information_model);
+        return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
+      }
+    }
+}
+
 TrinaryOpNode::TrinaryOpNode(DataTree &datatree_arg, const NodeID arg1_arg,
                              TrinaryOpcode op_code_arg, const NodeID arg2_arg, const NodeID arg3_arg) :
     ExprNode(datatree_arg),
@@ -3018,6 +3103,15 @@ TrinaryOpNode::substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode
   return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
 }
 
+NodeID
+TrinaryOpNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const
+{
+  NodeID arg1subst = arg1->substituteExpectation(subst_table, neweqs, partial_information_model);
+  NodeID arg2subst = arg2->substituteExpectation(subst_table, neweqs, partial_information_model);
+  NodeID arg3subst = arg3->substituteExpectation(subst_table, neweqs, partial_information_model);
+  return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
+}
+
 UnknownFunctionNode::UnknownFunctionNode(DataTree &datatree_arg,
     int symb_id_arg,
     const vector<NodeID> &arguments_arg) :
@@ -3203,3 +3297,10 @@ UnknownFunctionNode::substituteExoLag(subst_table_t &subst_table, vector<BinaryO
   cerr << "UnknownFunctionNode::substituteExoLag: not implemented!" << endl;
   exit(EXIT_FAILURE);
 }
+
+NodeID
+UnknownFunctionNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const
+{
+  cerr << "UnknownFunctionNode::substituteExpectation: not implemented!" << endl;
+  exit(EXIT_FAILURE);
+}
diff --git a/ExprNode.hh b/ExprNode.hh
index e6ae9c04e89c006a42f654c2ab7d040f26943541..2a7d11a82c8a3875f66bea372cc5a1b4e5eda91f 100644
--- a/ExprNode.hh
+++ b/ExprNode.hh
@@ -296,6 +296,13 @@ public:
     \param[out] neweqs Equations to be added to the model to match the creation of auxiliary variables.
   */
   virtual NodeID substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
+
+  //! Constructs a new expression where the expectation operator has been replaced by auxiliary variables
+  /*!
+    \param[in,out] subst_table Map used to store expressions that have already be substituted and their corresponding variable, in order to avoid creating two auxiliary variables for the same sub-expr.
+    \param[out] neweqs Equations to be added to the model to match the creation of auxiliary variables.
+  */
+  virtual NodeID substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const = 0;
 };
 
 //! Object used to compare two nodes (using their indexes)
@@ -317,6 +324,7 @@ private:
   virtual NodeID computeDerivative(int deriv_id);
 public:
   NumConstNode(DataTree &datatree_arg, int id_arg);
+  int get_id() const { return id; };
   virtual void prepareForDerivation();
   virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_type &temporary_terms) const;
   virtual void collectVariables(SymbolType type_arg, set<pair<int, int> > &result) const;
@@ -333,6 +341,7 @@ public:
   virtual NodeID substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual NodeID substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual NodeID substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual NodeID substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
 };
 
 //! Symbol or variable node
@@ -370,6 +379,7 @@ public:
   virtual NodeID substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual NodeID substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual NodeID substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual NodeID substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
 };
 
 //! Unary operator node
@@ -415,6 +425,7 @@ public:
   virtual NodeID substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual NodeID substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual NodeID substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual NodeID substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
 };
 
 //! Binary operator node
@@ -465,6 +476,7 @@ public:
   virtual NodeID substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual NodeID substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual NodeID substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual NodeID substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
 };
 
 //! Trinary operator node
@@ -509,6 +521,7 @@ public:
   virtual NodeID substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual NodeID substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual NodeID substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual NodeID substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
 };
 
 //! Unknown function node
@@ -545,6 +558,7 @@ public:
   virtual NodeID substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual NodeID substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual NodeID substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual NodeID substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
 };
 
 #endif
diff --git a/ModFile.cc b/ModFile.cc
index fb92b5233f9a3a80733a15e6a21c6c6967091fc8..d7f618e11cc01989b2afb8fc0a356f4f583d94b9 100644
--- a/ModFile.cc
+++ b/ModFile.cc
@@ -133,12 +133,15 @@ ModFile::checkPass()
 void
 ModFile::transformPass()
 {
-  // In stochastic models, create auxiliary vars for leads and lags greater than 2
   if (mod_file_struct.stoch_simul_present
       || mod_file_struct.estimation_present
       || mod_file_struct.osr_present
       || mod_file_struct.ramsey_policy_present)
     {
+      // In stochastic models, create auxiliary vars for Expecatation operator
+      dynamic_model.substituteExpectation(mod_file_struct.partial_information);
+
+      // In stochastic models, create auxiliary vars for leads and lags greater than 2
       dynamic_model.substituteEndoLeadGreaterThanTwo();
       dynamic_model.substituteExoLead();
       dynamic_model.substituteEndoLagGreaterThanTwo();
@@ -168,6 +171,9 @@ ModFile::transformPass()
 void
 ModFile::computingPass(bool no_tmp_terms)
 {
+
+  //  expressions_tree.replace_oExpectation_in_datatree();
+
   // Mod file may have no equation (for example in a standalone BVAR estimation)
   bool dynamic_model_needed = mod_file_struct.simul_present || mod_file_struct.check_present || mod_file_struct.stoch_simul_present
     || mod_file_struct.estimation_present|| mod_file_struct.osr_present
diff --git a/ModelTree.cc b/ModelTree.cc
index ac19cc4ca55fb9b556e9229136c128cdb47363e2..bf8ecb96b0dcd30930eadc9a39e434eac7b9a113 100644
--- a/ModelTree.cc
+++ b/ModelTree.cc
@@ -212,18 +212,43 @@ ModelTree::writeModelEquations(ostream &output, ExprNodeOutputType output_type)
   for (int eq = 0; eq < (int) equations.size(); eq++)
     {
       BinaryOpNode *eq_node = equations[eq];
-
       NodeID lhs = eq_node->get_arg1();
-      output << "lhs =";
-      lhs->writeOutput(output, output_type, temporary_terms);
-      output << ";" << endl;
-
       NodeID rhs = eq_node->get_arg2();
-      output << "rhs =";
-      rhs->writeOutput(output, output_type, temporary_terms);
-      output << ";" << endl;
 
-      output << "residual" << LEFT_ARRAY_SUBSCRIPT(output_type) << eq + ARRAY_SUBSCRIPT_OFFSET(output_type) << RIGHT_ARRAY_SUBSCRIPT(output_type) << "= lhs-rhs;" << endl;
+      // Test if the right hand side of the equation is empty.
+      double vrhs = 1.0;
+      try
+        {
+          vrhs = rhs->eval(eval_context_type());
+        }
+      catch(ExprNode::EvalException &e)
+        {
+        }
+
+      if (vrhs!=0)// The right hand side of the equation is not empty ==> residual=lhs-rhs;
+        {
+          output << "lhs =";
+          lhs->writeOutput(output, output_type, temporary_terms);
+          output << ";" << endl;
+
+          output << "rhs =";
+          rhs->writeOutput(output, output_type, temporary_terms);
+          output << ";" << endl;
+
+          output << "residual" << LEFT_ARRAY_SUBSCRIPT(output_type) 
+                 << eq + ARRAY_SUBSCRIPT_OFFSET(output_type) 
+                 << RIGHT_ARRAY_SUBSCRIPT(output_type) 
+                 << "= lhs-rhs;" << endl;
+        }
+      else// The right hand side of the equation is empty ==> residual=lhs;
+        {
+          output << "residual" << LEFT_ARRAY_SUBSCRIPT(output_type) 
+                 << eq + ARRAY_SUBSCRIPT_OFFSET(output_type) 
+                 << RIGHT_ARRAY_SUBSCRIPT(output_type) 
+                 << " = "; 
+          lhs->writeOutput(output, output_type, temporary_terms);
+          output << ";" << endl;      
+        }
     }
 }
 
diff --git a/NumericalInitialization.cc b/NumericalInitialization.cc
index 64a14b49adb852009e5d63289c32fc0a8d9a2f84..f90b2405ddd8fe9caee568173809cf157e70ffe8 100644
--- a/NumericalInitialization.cc
+++ b/NumericalInitialization.cc
@@ -142,6 +142,16 @@ EndValStatement::EndValStatement(const init_values_type &init_values_arg,
 }
 
 
+void
+EndValStatement::checkPass(ModFileStructure &mod_file_struct)
+{
+  if (mod_file_struct.shocks_present)
+    {
+      cerr << "ERROR: Putting a \"shocks\" block before an \"endval\" block is not permitted. Please swap the two blocks. This limitation will be removed in the next major release of Dynare." << endl;
+      exit(EXIT_FAILURE);
+    }
+}
+
 void
 EndValStatement::writeOutput(ostream &output, const string &basename) const
 {
diff --git a/NumericalInitialization.hh b/NumericalInitialization.hh
index 8379941870f993ab784cef4b8f0fc2876a42cc3d..f01e6436ab0a8918d5ca4d922385eadfa5428b03 100644
--- a/NumericalInitialization.hh
+++ b/NumericalInitialization.hh
@@ -79,6 +79,8 @@ class EndValStatement : public InitOrEndValStatement
 public:
   EndValStatement(const init_values_type &init_values_arg,
                   const SymbolTable &symbol_table_arg);
+  //! Workaround for trac ticket #35
+  virtual void checkPass(ModFileStructure &mod_file_struct);
   virtual void writeOutput(ostream &output, const string &basename) const;
 };
 
diff --git a/ParsingDriver.cc b/ParsingDriver.cc
index f5e4d3d347aaebb4c08f5adedbaed11538481a0f..45d3b2ea752ff93a81a3fa0159fe2ef91e6189b5 100644
--- a/ParsingDriver.cc
+++ b/ParsingDriver.cc
@@ -224,6 +224,8 @@ ParsingDriver::add_expression_variable(string *name)
 void
 ParsingDriver::periods(string *periods)
 {
+  warning("periods: this command is now deprecated and may be removed in a future version of Dynare. Please of the \"periods\" option of \"simul\" command instead.");
+
   int periods_val = atoi(periods->c_str());
   mod_file->addStatement(new PeriodsStatement(periods_val));
   delete periods;
@@ -419,13 +421,8 @@ ParsingDriver::end_shocks()
 void
 ParsingDriver::end_mshocks()
 {
-  mod_file->addStatement(new MShocksStatement(det_shocks, var_shocks, std_shocks,
-                                              covar_shocks, corr_shocks, mod_file->symbol_table));
+  mod_file->addStatement(new MShocksStatement(det_shocks, mod_file->symbol_table));
   det_shocks.clear();
-  var_shocks.clear();
-  std_shocks.clear();
-  covar_shocks.clear();
-  corr_shocks.clear();
 }
 
 void
@@ -572,6 +569,8 @@ ParsingDriver::add_value(string *p1)
 void
 ParsingDriver::do_sigma_e()
 {
+  warning("Sigma_e: this command is now deprecated and may be removed in a future version of Dynare. Please use the \"shocks\" command instead.");
+
   try
     {
       mod_file->addStatement(new SigmaeStatement(sigmae_matrix));
@@ -1270,6 +1269,14 @@ ParsingDriver::add_power(NodeID arg1, NodeID arg2)
   return data_tree->AddPower(arg1, arg2);
 }
 
+NodeID
+ParsingDriver::add_expectation(string *arg1, NodeID arg2)
+{
+  NodeID expectationNode = data_tree->AddExpectation(atoi(arg1->c_str()), arg2);
+  delete arg1;
+  return expectationNode;
+}
+
 NodeID
 ParsingDriver::add_exp(NodeID arg1)
 {
diff --git a/ParsingDriver.hh b/ParsingDriver.hh
index decd321f4755b3fb01cfe629f9e0c1597dbbe6ba..e2e35a91d99883fb5e2f248b08ffd549cf1b0495 100644
--- a/ParsingDriver.hh
+++ b/ParsingDriver.hh
@@ -392,6 +392,8 @@ public:
   NodeID add_different(NodeID arg1, NodeID arg2);
   //! Writes token "arg1^arg2" to model tree
   NodeID add_power(NodeID arg1,  NodeID arg2);
+  //! Writes token "E(arg1)(arg2)" to model tree
+  NodeID add_expectation(string *arg1,  NodeID arg2);
   //! Writes token "exp(arg1)" to model tree
   NodeID add_exp(NodeID arg1);
   //! Writes token "log(arg1)" to model tree
diff --git a/Shocks.cc b/Shocks.cc
index ca974a60392c2b947d85715f60357e8a40e5410a..79e19de7e8141e918ab9258ddbb4376a345b785f 100644
--- a/Shocks.cc
+++ b/Shocks.cc
@@ -27,17 +27,9 @@ using namespace std;
 
 AbstractShocksStatement::AbstractShocksStatement(bool mshocks_arg,
                                                  const det_shocks_type &det_shocks_arg,
-                                                 const var_and_std_shocks_type &var_shocks_arg,
-                                                 const var_and_std_shocks_type &std_shocks_arg,
-                                                 const covar_and_corr_shocks_type &covar_shocks_arg,
-                                                 const covar_and_corr_shocks_type &corr_shocks_arg,
                                                  const SymbolTable &symbol_table_arg) :
   mshocks(mshocks_arg),
   det_shocks(det_shocks_arg),
-  var_shocks(var_shocks_arg),
-  std_shocks(std_shocks_arg),
-  covar_shocks(covar_shocks_arg),
-  corr_shocks(corr_shocks_arg),
   symbol_table(symbol_table_arg)
 {
 }
@@ -83,7 +75,47 @@ AbstractShocksStatement::writeDetShocks(ostream &output) const
 }
 
 void
-AbstractShocksStatement::writeVarAndStdShocks(ostream &output) const
+AbstractShocksStatement::checkPass(ModFileStructure &mod_file_struct)
+{
+  mod_file_struct.shocks_present = true;
+}
+
+
+ShocksStatement::ShocksStatement(const det_shocks_type &det_shocks_arg,
+                                 const var_and_std_shocks_type &var_shocks_arg,
+                                 const var_and_std_shocks_type &std_shocks_arg,
+                                 const covar_and_corr_shocks_type &covar_shocks_arg,
+                                 const covar_and_corr_shocks_type &corr_shocks_arg,
+                                 const SymbolTable &symbol_table_arg) :
+  AbstractShocksStatement(false, det_shocks_arg, symbol_table_arg),
+  var_shocks(var_shocks_arg),
+  std_shocks(std_shocks_arg),
+  covar_shocks(covar_shocks_arg),
+  corr_shocks(corr_shocks_arg)
+{
+}
+
+void
+ShocksStatement::writeOutput(ostream &output, const string &basename) const
+{
+  output << "%" << endl
+         << "% SHOCKS instructions" << endl
+         << "%" << endl;
+
+  // Write instruction that initializes a shock
+  output << "make_ex_;" << endl;
+
+  writeDetShocks(output);
+  writeVarAndStdShocks(output);
+  writeCovarAndCorrShocks(output);
+  if (covar_shocks.size()+corr_shocks.size() > 0)
+    output << "M_.sigma_e_is_diagonal = 0;" << endl;
+  else
+    output << "M_.sigma_e_is_diagonal = 1;" << endl;
+}
+
+void
+ShocksStatement::writeVarAndStdShocks(ostream &output) const
 {
   var_and_std_shocks_type::const_iterator it;
 
@@ -107,7 +139,7 @@ AbstractShocksStatement::writeVarAndStdShocks(ostream &output) const
 }
 
 void
-AbstractShocksStatement::writeCovarAndCorrShocks(ostream &output) const
+ShocksStatement::writeCovarAndCorrShocks(ostream &output) const
 {
   covar_and_corr_shocks_type::const_iterator it;
 
@@ -135,41 +167,9 @@ AbstractShocksStatement::writeCovarAndCorrShocks(ostream &output) const
     }
 }
 
-
-ShocksStatement::ShocksStatement(const det_shocks_type &det_shocks_arg,
-                                 const var_and_std_shocks_type &var_shocks_arg,
-                                 const var_and_std_shocks_type &std_shocks_arg,
-                                 const covar_and_corr_shocks_type &covar_shocks_arg,
-                                 const covar_and_corr_shocks_type &corr_shocks_arg,
-                                 const SymbolTable &symbol_table_arg) :
-  AbstractShocksStatement(false, det_shocks_arg, var_shocks_arg, std_shocks_arg,
-                          covar_shocks_arg, corr_shocks_arg, symbol_table_arg)
-{
-}
-
-void
-ShocksStatement::writeOutput(ostream &output, const string &basename) const
-{
-  output << "%" << endl
-         << "% SHOCKS instructions" << endl
-         << "%" << endl;
-
-  // Write instruction that initializes a shock
-  output << "make_ex_;" << endl;
-
-  writeDetShocks(output);
-  writeVarAndStdShocks(output);
-  writeCovarAndCorrShocks(output);
-}
-
 MShocksStatement::MShocksStatement(const det_shocks_type &det_shocks_arg,
-                                   const var_and_std_shocks_type &var_shocks_arg,
-                                   const var_and_std_shocks_type &std_shocks_arg,
-                                   const covar_and_corr_shocks_type &covar_shocks_arg,
-                                   const covar_and_corr_shocks_type &corr_shocks_arg,
                                    const SymbolTable &symbol_table_arg) :
-  AbstractShocksStatement(true, det_shocks_arg, var_shocks_arg, std_shocks_arg,
-                          covar_shocks_arg, corr_shocks_arg, symbol_table_arg)
+  AbstractShocksStatement(true, det_shocks_arg, symbol_table_arg)
 {
 }
 
@@ -177,15 +177,13 @@ void
 MShocksStatement::writeOutput(ostream &output, const string &basename) const
 {
   output << "%" << endl
-         << "% SHOCKS instructions" << endl
+         << "% MSHOCKS instructions" << endl
          << "%" << endl;
 
   // Write instruction that initializes a shock
   output << "make_ex_;" << endl;
 
   writeDetShocks(output);
-  writeVarAndStdShocks(output);
-  writeCovarAndCorrShocks(output);
 }
 
 ConditionalForecastPathsStatement::ConditionalForecastPathsStatement(const AbstractShocksStatement::det_shocks_type &paths_arg, const SymbolTable &symbol_table_arg) :
diff --git a/Shocks.hh b/Shocks.hh
index b396ae2a49f86054366a0e173e21b4d949d26900..a06e06c2d2f41cba30be33fc507e740131c45cfc 100644
--- a/Shocks.hh
+++ b/Shocks.hh
@@ -40,30 +40,30 @@ public:
     NodeID value;
   };
   typedef map<string, vector<DetShockElement> > det_shocks_type;
-  typedef map<string, NodeID> var_and_std_shocks_type;
-  typedef map<pair<string, string>, NodeID> covar_and_corr_shocks_type;
+  //! Workaround for trac ticket #35
+  virtual void checkPass(ModFileStructure &mod_file_struct);
 protected:
   //! Is this statement a "mshocks" statement ? (instead of a "shocks" statement)
   const bool mshocks;
   const det_shocks_type det_shocks;
-  const var_and_std_shocks_type var_shocks, std_shocks;
-  const covar_and_corr_shocks_type covar_shocks, corr_shocks;
   const SymbolTable &symbol_table;
   void writeDetShocks(ostream &output) const;
-  void writeVarAndStdShocks(ostream &output) const;
-  void writeCovarAndCorrShocks(ostream &output) const;
 
   AbstractShocksStatement(bool mshocks_arg,
                           const det_shocks_type &det_shocks_arg,
-                          const var_and_std_shocks_type &var_shocks_arg,
-                          const var_and_std_shocks_type &std_shocks_arg,
-                          const covar_and_corr_shocks_type &covar_shocks_arg,
-                          const covar_and_corr_shocks_type &corr_shocks_arg,
                           const SymbolTable &symbol_table_arg);
 };
 
 class ShocksStatement : public AbstractShocksStatement
 {
+public:
+  typedef map<string, NodeID> var_and_std_shocks_type;
+  typedef map<pair<string, string>, NodeID> covar_and_corr_shocks_type;
+private:
+  const var_and_std_shocks_type var_shocks, std_shocks;
+  const covar_and_corr_shocks_type covar_shocks, corr_shocks;
+  void writeVarAndStdShocks(ostream &output) const;
+  void writeCovarAndCorrShocks(ostream &output) const;
 public:
   ShocksStatement(const det_shocks_type &det_shocks_arg,
                   const var_and_std_shocks_type &var_shocks_arg,
@@ -78,10 +78,6 @@ class MShocksStatement : public AbstractShocksStatement
 {
 public:
   MShocksStatement(const det_shocks_type &det_shocks_arg,
-                   const var_and_std_shocks_type &var_shocks_arg,
-                   const var_and_std_shocks_type &std_shocks_arg,
-                   const covar_and_corr_shocks_type &covar_shocks_arg,
-                   const covar_and_corr_shocks_type &corr_shocks_arg,
                    const SymbolTable &symbol_table_arg);
   virtual void writeOutput(ostream &output, const string &basename) const;
 };
diff --git a/Statement.cc b/Statement.cc
index 917f645b4ce6d3ea5a3401b06826d8f93be2c7b5..97d4ad64a55beccfd07ae5ac4bcc289ffaa88bd5 100644
--- a/Statement.cc
+++ b/Statement.cc
@@ -30,7 +30,8 @@ ModFileStructure::ModFileStructure() :
   bvar_density_present(false),
   bvar_forecast_present(false),
   identification_present(false),
-  partial_information(false)
+  partial_information(false),
+  shocks_present(false)
 {
 }
 
diff --git a/Statement.hh b/Statement.hh
index 43595fb344a65ab49d87c55a0d697f1052c55f44..cfeb071dfd252a68185d50d16ca687b3f0df7611 100644
--- a/Statement.hh
+++ b/Statement.hh
@@ -60,6 +60,9 @@ public:
   bool identification_present;
   //! Whether the option partial_information is given to stoch_simul/estimation/osr/ramsey_policy
   bool partial_information;
+  //! Whether a shocks or mshocks block is present
+  /*! Used for the workaround for trac ticket #35 */
+  bool shocks_present;
 };
 
 class Statement
diff --git a/StaticModel.cc b/StaticModel.cc
index d848b7390b58234205a1cd5747a25235acb79a41..6a55954490bd2cccb50eb81de98f8c02764a67ab 100644
--- a/StaticModel.cc
+++ b/StaticModel.cc
@@ -92,9 +92,7 @@ StaticModel::writeStaticMFile(ostream &output, const string &func_name) const
       int symb_id = it->first.second;
       NodeID d1 = it->second;
 
-      ostringstream g1;
-      g1 << "  g1(" << eq+1 << "," << symbol_table.getTypeSpecificID(symb_id)+1 << ")";
-      output << g1.str() << "=" << g1.str() << "+";
+      output << "  g1(" << eq+1 << "," << symbol_table.getTypeSpecificID(symb_id)+1 << ")=";
       d1->writeOutput(output, oMatlabStaticModel, temporary_terms);
       output << ";" << endl;
     }
diff --git a/SymbolTable.cc b/SymbolTable.cc
index 8fc5195dcecac172abc9f54b956c8abbaf1423e8..16e68a2bc052dbfd9624619bebc953598c735b29 100644
--- a/SymbolTable.cc
+++ b/SymbolTable.cc
@@ -302,3 +302,27 @@ SymbolTable::addExoLagAuxiliaryVar(int orig_symb_id, int orig_lag) throw (Frozen
 {
   return addLagAuxiliaryVarInternal(false, orig_symb_id, orig_lag);
 }
+
+int
+SymbolTable::addExpectationAuxiliaryVar(int arg1, int arg2) throw (FrozenException)
+{
+  ostringstream varname;
+  varname << "AUXE_" << arg1 << "_" << arg2;
+  int symb_id;
+  try
+    {
+      symb_id = addSymbol(varname.str(), eEndogenous);
+    }
+  catch(AlreadyDeclaredException &e)
+    {
+      cerr << "ERROR: you should rename your variable called " << varname.str() << ", this name is internally used by Dynare" << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  AuxVarInfo avi;
+  avi.symb_id = symb_id;
+  avi.type = avExpectation;
+  aux_vars.push_back(avi);
+
+  return symb_id;
+}
diff --git a/SymbolTable.hh b/SymbolTable.hh
index 89b11cac25b855e6e383082550d6bd9f433df151..461ae20c03299f14cb8d499edb209a2b64767100 100644
--- a/SymbolTable.hh
+++ b/SymbolTable.hh
@@ -32,10 +32,11 @@ using namespace std;
 //! Types of auxiliary variables
 enum aux_var_t
   {
-    avEndoLead = 0, //!< Substitute for endo leads >= 2
-    avEndoLag = 1,  //!< Substitute for endo lags >= 2
-    avExoLead = 2,  //!< Substitute for exo leads >= 2
-    avExoLag = 3,   //!< Substitute for exo lags >= 2
+    avEndoLead = 0,   //!< Substitute for endo leads >= 2
+    avEndoLag = 1,    //!< Substitute for endo lags >= 2
+    avExoLead = 2,    //!< Substitute for exo leads >= 2
+    avExoLag = 3,     //!< Substitute for exo lags >= 2
+    avExpectation = 4 //!< Substitute for Expectation Operator
   };
 
 //! Information on some auxiliary variables
@@ -171,6 +172,11 @@ public:
     \param[in] orig_lag lag value such that this new variable will be equivalent to orig_symb_id(orig_lag)
     \return the symbol ID of the new symbol */
   int addExoLagAuxiliaryVar(int orig_symb_id, int orig_lag) throw (FrozenException);
+  //! Adds an auxiliary variable for the expectations operator
+  /*!
+    \param[in] indeces Used to construct the variable name
+    \return the symbol ID of the new symbol */
+  int addExpectationAuxiliaryVar(int arg1, int arg2) throw (FrozenException);
   //! Tests if symbol already exists
   inline bool exists(const string &name) const;
   //! Get symbol name (by ID)