diff --git a/parser.src/ComputingTasks.cc b/parser.src/ComputingTasks.cc
index 8e375181ea04e2df75495c4b4b5fd794fc3df809..830fffce90ccc7a401ef8457e790fcf9f732c2dc 100644
--- a/parser.src/ComputingTasks.cc
+++ b/parser.src/ComputingTasks.cc
@@ -7,569 +7,626 @@
 #include <sstream>
 
 using namespace std;
-//------------------------------------------------------------------------------
+
 #include "ComputingTasks.hh"
 #include "Interface.hh"
-//------------------------------------------------------------------------------
-//ostringstream	ComputingTasks::output;
-//------------------------------------------------------------------------------
-ComputingTasks::ComputingTasks(const SymbolTable &symbol_table_arg) : symbol_table(symbol_table_arg)
+#include "Statement.hh"
+
+SteadyStatement::SteadyStatement(const OptionsList &options_list_arg) :
+  options_list(options_list_arg)
+{
+}
+
+void
+SteadyStatement::writeOutput(ostream &output) const
+{
+  options_list.writeOutput(output);
+  output << "steady;\n";
+}
+
+CheckStatement::CheckStatement(const OptionsList &options_list_arg) :
+  options_list(options_list_arg)
+{
+}
+
+void
+CheckStatement::writeOutput(ostream &output) const
 {
-  // Empty
+  options_list.writeOutput(output);
+  output << "check;\n";
 }
 
-//------------------------------------------------------------------------------
-ComputingTasks::~ComputingTasks()
+SimulStatement::SimulStatement(const OptionsList &options_list_arg) :
+  options_list(options_list_arg)
 {
-  // Empty
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::setOutput(ostringstream* iOutput)
+void
+SimulStatement::writeOutput(ostream &output) const
 {
-  output = iOutput;
+  options_list.writeOutput(output);
+  output << "simul(oo_.dr);\n";
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::set(void)
+StochSimulStatement::StochSimulStatement(const TmpSymbolTable &tmp_symbol_table_arg,
+                                         const OptionsList &options_list_arg) :
+  tmp_symbol_table(tmp_symbol_table_arg),
+  options_list(options_list_arg)
 {
-  // Empty
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::setSteady(void)
+void
+StochSimulStatement::writeOutput(ostream &output) const
 {
-  *output << "steady;\n";
+  options_list.writeOutput(output);
+  tmp_symbol_table.writeOutput("var_list_", output);
+  output << "stoch_simul(var_list_);\n";
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::setCheck(void)
+EstimationStatement::EstimationStatement(const TmpSymbolTable &tmp_symbol_table_arg,
+                                         const OptionsList &options_list_arg) :
+  tmp_symbol_table(tmp_symbol_table_arg),
+  options_list(options_list_arg)
 {
-  *output << "check;\n";
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::setSimul(void)
+void
+EstimationStatement::writeOutput(ostream &output) const
 {
-  *output << "simul(oo_.dr);\n";
+  options_list.writeOutput(output);
+  tmp_symbol_table.writeOutput("var_list_", output);
+  output << "dynare_estimation(var_list_);\n";
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::setStochSimul(string tmp1)
+PriorAnalysisStatement::PriorAnalysisStatement(const TmpSymbolTable &tmp_symbol_table_arg,
+                                               const OptionsList &options_list_arg) :
+  tmp_symbol_table(tmp_symbol_table_arg),
+  options_list(options_list_arg)
+{
+}
+
+void
+PriorAnalysisStatement::writeOutput(ostream &output) const
+{
+  options_list.writeOutput(output);
+  tmp_symbol_table.writeOutput("var_list_", output);
+  output << "prior_analysis(var_list_);\n";
+}
+
+PosteriorAnalysisStatement::PosteriorAnalysisStatement(const TmpSymbolTable &tmp_symbol_table_arg,
+                                                       const OptionsList &options_list_arg) :
+  tmp_symbol_table(tmp_symbol_table_arg),
+  options_list(options_list_arg)
+{
+}
+
+void
+PosteriorAnalysisStatement::writeOutput(ostream &output) const
+{
+  options_list.writeOutput(output);
+  tmp_symbol_table.writeOutput("var_list_", output);
+  output << "posterior_analysis(var_list_);\n";
+}
 
+RplotStatement::RplotStatement(const TmpSymbolTable &tmp_symbol_table_arg,
+                               const OptionsList &options_list_arg) :
+  tmp_symbol_table(tmp_symbol_table_arg),
+  options_list(options_list_arg)
 {
-  *output << tmp1;
-  *output << "stoch_simul(var_list_);\n";
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::setOption(string name, string value)
+void
+RplotStatement::writeOutput(ostream &output) const
 {
-  *output << "options_." << name << " = " << value << ";\n";
+  options_list.writeOutput(output);
+  tmp_symbol_table.writeOutput("var_list_", output);
+  output << "rplot(var_list_);\n";
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::setOption(string name, string value1, string value2)
+UnitRootVarsStatement::UnitRootVarsStatement(const TmpSymbolTable &tmp_symbol_table_arg) :
+  tmp_symbol_table(tmp_symbol_table_arg)
 {
-  *output << "options_." << name << " = [" << value1 << "; " << value2 << "];\n";
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::runEstimation(string tmp1)
+void
+UnitRootVarsStatement::writeOutput(ostream &output) const
 {
-  *output << tmp1;
-  *output << "dynare_estimation(var_list_);\n";
+  tmp_symbol_table.writeOutput("options_.unit_root_vars", output);
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::runPriorAnalysis(string tmp1)
+PeriodsStatement::PeriodsStatement(int periods_arg) : periods(periods_arg)
 {
-  *output << tmp1;
-  *output << "prior_analysis(var_list_);\n";
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::runPosteriorAnalysis(string tmp1)
+void
+PeriodsStatement::writeOutput(ostream &output) const
 {
-  *output << tmp1;
-  *output << "posterior_analysis(var_list_);\n";
+  output << "options_.periods = " << periods << ";" << endl;
+  output << "options_.simul = 1;" << endl;
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::runRplot(string tmp1)
+DsampleStatement::DsampleStatement(int val1_arg) : val1(val1_arg), val2(-1)
 {
-  *output << tmp1;
-  *output << "rplot(var_list_);\n";
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::setEstimationInit(void)
+DsampleStatement::DsampleStatement(int val1_arg, int val2_arg) : val1(val1_arg), val2(val2_arg)
 {
-  *output << "global estim_params_\n";
-  *output << "var_list_ = [];\n";
-  *output << "estim_params_.var_exo = [];\n";
-  *output << "estim_params_.var_endo = [];\n";
-  *output << "estim_params_.corrx = [];\n";
-  *output << "estim_params_.corrn = [];\n";
-  *output << "estim_params_.param_names = [];\n";
-  *output << "estim_params_.user_param_names = [];\n";
-  *output << "estim_params_.param_vals = [];\n";
-  *output << "M_.H = 0;\n";
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::setOptimOptions(string str1, string str2, int task)
+void
+DsampleStatement::writeOutput(ostream &output) const
 {
-  string optim_string;
-  int start;
-  switch(task)
+  if (val2 < 0)
+    output << "options_.dsample = " << val1 << ";" << endl;
+  else
+    output << "options_.dsample = [" << val1 << "; " << val2 << "];" << endl;
+}
+
+VarobsStatement::VarobsStatement(const TmpSymbolTable &tmp_symbol_table_arg) :
+  tmp_symbol_table(tmp_symbol_table_arg)
+{
+}
+
+void
+VarobsStatement::writeOutput(ostream &output) const
+{
+  tmp_symbol_table.writeOutput("options_.varobs", output);
+}
+
+EstimatedParamsStatement::EstimatedParamsStatement(const vector<EstimationParams> &estim_params_list_arg,
+                                                   const SymbolTable &symbol_table_arg) :
+  estim_params_list(estim_params_list_arg),
+  symbol_table(symbol_table_arg)
+{
+}
+
+void
+EstimatedParamsStatement::writeOutput(ostream &output) const
+{
+  output << "global estim_params_\n";
+  output << "var_list_ = [];\n";
+  output << "estim_params_.var_exo = [];\n";
+  output << "estim_params_.var_endo = [];\n";
+  output << "estim_params_.corrx = [];\n";
+  output << "estim_params_.corrn = [];\n";
+  output << "estim_params_.param_names = [];\n";
+  output << "estim_params_.user_param_names = [];\n";
+  output << "estim_params_.param_vals = [];\n";
+  output << "M_.H = 0;\n";
+
+  vector<EstimationParams>::const_iterator it;
+
+  for(it = estim_params_list.begin(); it != estim_params_list.end(); it++)
     {
-    case 1:
-      optim_string = "options_.optim_opt = '";
-      start = 0;
-      return;
-    case 2:
-      if (start > 0)
-        {
-          optim_string += ",";
-        }
-      else
+      if (symbol_table.isReferenced(it->name) == eNotReferenced
+          && it->name != "dsge_prior_weight")
+        return;
+      switch(it->type)
         {
-          start = 1;
+        case 1:
+          if (symbol_table.getType(it->name) == eExogenous)
+            output << "estim_params_.var_exo = [estim_params_.var_exo; ";
+          else if (symbol_table.getType(it->name) == eEndogenous)
+            output << "estim_params_.var_endo = [estim_params_.var_endo; ";
+          output << symbol_table.getID(it->name)+1;
+          break;
+        case 2:
+          output << "estim_params_.param_vals = [estim_params_.param_vals; ";
+          output << symbol_table.getID(it->name)+1;
+          break;
+        case 3:
+          if (symbol_table.getType(it->name) == eExogenous)
+            output << "estim_params_.corrx = [estim_params_.corrx; ";
+          else if (symbol_table.getType(it->name) == eEndogenous)
+            output << "estim_params_.corrn = [estim_params_.corrn; ";
+          output << symbol_table.getID(it->name)+1;
+          output << " " << symbol_table.getID(it->name2)+1;
+          break;
         }
-      optim_string += "''";
-      optim_string += str1;
-      optim_string += "'',";
-      if (str2[0] >= 'A' && str2[0] <= 'z')
-        {
-          optim_string += "''";
-          optim_string += str2;
-          optim_string += "''";
-        }
-      else
-        {
-          optim_string += str2;
-        }
-      return;
-    case 3:
-      optim_string += "';\n";
-      *output << optim_string;
+      output << " " << it->init_val << " " <<  it->low_bound
+             << " " << it->up_bound << " " <<  it->prior
+             << " " << it->mean << " " <<  it->std
+             << " " << it->p3 << " " <<  it->p4  << " " <<  it->jscale << "];\n";
     }
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::setEstimatedElements(void)
+EstimatedParamsInitStatement::EstimatedParamsInitStatement(const vector<EstimationParams> &estim_params_list_arg,
+                                                           const SymbolTable &symbol_table_arg) :
+  estim_params_list(estim_params_list_arg),
+  symbol_table(symbol_table_arg)
 {
-  if (!symbol_table.Exist(EstimParams->name))
-    {
-      string msg = "Unknown symbol: "+EstimParams->name;
-      error(msg.c_str());
-    }
-  if (symbol_table.isReferenced(EstimParams->name) == eNotReferenced & EstimParams->name != "dsge_prior_weight")
-    {
-      return;
-    }
-  if ((EstimParams->init_val).size() == 0)
-    {
-      EstimParams->init_val = EstimParams->mean;
-    }
-  switch(EstimParams->type)
+}
+
+void
+EstimatedParamsInitStatement::writeOutput(ostream &output) const
+{
+  vector<EstimationParams>::const_iterator it;
+
+  for(it = estim_params_list.begin(); it != estim_params_list.end(); it++)
     {
-    case 1:
-      if (symbol_table.getType(EstimParams->name) == eExogenous)
-        {
-          *output << "estim_params_.var_exo = [estim_params_.var_exo; ";
-        }
-      else if (symbol_table.getType(EstimParams->name) == eEndogenous)
-        {
-          *output << "estim_params_.var_endo = [estim_params_.var_endo; ";
-        }
-      *output << symbol_table.getID(EstimParams->name)+1;
-      break;
-    case 2:
-      *output << "estim_params_.param_vals = [estim_params_.param_vals; ";
-      *output << symbol_table.getID(EstimParams->name)+1;
-      break;
-    case 3:
-      if (symbol_table.getType(EstimParams->name) == eExogenous)
+      if (symbol_table.isReferenced(it->name) == eNotReferenced)
+        return;
+      if (it->type < 3)
         {
-          *output << "estim_params_.corrx = [estim_params_.corrx; ";
+          if (symbol_table.getType(it->name) == eExogenous)
+            {
+              output << "tmp1 = find(estim_params_.var_exo(:,1)==" << symbol_table.getID(it->name)+1 << ");\n";
+              output << "estim_params_.var_exo(tmp1,2) = " << it->init_val << ";\n";
+            }
+          else if (symbol_table.getType(it->name) == eEndogenous)
+            {
+              output << "tmp1 = find(estim_params_.var_endo(:,1)==" << symbol_table.getID(it->name)+1 << ");\n";
+              output << "estim_params_.var_endo(tmp1,2) = " << it->init_val << ";\n";
+            }
+          else if (symbol_table.getType(it->name) == eParameter)
+            {
+              output << "tmp1 = find(estim_params_.param_vals(:,1)==" << symbol_table.getID(it->name)+1 << ");\n";
+              output << "estim_params_.param_vals(tmp1,2) = " << it->init_val << ";\n";
+            }
         }
-      else if (symbol_table.getType(EstimParams->name) == eEndogenous)
+      else
         {
-          *output << "estim_params_.corrn = [estim_params_.corrn; ";
+          if (symbol_table.getType(it->name) == eExogenous)
+            {
+              output << "tmp1 = find((estim_params_.corrx(:,1)==" << symbol_table.getID(it->name)+1 << ")) & (estim_params_.corrx(:,2)==" << symbol_table.getID(it->name2)+1 << ");\n";
+              output << "estim_params_.corrx(tmp1,3) = " << it->init_val << ";\n";
+            }
+          else if (symbol_table.getType(it->name) == eEndogenous)
+            {
+              output << "tmp1 = find((estim_params_.corrn(:,1)==" << symbol_table.getID(it->name)+1 << ")) & (estim_params_.corrn(:,2)==" << symbol_table.getID(it->name2)+1 << ";\n";
+              output << "estim_params_.corrx(tmp1,3) = " << it->init_val << ";\n";
+            }
         }
-      *output << symbol_table.getID(EstimParams->name)+1;
-      *output << " " << symbol_table.getID(EstimParams->name2)+1;
-      break;
     }
-  *output << " " << EstimParams->init_val << " " <<  EstimParams->low_bound << " " <<
-    EstimParams->up_bound << " " <<  EstimParams->prior << " ";
-  *output <<  EstimParams->mean << " " <<  EstimParams->std << " " <<
-    EstimParams->p3 << " " <<  EstimParams->p4  << " " <<  EstimParams->jscale << "];\n";
-  EstimParams->clear();
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::setEstimatedInitElements(void)
+EstimatedParamsBoundsStatement::EstimatedParamsBoundsStatement(const vector<EstimationParams> &estim_params_list_arg,
+                                                               const SymbolTable &symbol_table_arg) :
+  estim_params_list(estim_params_list_arg),
+  symbol_table(symbol_table_arg)
 {
-  if (!symbol_table.Exist(EstimParams->name))
-    {
-      string msg = "Unknown symbol: "+EstimParams->name;
-      error(msg.c_str());
-    }
-  if (symbol_table.isReferenced(EstimParams->name) == eNotReferenced)
-    {
-      return;
-    }
-  if ((EstimParams->init_val).size() == 0)
-    {
-      EstimParams->init_val = EstimParams->mean;
-    }
-  if (EstimParams->type < 3)
+}
+
+void
+EstimatedParamsBoundsStatement::writeOutput(ostream &output) const
+{
+  vector<EstimationParams>::const_iterator it;
+
+  for(it = estim_params_list.begin(); it != estim_params_list.end(); it++)
     {
-      if (symbol_table.getType(EstimParams->name) == eExogenous)
+      if (symbol_table.isReferenced(it->name) == eNotReferenced)
+        return;
+      if (it->type < 3)
         {
-          *output << "tmp1 = find(estim_params_.var_exo(:,1)==" << symbol_table.getID(EstimParams->name)+1 << ");\n";
-          *output << "estim_params_.var_exo(tmp1,2) = " << EstimParams->init_val << ";\n";
+          if (symbol_table.getType(it->name) == eExogenous)
+            {
+              output << "tmp1 = find(estim_params_.var_exo(:,1)==" << symbol_table.getID(it->name)+1 << ");\n";
+              output << "estim_params_.var_exo(tmp1,3) = " << it->low_bound << ";\n";
+              output << "estim_params_.var_exo(tmp1,4) = " << it->up_bound << ";\n";
+            }
+          else if (symbol_table.getType(it->name) == eEndogenous)
+            {
+              output << "tmp1 = find(estim_params_.var_endo(:,1)==" << symbol_table.getID(it->name)+1 << ");\n";
+              output << "estim_params_.var_endo(tmp1,3) = " << it->low_bound << ";\n";
+              output << "estim_params_.var_endo(tmp1,4) = " << it->up_bound << ";\n";
+            }
+          else if (symbol_table.getType(it->name) == eParameter)
+            {
+              output << "tmp1 = find(estim_params_.param_vals(:,1)==" << symbol_table.getID(it->name)+1 << ");\n";
+              output << "estim_params_.param_vals(tmp1,3) = " << it->low_bound << ";\n";
+              output << "estim_params_.param_vals(tmp1,4) = " << it->up_bound << ";\n";
+            }
         }
-      else if (symbol_table.getType(EstimParams->name) == eEndogenous)
-        {
-          *output << "tmp1 = find(estim_params_.var_endo(:,1)==" << symbol_table.getID(EstimParams->name)+1 << ");\n";
-          *output << "estim_params_.var_endo(tmp1,2) = " << EstimParams->init_val << ";\n";
-        }
-      else if (symbol_table.getType(EstimParams->name) == eParameter)
+      else
         {
-          *output << "tmp1 = find(estim_params_.param_vals(:,1)==" << symbol_table.getID(EstimParams->name)+1 << ");\n";
-          *output << "estim_params_.param_vals(tmp1,2) = " << EstimParams->init_val << ";\n";
+          if (symbol_table.getType(it->name) == eExogenous)
+            {
+              output << "tmp1 = find((estim_params_.corrx(:,1)==" << symbol_table.getID(it->name)+1 << ")) & (estim_params_.corrx(:,2)==" << symbol_table.getID(it->name2)+1 << ");\n";
+              output << "estim_params_.corrx(tmp1,4) = " << it->low_bound << ";\n";
+              output << "estim_params_.corrx(tmp1,5) = " << it->up_bound << ";\n";
+            }
+          else if (symbol_table.getType(it->name) == eEndogenous)
+            {
+              output << "tmp1 = find((estim_params_.corrn(:,1)==" << symbol_table.getID(it->name)+1 << ")) & (estim_params_.corrn(:,2)==" << symbol_table.getID(it->name2)+1 << ";\n";
+              output << "estim_params_.corrx(tmp1,4) = " << it->low_bound << ";\n";
+              output << "estim_params_.corrx(tmp1,5) = " << it->up_bound << ";\n";
+            }
         }
     }
-  else
+}
+
+ObservationTrendsStatement::ObservationTrendsStatement(const trend_elements_type &trend_elements_arg,
+                                                       const SymbolTable &symbol_table_arg) :
+  trend_elements(trend_elements_arg),
+  symbol_table(symbol_table_arg)
+{
+}
+
+void
+ObservationTrendsStatement::writeOutput(ostream &output) const
+{
+  output << "options_.trend_coeff_ = {};" << endl;
+
+  trend_elements_type::const_iterator it;
+
+  for(it = trend_elements.begin(); it != trend_elements.end(); it++)
     {
-      if (symbol_table.getType(EstimParams->name) == eExogenous)
+      Type type = symbol_table.getType(it->first);
+      if (type == eEndogenous)
         {
-          *output << "tmp1 = find((estim_params_.corrx(:,1)==" << symbol_table.getID(EstimParams->name)+1 << ")) & (estim_params_.corrx(:,2)==" << symbol_table.getID(EstimParams->name2)+1 << ");\n";
-          *output << "estim_params_.corrx(tmp1,3) = " << EstimParams->init_val << ";\n";
-        }
-      else if (symbol_table.getType(EstimParams->name) == eEndogenous)
-        {
-          *output << "tmp1 = find((estim_params_.corrn(:,1)==" << symbol_table.getID(EstimParams->name)+1 << ")) & (estim_params_.corrn(:,2)==" << symbol_table.getID(EstimParams->name2)+1 << ";\n";
-          *output << "estim_params_.corrx(tmp1,3) = " << EstimParams->init_val << ";\n";
+          output << "tmp1 = strmatch('" << it->first << "',options_.varobs,'exact');\n";
+          output << "options_.trend_coeffs{tmp1} = '" << it->second << "';\n";
         }
+      else
+        cout << "Error : Non-variable symbol used in TREND_COEFF: " << it->first << endl;
     }
-  EstimParams->clear();
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::setEstimatedBoundsElements(void)
+CalibVarStatement::CalibVarStatement(const calib_var_type &calib_var_arg,
+                                     const calib_covar_type &calib_covar_arg,
+                                     const calib_ac_type &calib_ac_arg,
+                                     const SymbolTable &symbol_table_arg) :
+  calib_var(calib_var_arg),
+  calib_covar(calib_covar_arg),
+  calib_ac(calib_ac_arg),
+  symbol_table(symbol_table_arg)
 {
-  if (!symbol_table.Exist(EstimParams->name))
-    {
-      string msg = "Unknown symbol: "+EstimParams->name;
-      error(msg.c_str());
-    }
-  if (symbol_table.isReferenced(EstimParams->name) == eNotReferenced)
-    {
-      return;
-    }
-  if ((EstimParams->init_val).size() == 0)
+}
+
+void
+CalibVarStatement::writeOutput(ostream &output) const
+{
+
+  output << interfaces::comment() << "\n" << interfaces::comment() << "CALIB_VAR \n"
+         << interfaces::comment() << "\n";
+
+  for(int i = 1; i < 4 ; i++)
     {
-      EstimParams->init_val = EstimParams->mean;
+      output << "calib_var_index{" << i << "} = [];\n";
+      output << "calib_targets{" << i << "} = [];\n";
+      output << "calib_weights{" << i << "}=[];\n";
     }
-  if (EstimParams->type < 3)
+
+  // Print calibration variances
+  for(calib_var_type::const_iterator it = calib_var.begin();
+      it != calib_var.end(); it++)
     {
-      if (symbol_table.getType(EstimParams->name) == eExogenous)
-        {
-          *output << "tmp1 = find(estim_params_.var_exo(:,1)==" << symbol_table.getID(EstimParams->name)+1 << ");\n";
-          *output << "estim_params_.var_exo(tmp1,3) = " << EstimParams->low_bound << ";\n";
-          *output << "estim_params_.var_exo(tmp1,4) = " << EstimParams->up_bound << ";\n";
-        }
-      else if (symbol_table.getType(EstimParams->name) == eEndogenous)
+      const string &name = it->first;
+      const string &weight = it->second.first;
+      const string &expression = it->second.second;
+
+      int id = symbol_table.getID(name) + 1;
+      if (symbol_table.getType(name) == eEndogenous)
         {
-          *output << "tmp1 = find(estim_params_.var_endo(:,1)==" << symbol_table.getID(EstimParams->name)+1 << ");\n";
-          *output << "estim_params_.var_endo(tmp1,3) = " << EstimParams->low_bound << ";\n";
-          *output << "estim_params_.var_endo(tmp1,4) = " << EstimParams->up_bound << ";\n";
+          output << "calib_var_index{1} = [calib_var_index{1};" <<  id << "," << id << "];\n";
+          output << "calib_weights{1} = [calib_weights{1}; " << weight << "];\n";
+          output << "calib_targets{1} =[calib_targets{1}; " << expression << "];\n";
         }
-      else if (symbol_table.getType(EstimParams->name) == eParameter)
+      else if (symbol_table.getType(name) == eExogenous)
         {
-          *output << "tmp1 = find(estim_params_.param_vals(:,1)==" << symbol_table.getID(EstimParams->name)+1 << ");\n";
-          *output << "estim_params_.param_vals(tmp1,3) = " << EstimParams->low_bound << ";\n";
-          *output << "estim_params_.param_vals(tmp1,4) = " << EstimParams->up_bound << ";\n";
+          output << "calib_var_index{3} = [calib_var_index{3};" <<  id << "," << id << "];\n";
+          output << "calib_weights{3} = [calib_weights{3}; " << weight << "];\n";
+          output << "calib_targets{3} =[calib_targets{3}; " << expression << "];\n";
         }
     }
-  else
+
+  // Print calibration covariances
+  for(calib_covar_type::const_iterator it = calib_covar.begin();
+      it != calib_covar.end(); it++)
     {
-      if (symbol_table.getType(EstimParams->name) == eExogenous)
+      const string &name1 = it->first.first;
+      const string &name2 = it->first.second;
+      const string &weight = it->second.first;
+      const string &expression = it->second.second;
+
+      int id1 = symbol_table.getID(name1) + 1;
+      int id2 = symbol_table.getID(name2) + 1;
+      if (symbol_table.getType(name1) == eEndogenous)
         {
-          *output << "tmp1 = find((estim_params_.corrx(:,1)==" << symbol_table.getID(EstimParams->name)+1 << ")) & (estim_params_.corrx(:,2)==" << symbol_table.getID(EstimParams->name2)+1 << ");\n";
-          *output << "estim_params_.corrx(tmp1,4) = " << EstimParams->low_bound << ";\n";
-          *output << "estim_params_.corrx(tmp1,5) = " << EstimParams->up_bound << ";\n";
+          output << "calib_var_index{1} = [calib_var_index{1};" <<  id1 << "," << id2 << "];\n";
+          output << "calib_weights{1} = [calib_weights{1}; " << weight << "];\n";
+          output << "calib_targets{1} =[calib_targets{1}; " << expression << "];\n";
         }
-      else if (symbol_table.getType(EstimParams->name) == eEndogenous)
+      else if (symbol_table.getType(name1) == eExogenous)
         {
-          *output << "tmp1 = find((estim_params_.corrn(:,1)==" << symbol_table.getID(EstimParams->name)+1 << ")) & (estim_params_.corrn(:,2)==" << symbol_table.getID(EstimParams->name2)+1 << ";\n";
-          *output << "estim_params_.corrx(tmp1,4) = " << EstimParams->low_bound << ";\n";
-          *output << "estim_params_.corrx(tmp1,5) = " << EstimParams->up_bound << ";\n";
+          output << "calib_var_index{3} = [calib_var_index{3};" <<  id1 << "," << id2 << "];\n";
+          output << "calib_weights{3} = [calib_weights{3}; " << weight << "];\n";
+          output << "calib_targets{3} =[calib_targets{3}; " << expression << "];\n";
         }
     }
-  EstimParams->clear();
-}
 
-//-----------------------------------------------------------------------
-void ComputingTasks::set_trend_element (string name, string expression)
-{
-  //Testing if symbol exists
-  if (!symbol_table.Exist(name))
-    {
-      string msg = "Unknown variable: " + name;
-      (* error) (msg.c_str());
-    }
-  Type type = symbol_table.getType(name);
-  // int id = symbol_table.getID(name);
-  if (type == eEndogenous)
-    {
-      *output << "tmp1 = strmatch('" << name << "',options_.varobs,'exact');\n";
-      *output << "options_.trend_coeffs{tmp1} = '" << expression << "';\n";
-    }
-  else
+  // Print calibration autocorrelations
+  int max_iar = 3;
+
+  for(calib_ac_type::const_iterator it = calib_ac.begin();
+      it != calib_ac.end(); it++)
     {
-      cout << "Error : Non-variable symbol used in TREND_COEFF: " << name << endl;
+      const string &name = it->first.first;
+      int iar = it->first.second + 3;
+      const string &weight = it->second.first;
+      const string &expression = it->second.second;
+
+      int id = symbol_table.getID(name) + 1;
+
+      if (iar > max_iar)
+        {
+          // Create new variables
+          for(int i = max_iar + 1; i <= iar; i++)
+            {
+              output << "calib_var_index{" << i << "} = [];\n";
+              output << "calib_targets{" << i << "} = [];\n";
+              output << "calib_weights{" << i << "}=[];\n";
+            }
+          max_iar = iar;
+        }
+
+      output << "calib_var_index{" << iar << "} = [calib_var_index{" << iar << "};" <<  id << "];\n";
+      output << "calib_weights{" << iar << "} = [calib_weights{" << iar << "}; " << weight << "];\n";
+      output << "calib_targets{" << iar << "} =[calib_targets{" << iar << "}; " << expression << "];\n";
     }
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::BeginCalibVar(void)
+CalibStatement::CalibStatement(int covar_arg) : covar(covar_arg)
 {
-
-  *output << interfaces::comment() << "\n" << interfaces::comment() << "CALIB_VAR \n"
-          << interfaces::comment() << "\n";
-  for(int i=1;i<4;++i)
-    {
-      *output << "calib_var_index{" << i << "} = [];\n";
-      *output << "calib_targets{" << i << "} = [];\n";
-      *output << "calib_weights{" << i << "}=[];\n";
-    }
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::setCalibVar(string name, string weight, string expression)
+void
+CalibStatement::writeOutput(ostream &output) const
 {
-  if (!symbol_table.Exist(name))
-    {
-      string msg = "calib_var: " + name + " doesn't exist";
-      error(msg.c_str());
-    }
-  int id = symbol_table.getID(name) + 1;
-  if (symbol_table.getType(name) == eEndogenous)
-    {
-      *output << "calib_var_index{1} = [calib_var_index{1};" <<  id << "," << id << "];\n";
-      *output << "calib_weights{1} = [calib_weights{1}; " << weight << "];\n";
-      *output << "calib_targets{1} =[calib_targets{1}; " << expression << "];\n";
-    }
-  else if (symbol_table.getType(name) == eExogenous)
-    {
-      *output << "calib_var_index{3} = [calib_var_index{3};" <<  id << "," << id << "];\n";
-      *output << "calib_weights{3} = [calib_weights{3}; " << weight << "];\n";
-      *output << "calib_targets{3} =[calib_targets{3}; " << expression << "];\n";
-    }
-  else
-    {
-      string msg = "calib_var: " + name + "isn't a endogenous or an exogenous variable";
-      error(msg.c_str());
-    }
+  output << "M_.Sigma_e=calib(calib_var_index,calib_targets,calib_weights," << covar << ",Sigma_e_);\n";
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::setCalibVar(string name1, string name2, string weight, string expression)
+OsrParamsStatement::OsrParamsStatement(const TmpSymbolTable &tmp_symbol_table_arg) :
+  tmp_symbol_table(tmp_symbol_table_arg)
 {
-  if (!symbol_table.Exist(name1))
-    {
-      string msg = "calib_var: " + name1 + " doesn't exist";
-      error(msg.c_str());
-    }
-  if (!symbol_table.Exist(name2))
-    {
-      string msg = "calib_var: " + name2 + " doesn't exist";
-      error(msg.c_str());
-    }
-  if (symbol_table.getType(name1) != symbol_table.getType(name2))
-    {
-      string msg = "calib_var: " + name1 + " and " + name2 + " don't have the same type";
-      error(msg.c_str());
-    }
-  int id1 = symbol_table.getID(name1) + 1;
-  int id2 = symbol_table.getID(name2) + 1;
-  if (symbol_table.getType(name1) == eEndogenous)
-    {
-      *output << "calib_var_index{1} = [calib_var_index{1};" <<  id1 << "," << id2 << "];\n";
-      *output << "calib_weights{1} = [calib_weights{1}; " << weight << "];\n";
-      *output << "calib_targets{1} =[calib_targets{1}; " << expression << "];\n";
-    }
-  else if (symbol_table.getType(name1) == eExogenous)
-    {
-      *output << "calib_var_index{3} = [calib_var_index{3};" <<  id1 << "," << id2 << "];\n";
-      *output << "calib_weights{3} = [calib_weights{3}; " << weight << "];\n";
-      *output << "calib_targets{3} =[calib_targets{3}; " << expression << "];\n";
-    }
-  else
-    {
-      string msg = "calib_var: " + name1 + " and " + name2 + "aren't endogenous or exogenous variables";
-      error(msg.c_str());
-    }
 }
 
-void ComputingTasks::setCalibAc(string name, string ar, string weight, string expression)
+void
+OsrParamsStatement::writeOutput(ostream &output) const
 {
-  int max_iar = 3;
-  if (!symbol_table.Exist(name))
-    {
-      string msg = "calib_var: " + name + " doesn't exist";
-      error(msg.c_str());
-    }
-  int id = symbol_table.getID(name) + 1;
-  int iar = atoi(ar.c_str())+3;
-  if (iar > max_iar)
-    {
-      // creates new variables
-      for(int i=max_iar+1; i <= iar; ++i)
-        {
-          *output << "calib_var_index{" << i << "} = [];\n";
-          *output << "calib_targets{" << i << "} = [];\n";
-          *output << "calib_weights{" << i << "}=[];\n";
-        }
-      max_iar = iar;
-    }
-  if (symbol_table.getType(name) == eEndogenous)
-    {
-      *output << "calib_var_index{" << iar << "} = [calib_var_index{" << iar << "};" <<  id << "];\n";
-      *output << "calib_weights{" << iar << "} = [calib_weights{" << iar << "}; " << weight << "];\n";
-      *output << "calib_targets{" << iar << "} =[calib_targets{" << iar << "}; " << expression << "];\n";
-    }
-  else
-    {
-      string msg = "calib_var: " + name + "isn't a endogenous variable";
-      error(msg.c_str());
-    }
+  tmp_symbol_table.writeOutput("osr_params_", output);
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::runCalib(int cova)
+OsrStatement::OsrStatement(const TmpSymbolTable &tmp_symbol_table_arg,
+                           const OptionsList &options_list_arg) :
+  tmp_symbol_table(tmp_symbol_table_arg),
+  options_list(options_list_arg)
 {
-  *output << "M_.Sigma_e=calib(calib_var_index,calib_targets,calib_weights," << cova << ",Sigma_e_);\n";
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::setOsrParams(string tmp)
+void
+OsrStatement::writeOutput(ostream &output) const
 {
-  *output << tmp;
+  options_list.writeOutput(output);
+  tmp_symbol_table.writeOutput("var_list_", output);
+  output << "osr(var_list_,osr_params_,obj_var_,optim_weights_);\n";
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::runOsr(string tmp1)
+OlrInstStatement::OlrInstStatement(const TmpSymbolTable &tmp_symbol_table_arg) :
+  tmp_symbol_table(tmp_symbol_table_arg)
 {
-  *output << tmp1;
-  *output << "osr(var_list_,osr_params_,obj_var_,optim_weights_);\n";
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::setOlrInst(string tmp)
+void
+OlrInstStatement::writeOutput(ostream &output) const
 {
-  *output << tmp;
+  tmp_symbol_table.writeOutput("options_.olr_inst", output);
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::runOlr(string tmp1)
+
+OlrStatement::OlrStatement(const TmpSymbolTable &tmp_symbol_table_arg,
+                           const OptionsList &options_list_arg) :
+  tmp_symbol_table(tmp_symbol_table_arg),
+  options_list(options_list_arg)
 {
-  *output << tmp1;
-  *output << "options_.olr = 1;\n";
-  *output << "options_.olr_w = optim_weights_;\n";
-  *output << "options_.olr_inst = olr_inst_;\n";
-  *output << "info = stoch_simul(var_list_);\n";
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::BeginOptimWeights(void)
+void
+OlrStatement::writeOutput(ostream &output) const
 {
-  *output << interfaces::comment() << "OPTIM_WEIGHTS\n\n";
-  *output << "optim_weights_ = sparse(M_.endo_nbr,M_.endo_nbr);\n";
-  *output << "obj_var_ = [];\n\n";
+  options_list.writeOutput(output);
+  tmp_symbol_table.writeOutput("var_list_", output);
+  output << "options_.olr = 1;\n";
+  output << "options_.olr_w = optim_weights_;\n";
+  output << "options_.olr_inst = olr_inst_;\n";
+  output << "info = stoch_simul(var_list_);\n";
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::setOptimWeights(string name, string exp)
+OptimWeightsStatement::OptimWeightsStatement(const var_weights_type &var_weights_arg,
+                                             const covar_weights_type &covar_weights_arg,
+                                             const SymbolTable &symbol_table_arg) :
+  var_weights(var_weights_arg),
+  covar_weights(covar_weights_arg),
+  symbol_table(symbol_table_arg)
 {
-  if (!symbol_table.Exist(name) || symbol_table.getType(name) != eEndogenous)
-    {
-      string msg = "optim_weights: " + name + " isn't an endogenous variable";
-      error(msg.c_str());
-    }
-  int id = symbol_table.getID(name) + 1;
-  *output <<  "optim_weights_(" << id << "," << id << ") = " << exp << ";\n";
-  *output << "obj_var_ = [obj_var_; " << id << "];\n";
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::setOptimWeights(string name1, string name2, string exp)
+void
+OptimWeightsStatement::writeOutput(ostream &output) const
 {
-  if (!symbol_table.Exist(name1) || symbol_table.getType(name1) != eEndogenous)
+  output << interfaces::comment() << "OPTIM_WEIGHTS\n\n";
+  output << "optim_weights_ = sparse(M_.endo_nbr,M_.endo_nbr);\n";
+  output << "obj_var_ = [];\n\n";
+
+  for(var_weights_type::const_iterator it = var_weights.begin();
+      it != var_weights.end(); it++)
     {
-      string msg = "optim_weights: " + name1 + " isn't an endogenous variable";
-      error(msg.c_str());
+      const string &name = it->first;
+      const string &value = it->second;
+      int id = symbol_table.getID(name) + 1;
+      output <<  "optim_weights_(" << id << "," << id << ") = " << value << ";\n";
+      output << "obj_var_ = [obj_var_; " << id << "];\n";
     }
-  if (!symbol_table.Exist(name2) || symbol_table.getType(name2) != eEndogenous)
+
+  for(covar_weights_type::const_iterator it = covar_weights.begin();
+      it != covar_weights.end(); it++)
     {
-      string msg = "optim_weights: " + name2 + " isn't an endogenous variable";
-      error(msg.c_str());
+      const string &name1 = it->first.first;
+      const string &name2 = it->first.second;
+      const string &value = it->second;
+      int id1 = symbol_table.getID(name1) + 1;
+      int id2 = symbol_table.getID(name2) + 1;
+      output <<  "optim_weights_(" << id1 << "," << id2 << ") = " << value << ";\n";
+      output << "obj_var_ = [obj_var_; " << id1 << " " << id2 << "];\n";
     }
-  int id1 = symbol_table.getID(name1) + 1;
-  int id2 = symbol_table.getID(name2) + 1;
-  *output <<  "optim_weights_(" << id1 << "," << id2 << ") = " << exp << ";\n";
-  *output << "obj_var_ = [obj_var_; " << id1 << " " << id2 << "];\n";
 }
 
-//------------------------------------------------------------------------------
-void ComputingTasks::runDynasave(string filename, string ext, string varlist)
+DynaSaveStatement::DynaSaveStatement(const TmpSymbolTable &tmp_symbol_table_arg,
+                                     const string &filename_arg, const string &ext_arg) :
+  tmp_symbol_table(tmp_symbol_table_arg),
+  filename(filename_arg),
+  ext(ext_arg)
 {
-  *output << varlist;
-  *output << "dynasave(" << filename;
-  if (ext.size() > 0)
-    {
-      *output << "," << ext;
-    }
-  *output << ",varlist_);\n";
 }
 
-void ComputingTasks::runDynatype(string filename, string ext, string varlist)
+void
+DynaSaveStatement::writeOutput(ostream &output) const
 {
-  *output << varlist;
-  *output << "dynatype(" << filename;
+  tmp_symbol_table.writeOutput("var_list_", output);
+  output << "dynasave(" << filename;
   if (ext.size() > 0)
-    {
-      *output << "," << ext;
-    }
-  *output << ",varlist_);\n";
+    output << "," << ext;
+  output << ",var_list_);\n";
 }
 
-void ComputingTasks::beginModelComparison(void)
+DynaTypeStatement::DynaTypeStatement(const TmpSymbolTable &tmp_symbol_table_arg,
+                                     const string &filename_arg, const string &ext_arg) :
+  tmp_symbol_table(tmp_symbol_table_arg),
+  filename(filename_arg),
+  ext(ext_arg)
 {
-  *output << "ModelNames_ = {};\n";
-  *output << "ModelPriors_ = {};\n";
 }
 
-void ComputingTasks::addMcFilename(string filename, string prior)
+void
+DynaTypeStatement::writeOutput(ostream &output) const
 {
-  *output << "ModelNames_ = { ModelNames_{:} '" << filename << "};\n";
-  *output << "ModelPriors_ = { ModelPriors_{:} '" << prior << "};\n";
+  tmp_symbol_table.writeOutput("var_list_", output);
+  output << "dynatype(" << filename;
+  if (ext.size() > 0)
+    output << "," << ext;
+  output << ",var_list_);\n";
 }
 
-void ComputingTasks::runModelComparison(void)
+ModelComparisonStatement::ModelComparisonStatement(const filename_list_type &filename_list_arg,
+                                                   const OptionsList &options_list_arg) :
+  filename_list(filename_list_arg),
+  options_list(options_list_arg)
 {
-  *output << "model_comparison(ModelNames_,ModelPriors_);\n";
 }
 
-/*
-  string ComputingTasks::get(void)
-  {
-  return output.str();
-  }
-*/
+void
+ModelComparisonStatement::writeOutput(ostream &output) const
+{
+  options_list.writeOutput(output);
+
+  output << "ModelNames_ = {};\n";
+  output << "ModelPriors_ = {};\n";
+
+  for(filename_list_type::const_iterator it = filename_list.begin();
+      it != filename_list.end(); it++)
+    {
+      output << "ModelNames_ = { ModelNames_{:} '" << it->first << "};\n";
+      output << "ModelPriors_ = { ModelPriors_{:} '" << it->second << "};\n";
+    }
+  output << "model_comparison(ModelNames_,ModelPriors_);\n";
+}
diff --git a/parser.src/DynareBison.yy b/parser.src/DynareBison.yy
index dac38c54a1a029b7bb32305c7b3e6786bb419ebb..eb9d7e7c7d7b67d8605fca08f8ceb73e0f92b18c 100644
--- a/parser.src/DynareBison.yy
+++ b/parser.src/DynareBison.yy
@@ -100,7 +100,7 @@ typedef pair<int, Type> ExpObj;
  	| initval
  	| endval
  	| histval
- 	| equality_expression
+ 	| init_param
  	| shocks
  	| mshocks
  	| sigma_e
@@ -140,8 +140,8 @@ typedef pair<int, Type> ExpObj;
  	;
 
  	
- dsample : DSAMPLE INT_NUMBER ';' {driver.option_num("dsample", $2);}
-         | DSAMPLE INT_NUMBER INT_NUMBER ';' {driver.option_num("dsample", $2, $3);}
+ dsample : DSAMPLE INT_NUMBER ';' { driver.dsample($2);}
+         | DSAMPLE INT_NUMBER INT_NUMBER ';' {driver.dsample($2, $3);}
          ; 
 
  rplot : RPLOT tmp_var_list ';' {driver.rplot();}
@@ -225,18 +225,16 @@ typedef pair<int, Type> ExpObj;
  periods 
  	: PERIODS INT_NUMBER ';'
  		{
-		 driver.option_num("periods", $2);
-		 driver.option_num("simul", "1");
+      driver.periods($2);
 		}
     | PERIODS EQUAL INT_NUMBER ';'
  		{
-		 driver.option_num("periods", $3);
-		 driver.option_num("simul", "1");
+      driver.periods($3);
 		}
     ;
 
  		   
- equality_expression
+ init_param
  	: NAME EQUAL expression ';'
     {driver.init_param($1, $3);} 
 	;
@@ -297,12 +295,12 @@ typedef pair<int, Type> ExpObj;
           {$$ = driver.add_expression_token($1, $3, token::COMMA);}
 
  initval
- 	: INITVAL ';' {driver.begin_initval();} initval_list END
+ 	: INITVAL ';' initval_list END
  		{driver.end_initval();}
  	;
  	
  endval
- 	: ENDVAL ';' {driver.begin_endval();} initval_list END
+ 	: ENDVAL ';' initval_list END
  		{driver.end_endval();}
  	; 	
 
@@ -317,7 +315,9 @@ typedef pair<int, Type> ExpObj;
  	; 	
  	
  histval
- 	: HISTVAL ';' {driver.begin_histval();} histval_list END 
+  : HISTVAL ';' histval_list END
+    { driver.end_histval(); }
+  ;
 
  histval_list
  	: histval_list histval_elem
@@ -331,9 +331,9 @@ typedef pair<int, Type> ExpObj;
 	
  model
  	: MODEL ';' equation_list END 
- 	| MODEL '(' LINEAR ')' ';' {driver.option_num("linear","1");} 
+ 	| MODEL '(' o_linear ')' ';'
 		equation_list END
- 	| MODEL '(' USE_DLL ')' ';' {driver.use_dll();} 
+ 	| MODEL '(' USE_DLL ')' ';' {driver.use_dll();}
 		equation_list END
  	;
 
@@ -405,11 +405,11 @@ typedef pair<int, Type> ExpObj;
 	;
 	
  shocks
- 	: SHOCKS  ';' {driver.begin_shocks();} shock_list END {driver.end_shocks();} 
+ 	: SHOCKS  ';' shock_list END {driver.end_shocks();}
  	;
 
  mshocks
- 	: MSHOCKS  ';' {driver.begin_mshocks();} shock_list END {driver.end_shocks();} 
+  : MSHOCKS  ';' shock_list END {driver.end_mshocks();}
  	;
 
  shock_list 
@@ -612,14 +612,15 @@ typedef pair<int, Type> ExpObj;
  	;
 
  estimated_params 
-	: ESTIMATED_PARAMS ';' {driver.estimation_init();} estimated_list END
+	: ESTIMATED_PARAMS ';' estimated_list END
+    { driver.estimated_params(); }
 	;
 	
  estimated_list 
 	: estimated_list estimated_elem 
-		{driver.set_estimated_elements();}
+		{driver.add_estimated_params_element();}
 	| estimated_elem 
-		{driver.set_estimated_elements();}
+		{driver.add_estimated_params_element();}
 	;
 
  estimated_elem 
@@ -723,12 +724,13 @@ typedef pair<int, Type> ExpObj;
 	;
 
  estimated_params_init: ESTIMATED_PARAMS_INIT ';' estimated_init_list END
+                        { driver.estimated_params_init(); }
                       ;
 
  estimated_init_list : estimated_init_list estimated_init_elem
-                       {driver.set_estimated_init_elements();}
+                       {driver.add_estimated_params_element();}
                      | estimated_init_elem
-                       {driver.set_estimated_init_elements();}
+                       {driver.add_estimated_params_element();}
                      ;
 
  estimated_init_elem : STDERR NAME COMMA value ';'
@@ -757,12 +759,13 @@ typedef pair<int, Type> ExpObj;
                      ;
 
  estimated_params_bounds: ESTIMATED_PARAMS_BOUNDS ';' estimated_bounds_list END
+                          { driver.estimated_params_bounds(); }
                       ;
 
  estimated_bounds_list : estimated_bounds_list estimated_bounds_elem
-                       {driver.set_estimated_bounds_elements();}
+                       {driver.add_estimated_params_element();}
                      | estimated_bounds_elem
-                       {driver.set_estimated_bounds_elements();}
+                       {driver.add_estimated_params_element();}
                      ;
 
  estimated_bounds_elem : STDERR NAME COMMA value COMMA value ';'
@@ -914,8 +917,8 @@ typedef pair<int, Type> ExpObj;
 	;
 
  list_optim_option
- 	: '\'' NAME '\'' COMMA '\'' NAME '\'' {driver.optim_options($2, $6, 2);}
-	| '\'' NAME '\'' COMMA value {driver.optim_options($2, $5, 2);}
+ 	: '\'' NAME '\'' COMMA '\'' NAME '\'' {driver.optim_options_string($2, $6);}
+	| '\'' NAME '\'' COMMA value {driver.optim_options_num($2, $5);}
 	;
 
  optim_options
@@ -929,7 +932,8 @@ typedef pair<int, Type> ExpObj;
 	;
 
  observation_trends
-        : OBSERVATION_TRENDS ';' {driver.set_trend_init();} trend_list END
+        : OBSERVATION_TRENDS ';' trend_list END
+{ driver.set_trends(); }
 	;
 
  trend_list 
@@ -944,7 +948,8 @@ typedef pair<int, Type> ExpObj;
  unit_root_vars : UNIT_ROOT_VARS tmp_var_list ';' {driver.set_unit_root_vars();}
                 ;
 
- optim_weights : OPTIM_WEIGHTS ';' {driver.begin_optim_weights();} optim_weights_list END
+ optim_weights : OPTIM_WEIGHTS ';' optim_weights_list END
+                 { driver.optim_weights(); }
                ;
 
  optim_weights_list : optim_weights_list NAME expression ';' 
@@ -983,7 +988,8 @@ typedef pair<int, Type> ExpObj;
  olr_inst : OLR_INST tmp_var_list ';' {driver.set_olr_inst();}
           ;
 
- calib_var : CALIB_VAR ';' {driver.begin_calib_var();} calib_var_list END
+ calib_var : CALIB_VAR ';' calib_var_list END
+             { driver.run_calib_var(); }
            ;
 
  calib_var_list : calib_var_list calib_arg1
@@ -991,7 +997,7 @@ typedef pair<int, Type> ExpObj;
                 ;
 
  calib_arg1 : NAME calib_arg2 EQUAL expression ';' {driver.set_calib_var($1, $2, $4);}
-            | NAME COMMA NAME calib_arg2 EQUAL expression ';' {driver.set_calib_var($1, $3, $4, $6);}
+            | NAME COMMA NAME calib_arg2 EQUAL expression ';' {driver.set_calib_covar($1, $3, $4, $6);}
             | AUTOCORR NAME '(' INT_NUMBER ')' calib_arg2 EQUAL expression ';' {driver.set_calib_ac($2, $4, $6, $8);}
             ;
 
@@ -1018,8 +1024,8 @@ typedef pair<int, Type> ExpObj;
           | DYNASAVE '(' NAME '.' NAME ')' tmp_var_list ';' {driver.run_dynasave($3, $5);}
           | DYNASAVE NAME '.' NAME ';' {driver.run_dynasave($2, $4);};
 
- model_comparison : MODEL_COMPARISON '(' model_comparison_options ')' {driver.begin_model_comparison();} 
-                       filename_list ';' {driver.run_model_comparison();}
+ model_comparison : MODEL_COMPARISON '(' model_comparison_options ')' filename_list ';'
+                    {driver.run_model_comparison();}
                   ;
 
  model_comparison_options: model_comparison_options COMMA model_comparison_option
@@ -1051,7 +1057,7 @@ typedef pair<int, Type> ExpObj;
  o_dr_algo: DR_ALGO EQUAL INT_NUMBER {driver.option_num("dr_algo", $3);};
  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_linear: LINEAR {driver.option_num("linear", "1");};
+ o_linear: LINEAR {driver.linear();};
  o_order: ORDER EQUAL INT_NUMBER {driver.option_num("order", $3);};
  o_replic: REPLIC EQUAL INT_NUMBER {driver.option_num("replic", $3);};
  o_drop: DROP EQUAL INT_NUMBER {driver.option_num("drop", $3);};
diff --git a/parser.src/DynareFlex.ll b/parser.src/DynareFlex.ll
index f62147895a3eeb20cf765041c1676ee892212836..a9f4773213d66b6364bfc75d7c2cabb0368ab82d 100644
--- a/parser.src/DynareFlex.ll
+++ b/parser.src/DynareFlex.ll
@@ -67,7 +67,6 @@ int sigma_e = 0;
 <INITIAL>rplot	 	{BEGIN DYNARE_STATEMENT; return token::RPLOT;}
 <INITIAL>osr_params 	{BEGIN DYNARE_STATEMENT; return token::OSR_PARAMS;}
 <INITIAL>osr	 	{BEGIN DYNARE_STATEMENT; return token::OSR;}
-<INITIAL>calib_var 	{BEGIN DYNARE_STATEMENT; return token::CALIB_VAR;}
 <INITIAL>dynatype	{BEGIN DYNARE_STATEMENT; return token::DYNATYPE;}
 <INITIAL>dynasave 	{BEGIN DYNARE_STATEMENT; return token::DYNASAVE;}
 <INITIAL>olr	 	{BEGIN DYNARE_STATEMENT; return token::OLR;}
@@ -80,6 +79,7 @@ int sigma_e = 0;
 <INITIAL>stoch_simul {BEGIN DYNARE_STATEMENT; return token::STOCH_SIMUL;}
 <INITIAL>dsample {BEGIN DYNARE_STATEMENT; return token::DSAMPLE;}
 <INITIAL>Sigma_e {BEGIN DYNARE_STATEMENT; sigma_e = 1; return token::SIGMA_E;}
+<INITIAL>calib {BEGIN DYNARE_STATEMENT; return token::CALIB;}
 
  /* End of a Dynare statement */
 <DYNARE_STATEMENT>; {
@@ -102,6 +102,7 @@ int sigma_e = 0;
 <INITIAL>estimated_params_bounds 	{BEGIN DYNARE_BLOCK; return token::ESTIMATED_PARAMS_BOUNDS;}
 <INITIAL>observation_trends {BEGIN DYNARE_BLOCK; return token::OBSERVATION_TRENDS;}
 <INITIAL>optim_weights {BEGIN DYNARE_BLOCK; return token::OPTIM_WEIGHTS;}
+<INITIAL>calib_var 	{BEGIN DYNARE_BLOCK; return token::CALIB_VAR;}
 
  /* End of a Dynare block */
 <DYNARE_BLOCK>end[ \t\n]*; 	{BEGIN INITIAL; return token::END;}   
@@ -154,6 +155,7 @@ int sigma_e = 0;
 <DYNARE_STATEMENT>modifiedharmonicmean {return token::MODIFIEDHARMONICMEAN;}
 <DYNARE_STATEMENT>constant	{return token::CONSTANT;}
 <DYNARE_STATEMENT>noconstant	{return token::NOCONSTANT;}
+<DYNARE_STATEMENT>covar {return token::COVAR;}
 
 <DYNARE_STATEMENT>[\$][^$]*[\$] {
   strtok(yytext+1, "$");
@@ -179,6 +181,7 @@ int sigma_e = 0;
 <DYNARE_BLOCK>; {return yy::parser::token_type (yytext[0]);}
 <DYNARE_BLOCK># {return yy::parser::token_type (yytext[0]);}
 
+<DYNARE_BLOCK>autocorr {return token::AUTOCORR;}
 
  /* Inside Dynare statement */
 <DYNARE_STATEMENT>solve_algo {return token::SOLVE_ALGO;}
@@ -195,7 +198,6 @@ int sigma_e = 0;
 <DYNARE_STATEMENT>simul_seed {return token::SIMUL_SEED;}
 <DYNARE_STATEMENT>qz_criterium {return token::QZ_CRITERIUM;}
 <DYNARE_STATEMENT>simul {return token::SIMUL;}
-<DYNARE_STATEMENT>autocorr {return token::AUTOCORR;}
 <DYNARE_STATEMENT>olr_beta {return token::OLR_BETA;}
 <DYNARE_STATEMENT>xtick   		{return token::XTICK;}  	  
 <DYNARE_STATEMENT>xticklabel   		{return token::XTICKLABEL;}  	  
@@ -267,15 +269,17 @@ int sigma_e = 0;
     }
   else
     {
+      /* Enter a native block */
       BEGIN NATIVE;
-      driver.add_native(yytext);
+      yyless(0);
     }
 }
 
-<INITIAL>. {BEGIN NATIVE; driver.add_native(yytext); }
+ /* Enter a native block */
+<INITIAL>. { BEGIN NATIVE; yyless(0); }
 
- /* NATIVE Block */
-<NATIVE>.* {BEGIN INITIAL; driver.add_native(yytext); driver.add_native("\n"); }
+ /* Add the native statement */
+<NATIVE>.* { driver.add_native(yytext); BEGIN INITIAL; }
 
 <*>.      { driver.error("Unrecognized character: '" + string(yytext) + "'"); }
 %%
diff --git a/parser.src/DynareMain.cc b/parser.src/DynareMain.cc
index 25cad9020afaefd1870485ba80039e6c0e5b2e29..295be442e8a5ca25c42f936def7620ebdb21863a 100644
--- a/parser.src/DynareMain.cc
+++ b/parser.src/DynareMain.cc
@@ -6,7 +6,6 @@
 using namespace std;
 
 #include "ParsingDriver.hh"
-#include "OutputFile.hh"
 #include "ModFile.hh"
 
 /*!
@@ -17,9 +16,6 @@ using namespace std;
 int
 main(int argc, char** argv)
 {
-  OutputFile output_file;
-  ostringstream output;
-
   if (argc < 2)
     {
       cerr << "Missing model file" << endl;
@@ -29,8 +25,7 @@ main(int argc, char** argv)
 
   ParsingDriver p;
 
-  // Sets string output of parser
-  p.setoutput(&output);
+  bool clear_all = true;
 
   // Parse options
   for (int arg = 2; arg < argc; arg++)
@@ -42,7 +37,7 @@ main(int argc, char** argv)
         }
       else
         if (string(argv[arg]) == string("noclearall"))
-          output_file.clear_all = false;
+          clear_all = false;
     }
 
   cout << "Starting Dynare ..." << endl;
@@ -51,15 +46,11 @@ main(int argc, char** argv)
   // Launch parsing
   ModFile *mod_file = p.parse(argv[1]);
 
-  // Execute final instructions
-  p.finish();
+  // FIXME
+  string basename = argv[1];
+  basename.erase(basename.size() - 4, 4);
 
-  string name = argv[1];
-  name.erase(name.size() - 4,4);
-  // Opening and init main Output file (.m or .sci file)
-  output_file.Open(name, mod_file);
-  // Writing remaining string output to output file
-  output_file.Save(output, mod_file);
+  mod_file->writeOutputFiles(basename, clear_all);
 
   delete mod_file;
 
diff --git a/parser.src/Makefile b/parser.src/Makefile
index 11185979d17c9c494614eb6466c9037fbf6a0f4f..cf04d72ed05c6218470f0cb4a6aa7660ba1970ff 100644
--- a/parser.src/Makefile
+++ b/parser.src/Makefile
@@ -43,7 +43,6 @@ COMMON_OBJ=\
 	NumericalConstants.o\
 	NumericalInitialization.o\
 	OperatorTable.o\
-	OutputFile.o\
 	Shocks.o\
 	SigmaeInitialization.o\
 	SymbolTable.o\
@@ -51,7 +50,8 @@ COMMON_OBJ=\
 	VariableTable.o\
 	ParsingDriver.o\
 	DataTree.o \
-	ModFile.o
+	ModFile.o \
+	Statement.o
 
 MATLAB_OBJ = InterfaceMatlab.o
 
diff --git a/parser.src/ModFile.cc b/parser.src/ModFile.cc
index 4e8e6bff4b286826bf08e3b9695aa20d8b866b04..7999f4c30ea7371f7efd380fdc193beb0aecd9c6 100644
--- a/parser.src/ModFile.cc
+++ b/parser.src/ModFile.cc
@@ -1,9 +1,108 @@
 #include "ModFile.hh"
+#include "Interface.hh"
 
 ModFile::ModFile() : symbol_table(model_parameters),
                      variable_table(symbol_table, model_parameters),
-                     numerical_initialization(symbol_table, model_parameters),
-                     computing_tasks(symbol_table),
-                     model_tree(symbol_table, variable_table, model_parameters, num_constants)
+                     model_tree(symbol_table, variable_table, model_parameters, num_constants),
+                     order(-1), linear(-1)
 {
 }
+
+ModFile::~ModFile()
+{
+  for(vector<Statement *>::iterator it = statements.begin();
+      it != statements.end(); it++)
+    delete (*it);
+}
+
+void
+ModFile::addStatement(Statement *st)
+{
+  statements.push_back(st);
+}
+
+void
+ModFile::writeOutputFiles(const string &basename, bool clear_all)
+{
+  ofstream mOutputFile;
+
+  if (basename.size())
+    {
+      string fname(basename);
+      fname += interfaces::function_file_extension();
+      mOutputFile.open(fname.c_str(), ios::out | ios::binary);
+      if (!mOutputFile.is_open())
+        {
+          cerr << "Error: Can't open file " << fname
+               << " for writing" << endl;
+          exit(-1);
+        }
+    }
+  else
+    {
+      cerr << "Error: Missing file name" << endl;
+      exit(-1);
+    }
+
+  mOutputFile << interfaces::comment() << endl;
+  mOutputFile << interfaces::comment() << "Status : main Dynare file " << endl;
+  mOutputFile << interfaces::comment() << endl;
+  mOutputFile << interfaces::comment() << "Warning : this file is generated automatically by Dynare" << endl;
+  mOutputFile << interfaces::comment() << "          from model file (.mod)" << endl << endl;
+
+  if (clear_all)
+    mOutputFile << "clear all" << endl;
+  mOutputFile << "tic;" << endl;
+  mOutputFile << "global M_ oo_ exedet_ exdet_ recur_ recurs_ " << endl;
+  mOutputFile << "global options_ endval_" << endl;
+  mOutputFile << "global ys0_ recurs0_ ex0_ ct_" << endl;
+  mOutputFile << "options_ = [];" << endl;
+  mOutputFile << "M_.fname = '" << basename << "';" << endl;
+  mOutputFile << interfaces::comment() << endl;
+  mOutputFile << interfaces::comment() << "Some global variables initialisation" << endl;
+  mOutputFile << interfaces::comment() << endl;
+  mOutputFile << "global_initialization;" << endl;
+  mOutputFile << "diary off;" << endl << "warning off;" << endl << endl;
+  mOutputFile << interfaces::delete_file(basename + ".log") << ";" << endl;
+  mOutputFile << "warning on;" << endl << "warning backtrace;" << endl;
+  mOutputFile << "logname_ = '" << basename << ".log';" << endl;
+  mOutputFile << "diary '" << basename << ".log';" << endl;
+
+  if (model_tree.offset == 0)
+    {
+      mOutputFile << "if ";
+      mOutputFile << interfaces::file_exist(basename + "_static.c)") << endl;
+      mOutputFile << "   clear " << basename << "_static" << endl;
+      mOutputFile << "   " << interfaces::compile(basename +"_static.c") << endl;
+      mOutputFile << "end" << endl;
+      mOutputFile << "if ";
+      mOutputFile << interfaces::file_exist(basename + "_dynamic.c)") << endl;
+      mOutputFile << "   clear " << basename << "_dynamic" << endl;
+      mOutputFile << "   " + interfaces::compile(basename+"_dynamic.c") << endl;
+      mOutputFile << "end" << endl;
+    }
+  else
+    {
+      mOutputFile << "erase_compiled_function('" + basename +"_static');" << endl;
+      mOutputFile << "erase_compiled_function('" + basename +"_dynamic');" << endl;
+      mOutputFile << interfaces::load_model_function_files(basename);
+    }
+
+  symbol_table.writeOutput(mOutputFile);
+
+  if (linear == 1)
+    mOutputFile << "options_.linear = 1;" << endl;
+
+  model_tree.writeOutput(mOutputFile, basename, order, linear);
+
+  // Print statements
+  for(vector<Statement *>::iterator it = statements.begin();
+      it != statements.end(); it++)
+    (*it)->writeOutput(mOutputFile);
+
+  mOutputFile << "save('" << basename << "_results', 'oo_');" << endl;
+  mOutputFile << "diary off" << endl;
+
+  mOutputFile << endl << "disp(['Total computing time : ' sec2hms(round(toc)) ]);" << endl;
+  mOutputFile.close();
+}
diff --git a/parser.src/ModelTree.cc b/parser.src/ModelTree.cc
index 27dadaf5c36ed104e66030788ff575385ca8400b..f172e69d42e8286c068d9772ad3765ca23098c09 100644
--- a/parser.src/ModelTree.cc
+++ b/parser.src/ModelTree.cc
@@ -38,11 +38,11 @@ ModelTree::ModelTree(SymbolTable &symbol_table_arg, VariableTable &variable_tabl
                      ModelParameters &mod_param_arg, const NumericalConstants &num_constants_arg) :
   DataTree(symbol_table_arg, variable_table_arg),
   mod_param(mod_param_arg),
-  num_constants(num_constants_arg)
+  num_constants(num_constants_arg),
+  computeHessian(false),
+  computeJacobian(false),
+  computeJacobianExo(false)
 {
-  computeJacobian = false;
-  computeJacobianExo = false;
-  computeHessian = false;
 }
 
 //------------------------------------------------------------------------------
@@ -1249,8 +1249,8 @@ inline string ModelTree::getArgument(NodeID id, Type type, EquationType iEquatio
   return argument.str();
 }
 
-//------------------------------------------------------------------------------
-void ModelTree::ModelInitialization(void)
+void
+ModelTree::ModelInitialization(ostream &output)
 {
   // Exit if there is no equation in model file*/
   if (mod_param.eq_nbr == 0)
@@ -1353,13 +1353,6 @@ void ModelTree::ModelInitialization(void)
     }
 }
 
-//------------------------------------------------------------------------------
-string ModelTree::get()
-{
-  return output.str();
-}
-
-//------------------------------------------------------------------------------
 inline int ModelTree::optimize(NodeID node)
 {
   int cost;
@@ -1390,3 +1383,41 @@ inline int ModelTree::optimize(NodeID node)
   node->tmp_status = 0;
   return tmp_status;
 }
+
+void
+ModelTree::writeOutput(ostream &output, const string &basename, int order, int linear)
+{
+  // Setting flags to compute what is necessary
+  if (order == 1 || linear == 1)
+    {
+      computeJacobianExo = true;
+      computeJacobian = false;
+    }
+  else if (order != -1 && linear != -1)
+    {
+      computeHessian = true;
+      computeJacobianExo = true;
+    }
+
+  ModelInitialization(output);
+
+  if (computeHessian)
+    derive(2);
+  else
+    derive(1);
+
+  cout << "Processing outputs ..." << endl;
+  setStaticModel();
+  setDynamicModel();
+
+  if (offset == 0)
+    {
+      OpenCFiles(basename + "_static", basename + "_dynamic");
+      SaveCFiles();
+    }
+  else
+    {
+      OpenMFiles(basename + "_static", basename + "_dynamic");
+      SaveMFiles();
+    }
+}
diff --git a/parser.src/NumericalInitialization.cc b/parser.src/NumericalInitialization.cc
index 21648cddeb3cdaccdb0b31f12fc5c9e2c0d47a8e..fce5ea2d4e560f1acf6e7f12570834d41d84e2ef 100644
--- a/parser.src/NumericalInitialization.cc
+++ b/parser.src/NumericalInitialization.cc
@@ -1,204 +1,128 @@
-/*! \file
-  \version 1.0
-  \date 04/09/2004
-  \par This file implements the NumericalInitialization class methodes.
-*/
-//------------------------------------------------------------------------------
-using namespace std;
-//------------------------------------------------------------------------------
 #include "NumericalInitialization.hh"
 #include "Interface.hh"
 
-NumericalInitialization::NumericalInitialization(const SymbolTable &symbol_table_arg,
-                                                 const ModelParameters &mod_param_arg) :
-  symbol_table(symbol_table_arg),
-  mod_param(mod_param_arg)
+InitParamStatement::InitParamStatement(const string &param_name_arg,
+                                       const string &param_value_arg,
+                                       const SymbolTable &symbol_table_arg) :
+  param_name(param_name_arg),
+  param_value(param_value_arg),
+  symbol_table(symbol_table_arg)
 {
-  //Empty
 }
 
-//------------------------------------------------------------------------------
-NumericalInitialization::~NumericalInitialization()
+void
+InitParamStatement::writeOutput(ostream &output) const
 {
-  //Empty
+  int id = symbol_table.getID(param_name) + 1;
+  output << "M_.params( " << id << " ) = " << param_value << ";\n";
+  output << param_name << " = M_.params( " << id << " );\n";
 }
 
-//------------------------------------------------------------------------------
-void NumericalInitialization::setOutput(ostringstream* iOutput)
+InitOrEndValStatement::InitOrEndValStatement(const init_values_type &init_values_arg,
+                                             const SymbolTable &symbol_table_arg) :
+  init_values(init_values_arg),
+  symbol_table(symbol_table_arg)
 {
-  output = iOutput;
 }
 
-//------------------------------------------------------------------------------
-void  NumericalInitialization::SetConstant (string name, string expression)
+void
+InitOrEndValStatement::writeInitValues(ostream &output) const
 {
-
-  //Testing if symbol exists
-  if (!symbol_table.Exist(name))
-    {
-      string msg = "Unknown parameter: " + name;
-      (* error) (msg.c_str());
-    }
-  // Testing symbol type
-  if (symbol_table.getType(name) != eParameter)
+  for(init_values_type::const_iterator it = init_values.begin();
+      it != init_values.end(); it++)
     {
-      string msg = "Non-parameter used as a parameter: " + name;
-      (* error) (msg.c_str());
+      const string &name = it->first;
+      const string &expression = it->second;
+
+      Type type = symbol_table.getType(name);
+      int id = symbol_table.getID(name) + 1;
+
+      if (type == eEndogenous)
+        output << "oo_.steady_state( " << id << " ) = " << expression << ";\n";
+      else if (type == eExogenous)
+        output << "oo_.exo_steady_state( " << id << " ) = " << expression << ";\n";
+      else if (type == eExogenousDet)
+        output << "oo_.exo_det_steady_state( " << id << " ) = " << expression << ";\n";
     }
-  // Writing expression
-  *output << "M_.params( " << symbol_table.getID(name)+1 << " ) = " << expression << ";\n";
-  *output << name << " = M_.params( " << symbol_table.getID(name)+1 << " );\n";
-  // Deleting expression
-  //TODO
-
 }
 
-//------------------------------------------------------------------------------
-void  NumericalInitialization::BeginInitval (void)
+InitValStatement::InitValStatement(const init_values_type &init_values_arg,
+                                   const SymbolTable &symbol_table_arg,
+                                   const ModelParameters &mod_param_arg) :
+  InitOrEndValStatement(init_values_arg, symbol_table_arg),
+  mod_param(mod_param_arg)
 {
+}
 
-  // Writing a Matlab comment
-  *output << interfaces::comment() << "\n" << interfaces::comment() << "INITVAL instructions \n"
-          << interfaces::comment() << "\n";
+void
+InitValStatement::writeOutput(ostream &output) const
+{
+  output << interfaces::comment() << "\n" << interfaces::comment() << "INITVAL instructions \n"
+         << interfaces::comment() << "\n";
   // Writing initval block to set initial values for variables
-  *output << "options_.initval_file = 0;\nendval_=0;\n";
+  output << "options_.initval_file = 0;\nendval_=0;\n";
 
   if (mod_param.recur_nbr > 0)
-    *output << "recurs_ = zeros(" << mod_param.recur_nbr << ", 1);\n";
-}
+    output << "recurs_ = zeros(" << mod_param.recur_nbr << ", 1);\n";
 
-//------------------------------------------------------------------------------
-void  NumericalInitialization::SetInit (string name, string expression)
-{
+  writeInitValues(output);
 
-  //Testing if symbol exists
-  if (!symbol_table.Exist(name))
-    {
-      string msg = "Unknown variable: " + name;
-      (* error) (msg.c_str());
-    }
-  Type type = symbol_table.getType(name);
-  int id = symbol_table.getID(name);
-  // Writing instrcuction that set initial or terminal value
-  // for a variable
-  if (type == eEndogenous)
-    {
-      *output << "oo_.steady_state( " << id+1 << " ) = " << expression << ";\n";
-    }
-  else if (type == eExogenous)
-    {
-      *output << "oo_.exo_steady_state( " << id+1 << " ) = " << expression << ";\n";
-    }
-  else if (type == eExogenousDet)
-    {
-      *output << "oo_.exo_det_steady_state( " << id+1 << " ) = " << expression << ";\n";
-    }
-  // Testing if symbol is a variable (Exogenousous deterministic or recursive)
-  else if ( type != eRecursiveVariable )
-    {
-      cout << "Error : Non-variable symbol used in INITVAL: " << name << endl;
-    }
+  output << "oo_.y_simul=[oo_.steady_state*ones(1,M_.maximum_lag)];\n";
+  output << "if M_.exo_nbr > 0;\n";
+  output << "\too_.exo_simul = [ones(M_.maximum_lag,1)*oo_.exo_steady_state'];\n";
+  output <<"end;\n";
+  output << "if M_.exo_det_nbr > 0;\n";
+  output << "\too_.exo_det_simul = [ones(M_.maximum_lag,1)*oo_.exo_det_steady_state'];\n";
+  output <<"end;\n";
 }
 
-//------------------------------------------------------------------------------
-void  NumericalInitialization::EndInitval(void)
+
+EndValStatement::EndValStatement(const init_values_type &init_values_arg,
+                                 const SymbolTable &symbol_table_arg) :
+  InitOrEndValStatement(init_values_arg, symbol_table_arg)
 {
-  *output << "oo_.y_simul=[oo_.steady_state*ones(1,M_.maximum_lag)];\n";
-  *output << "if M_.exo_nbr > 0;\n";
-  *output << "\too_.exo_simul = [ones(M_.maximum_lag,1)*oo_.exo_steady_state'];\n";
-  *output <<"end;\n";
-  *output << "if M_.exo_det_nbr > 0;\n";
-  *output << "\too_.exo_det_simul = [ones(M_.maximum_lag,1)*oo_.exo_det_steady_state'];\n";
-  *output <<"end;\n";
 }
 
-//------------------------------------------------------------------------------
-void  NumericalInitialization::BeginEndval (void)
+
+void
+EndValStatement::writeOutput(ostream &output) const
 {
-  // Writing a Matlab comment
-  *output << interfaces::comment() << "\n" << interfaces::comment() << "ENDVAL instructions\n"
-          << interfaces::comment() << "\n";
+  output << interfaces::comment() << "\n" << interfaces::comment() << "ENDVAL instructions\n"
+         << interfaces::comment() << "\n";
   // Writing endval block to set terminal values for variables
-  *output << "ys0_= oo_.steady_state;\nex0_ = oo_.exo_steady_state;\nrecurs0_ = recurs_;\nendval_ = 1;\n";
+  output << "ys0_= oo_.steady_state;\nex0_ = oo_.exo_steady_state;\nrecurs0_ = recurs_;\nendval_ = 1;\n";
 
+  writeInitValues(output);
 }
 
-//------------------------------------------------------------------------------
-void  NumericalInitialization::EndEndval (void)
+HistValStatement::HistValStatement(const hist_values_type &hist_values_arg,
+                                   const SymbolTable &symbol_table_arg) :
+  hist_values(hist_values_arg),
+  symbol_table(symbol_table_arg)
 {
 }
 
-//------------------------------------------------------------------------------
-void  NumericalInitialization::BeginHistval (void)
+void
+HistValStatement::writeOutput(ostream &output) const
 {
-  // Writing a Matlab comment
-  *output << interfaces::comment() << "\n" << interfaces::comment() << "HISTVAL instructions\n"
-          << interfaces::comment() << "\n";
+  output << interfaces::comment() << "\n" << interfaces::comment() << "HISTVAL instructions\n"
+         << interfaces::comment() << "\n";
 
-}
-
-//------------------------------------------------------------------------------
-void  NumericalInitialization::SetHist (string name, int lag, string expression)
-{
-  //Testing if symbol exists
-  if (!symbol_table.Exist(name))
-    {
-      string msg = "Unknown parameter: " + name;
-      (* error) (msg.c_str());
-    }
-  Type  type = symbol_table.getType(name);
-  int   id = symbol_table.getID(name);
-  // Testing symbol type
-  if (type == eEndogenous)
-    {
-      *output << "oo_.endo_simul( " << id+1 << ", M_.maximum_lag + " << lag + 1 << ") = " << expression << ";\n";
-    }
-  else if (type == eExogenous)
-    {
-      *output << "oo_.exo_simul( M_.maximum_lag + " << lag + 1 << ", " << id+1 << " ) = " << expression << ";\n";
-    }
-  // Tetsting if symbol is a variable (Exogenousous deterministic or recursive)
-  else if (type != eExogenousDet)
-    {
-      *output << "oo_.exo_det_simul( M_.maximum_lag + " << lag + 1 << ", " << id+1 << " ) = " << expression << ";\n";
-    }
-  else if (type != eRecursiveVariable)
+  for(hist_values_type::const_iterator it = hist_values.begin();
+      it != hist_values.end(); it++)
     {
-      string msg = "Non-variable symbol : " + name;
-      (* error) (msg.c_str());
+      const string &name = it->first.first;
+      const int &lag = it->first.second;
+      const string &expression = it->second;
+
+      Type type = symbol_table.getType(name);
+      int id = symbol_table.getID(name) + 1;
+
+      if (type == eEndogenous)
+        output << "oo_.endo_simul( " << id << ", M_.maximum_lag + " << lag + 1 << ") = " << expression << ";\n";
+      else if (type == eExogenous)
+        output << "oo_.exo_simul( M_.maximum_lag + " << lag + 1 << ", " << id << " ) = " << expression << ";\n";
+      else if (type != eExogenousDet)
+        output << "oo_.exo_det_simul( M_.maximum_lag + " << lag + 1 << ", " << id << " ) = " << expression << ";\n";
     }
-  // Deleting expression
-  // TODO
-
-  /////////////////////////////////
-  /*
-    char buffer[200];
-    int offset;
-
-    offset = v->var_ptr-var_list;
-    if (v->endo_exo == 1)
-    {
-
-    sprintf(buffer,"oo_.y_simul(%d,M_.maximum_lag+(%s))=",v->nbr+1,lag);
-    }
-    else if (v->endo_exo == 0)
-    {
-    initval_check1[offset] = 1;
-    sprintf(buffer,"oo_.exo_simul(M_.maximum_lag+(%s),%d)=",lag,v->nbr+1);
-    }
-
-    str_output(buffer);
-
-    p_expression(q);
-    str_output(";\n");
-  */
 }
-
-//------------------------------------------------------------------------------
-/*
-  string NumericalInitialization::get(void)
-  {
-  return output.str();
-  }
-*/
-//------------------------------------------------------------------------------
diff --git a/parser.src/OutputFile.cc b/parser.src/OutputFile.cc
deleted file mode 100644
index a6904f3e84f88605456286f89c9629700136fa55..0000000000000000000000000000000000000000
--- a/parser.src/OutputFile.cc
+++ /dev/null
@@ -1,228 +0,0 @@
-/*! \file
-  \version 1.0
-  \date 04/09/2004
-  \par This file implements the OutputFile class methodes.
-*/
-//------------------------------------------------------------------------------
-#include <iostream>
-#include <sstream>
-using namespace std;
-//------------------------------------------------------------------------------
-#include "OutputFile.hh"
-#include "SymbolTable.hh"
-#include "ModelTree.hh"
-#include "Interface.hh"
-//------------------------------------------------------------------------------
-OutputFile::OutputFile()
-{
-  clear_all = true;
-}
-
-//------------------------------------------------------------------------------
-OutputFile::~OutputFile()
-{
-  // Empty
-}
-
-//------------------------------------------------------------------------------
-void OutputFile::Open(string iFileName, ModFile *mod_file)
-{
-  if (iFileName.size())
-    {
-      string fname(iFileName);
-      fname += interfaces::function_file_extension();
-      mOutputFile.open(fname.c_str(),ios::out|ios::binary);
-      if (!mOutputFile.is_open())
-        {
-          cout << "OutputFile::Open : Error : Can't open file " << iFileName
-               << " for writing\n";
-          exit(-1);
-        }
-    }
-  else
-    {
-      cout << "OutputFile::Open : Error : Missing file name\n";
-      exit(-1);
-    }
-  mOutputFile << interfaces::comment() << "\n" << interfaces::comment();
-  mOutputFile << "Status : main Dynare file \n";
-  mOutputFile << interfaces::comment() << "\n" << interfaces::comment();
-  mOutputFile << "Warning : this file is generated automatically by Dynare\n";
-  mOutputFile << interfaces::comment();
-  mOutputFile << "          from model file (.mod)\n\n";
-  if (clear_all)
-    mOutputFile << "clear all\n";
-  mOutputFile << "tic;\n";
-  mOutputFile << "global M_ oo_ exedet_ exdet_ recur_ recurs_ \n";
-  mOutputFile << "global options_ endval_\n";
-  mOutputFile << "global ys0_ recurs0_ ex0_ ct_\n";
-  mOutputFile << "options_ = [];\n";
-  mOutputFile << "M_.fname = '" << iFileName << "';\n";
-  mOutputFile << interfaces::comment() << "\n" << interfaces::comment();
-  mOutputFile << "Some global variables initialisation\n";
-  mOutputFile << interfaces::comment() << "\n";
-  mOutputFile << "global_initialization;\n";
-  mOutputFile << "diary off;\nwarning off;\n";
-  mOutputFile << "\n" << interfaces::delete_file(iFileName + ".log;\n");
-  mOutputFile << "warning on;\nwarning backtrace;\n";
-  mOutputFile << "logname_ = '" << iFileName << ".log';\n";
-  mOutputFile << "diary '" << iFileName << ".log';\n";
-  if (mod_file->model_tree.offset == 0)
-    {
-      mOutputFile << "if ";
-      mOutputFile << interfaces::file_exist(iFileName + "_static.c)")+"\n";
-      mOutputFile << "   clear " << iFileName << "_static\n";
-      mOutputFile << "   " + interfaces::compile(iFileName +"_static.c")+"\n";
-      mOutputFile << "end\n";
-      mOutputFile << "if ";
-      mOutputFile << interfaces::file_exist(iFileName + "_dynamic.c)")+"\n";
-      mOutputFile << "   clear " << iFileName << "_dynamic\n";
-      mOutputFile << "   " + interfaces::compile(iFileName+"_dynamic.c")+"\n";
-      mOutputFile << "end\n";
-    }
-  else
-    {
-      mOutputFile << "erase_compiled_function('" + iFileName +"_static');\n";
-      mOutputFile << "erase_compiled_function('" + iFileName +"_dynamic');\n";
-      mOutputFile << interfaces::load_model_function_files(iFileName);
-    }
-}
-
-//------------------------------------------------------------------------------
-void OutputFile::Save(ostringstream& iOutput, ModFile *mod_file)
-{
-  mOutputFile << mod_file->symbol_table.get();
-  mOutputFile << mod_file->model_tree.get();
-  mOutputFile << iOutput.str();
-  mOutputFile << "\ndisp(['Total computing time : ' sec2hms(round(toc)) ]);\n";
-  mOutputFile.close();
-}
-
-//------------------------------------------------------------------------------
-#ifdef TEST_OUTPUTFILE
-#include "NumericalInitialization.h"
-#include "ComputingTasks.h"
-#include "Expression.h"
-#include "Shocks.h"
-#include "SigmaeInitialization.h"
-#include "TmpSymbolTable.h"
-int main(void)
-{
-  OutputFile outputfile;
-  SymbolTable st;
-  Expression  exp;
-  NumericalConstants numconst;
-  NumericalInitialization numinit;
-
-  outputfile.Open("Test.m");
-
-  //0
-  SymbolTable::AddSymbolDeclar("a",eExogenous);
-  //1
-  SymbolTable::AddSymbolDeclar("b",eParameter);
-  //2
-  SymbolTable::AddSymbolDeclar("c",eExogenous);
-  //3
-  SymbolTable::AddSymbolDeclar("d",eExogenousDet);
-  //3
-  SymbolTable::AddSymbolDeclar("x",eParameter);
-  //3
-  SymbolTable::AddSymbolDeclar("y",eExogenous);
-
-  numconst.AddConstant("alpha");
-  exp.AddToken(0,eExogenous,EXP);//0
-                                 //1
-  exp.AddToken(0,eParameter,0,eExogenousDet,PLUS);
-  //2
-  exp.AddToken(0,eTempResult,UMINUS);
-  //3
-  exp.AddToken(1,eExogenous,1,eTempResult,TIMES);
-  //4
-  exp.AddToken(3,eTempResult,0,eNumericalConstant,TIMES);
-  //5
-  exp.AddToken(4,eTempResult,0,eTempResult,COMMA);
-  //6
-  exp.AddToken(5,eTempResult,0,eExogenous,COMMA);
-  //6
-  exp.AddToken(6,eTempResult,"function1");
-  exp.set();
-  //cout << exp.get();
-
-  numinit.SetConstant("x","1");
-  numinit.InitInitval();
-  numinit.SetInit("y",exp.get());
-  numinit.EndInitval();
-  numinit.InitEndval();
-  numinit.EndEndval();
-  numinit.InitHistval();
-  numinit.SetHist("y",3, exp.get());
-  //cout << numinit.get();
-
-  SigmaeInitialization siginit;
-
-  siginit.AddExpression("00");
-  siginit.EndOfRow();
-  siginit.AddExpression("10");
-  siginit.AddExpression("11");
-  siginit.EndOfRow();
-  siginit.AddExpression("20");
-  siginit.AddExpression("21");
-  siginit.AddExpression("22");
-  siginit.EndOfRow();
-  siginit.AddExpression("30");
-  siginit.AddExpression("31");
-  siginit.AddExpression("32");
-  siginit.AddExpression("33");
-  siginit.EndOfRow();
-  siginit.set();
-
-  TmpSymbolTable tmp_symbol_table1, tmp_symbol_table2;
-
-  ComputingTasks computing_tasks;
-
-  computing_tasks.set();
-  computing_tasks.SetSteady();
-  computing_tasks.SetCheck();
-  computing_tasks.SetSimul();
-
-  tmp_symbol_table1.AddTempSymbol("tmp1");
-  tmp_symbol_table1.AddTempSymbol("tmp2");
-  tmp_symbol_table1.AddTempSymbol("tmp3");
-  tmp_symbol_table1.set("var_list_");
-
-  computing_tasks.SetStochSimul(tmp_symbol_table1.get());
-
-  computing_tasks.SetOption("DROP", "500");
-  computing_tasks.RunEstimation();
-  computing_tasks.SetEstimationInit();
-  computing_tasks.SetEstimation("d", "init_val", "lo_bound", "up_bound", "prior", "p1", "p2", "p3","p4");
-  computing_tasks.SetEstimation("a", "c", "init_val", "lo_bound", "up_bound", "prior", "p1", "p2", "p3", "p4");
-  computing_tasks.SetCalibInit();
-  computing_tasks.SetCalibVariance();
-  computing_tasks.SetCalibCovariance();
-  computing_tasks.SetCalibAutoCorrelation();
-  computing_tasks.SetCalib();
-
-  tmp_symbol_table1.AddTempSymbol("tmp11");
-  tmp_symbol_table1.AddTempSymbol("tmp22");
-  tmp_symbol_table1.AddTempSymbol("tmp33");
-  tmp_symbol_table1.set("varl_list_");
-  computing_tasks.SetOsr(tmp_symbol_table1.get());
-
-  tmp_symbol_table1.AddTempSymbol("tmp11");
-  tmp_symbol_table1.AddTempSymbol("tmp22");
-  tmp_symbol_table1.AddTempSymbol("tmp33");
-  tmp_symbol_table1.set("var_list_");
-  tmp_symbol_table2.AddTempSymbol("tmp4");
-  tmp_symbol_table2.AddTempSymbol("tmp5");
-  tmp_symbol_table2.AddTempSymbol("tmp6");
-  tmp_symbol_table2.set("olr_inst_");
-  computing_tasks.SetOlr(tmp_symbol_table1.get(), tmp_symbol_table2.get());
-
-  computing_tasks.SetOptimWeightsInit();
-  computing_tasks.SetOptimWeights1();
-  computing_tasks.SetOptimWeights2();
-  outputfile.Save();
-}
-#endif
-//------------------------------------------------------------------------------
diff --git a/parser.src/ParsingDriver.cc b/parser.src/ParsingDriver.cc
index a480b43b771bf59fa27f9dc1c2717c609fa9e3ae..5d31156aa3d34f36a3e07c86f68888992d6c705e 100644
--- a/parser.src/ParsingDriver.cc
+++ b/parser.src/ParsingDriver.cc
@@ -1,21 +1,8 @@
 #include "ParsingDriver.hh"
+#include "Statement.hh"
 
 ParsingDriver::ParsingDriver() : trace_scanning(false), trace_parsing(false)
 {
-  mod_file = new ModFile();
-  mod_file->order = -1;
-  mod_file->linear = -1;
-
-  mod_file->model_tree.error = error;
-  mod_file->symbol_table.error = error;
-  mod_file->variable_table.error = error;
-  mod_file->shocks.error = error;
-  mod_file->numerical_initialization.error = error;
-  mod_file->computing_tasks.error = error;
-  tmp_symbol_table.error = error;
-
-  tmp_symbol_table.setGlobalSymbolTable(&mod_file->symbol_table);
-  expression.setNumericalConstants(&mod_file->num_constants);
 }
 
 ParsingDriver::~ParsingDriver()
@@ -53,12 +40,24 @@ ParsingDriver::get_expression(ExpObj *exp)
 ModFile *
 ParsingDriver::parse(const string &f)
 {
+  mod_file = new ModFile();
+
+  mod_file->model_tree.error = error;
+  mod_file->symbol_table.error = error;
+  mod_file->variable_table.error = error;
+
+  expression.setNumericalConstants(&mod_file->num_constants);
+  tmp_symbol_table = new TmpSymbolTable(mod_file->symbol_table);
+
   file = f;
   scan_begin();
   yy::parser parser(*this);
   parser.set_debug_level(trace_parsing);
   parser.parse();
   scan_end();
+
+  delete tmp_symbol_table;
+
   return mod_file;
 }
 
@@ -77,16 +76,6 @@ ParsingDriver::error(const string &m)
   exit(-1);
 }
 
-void
-ParsingDriver::setoutput(ostringstream *ostr)
-{
-  output = ostr;
-  mod_file->numerical_initialization.setOutput(ostr);
-  mod_file->shocks.setOutput(ostr);
-  mod_file->sigmae.setOutput(ostr);
-  mod_file->computing_tasks.setOutput(ostr);
-}
-
 void
 ParsingDriver::declare_endogenous(string *name, string *tex_name)
 {
@@ -213,148 +202,166 @@ ParsingDriver::add_expression_token(ExpObj *arg1, string *op_name)
 }
 
 void
-ParsingDriver::init_param(string *name, ExpObj *rhs)
+ParsingDriver::periods(string *periods)
 {
-  mod_file->numerical_initialization.SetConstant(*name, get_expression(rhs));
-  delete name;
-  delete rhs;
+  int periods_val = atoi(periods->c_str());
+  mod_file->addStatement(new PeriodsStatement(periods_val));
+  delete periods;
 }
 
 void
-ParsingDriver::init_val(string *name, ExpObj *rhs)
+ParsingDriver::dsample(string *arg1)
 {
-  mod_file->numerical_initialization.SetInit(*name, get_expression(rhs));
-  delete name;
-  delete rhs;
+  int arg1_val = atoi(arg1->c_str());
+  mod_file->addStatement(new DsampleStatement(arg1_val));
+  delete arg1;
 }
 
 void
-ParsingDriver::hist_val(string *name, string *lag, ExpObj *rhs)
+ParsingDriver::dsample(string *arg1, string *arg2)
 {
-  int ilag = atoi(lag->c_str());
-  mod_file->numerical_initialization.SetHist(*name, ilag, get_expression(rhs));
+  int arg1_val = atoi(arg1->c_str());
+  int arg2_val = atoi(arg2->c_str());
+  mod_file->addStatement(new DsampleStatement(arg1_val, arg2_val));
+  delete arg1;
+  delete arg2;
+}
+
+void
+ParsingDriver::init_param(string *name, ExpObj *rhs)
+{
+  check_symbol_existence(*name);
+  if (mod_file->symbol_table.getType(*name) != eParameter)
+    error(*name + " is not a parameter");
+
+  mod_file->addStatement(new InitParamStatement(*name, get_expression(rhs), mod_file->symbol_table));
   delete name;
-  delete lag;
   delete rhs;
 }
 
 void
-ParsingDriver::use_dll()
+ParsingDriver::init_val(string *name, ExpObj *rhs)
 {
-  // Seetting variable momber offset to use C outputs
-  mod_file->model_tree.offset = 0;
+  check_symbol_existence(*name);
+  Type type = mod_file->symbol_table.getType(*name);
+
+  if (type != eEndogenous
+      && type != eExogenous
+      && type != eExogenousDet)
+    error("initval/endval: " + *name + " should be an endogenous or exogenous variable");
+
+  init_values.push_back(pair<string, string>(*name, get_expression(rhs)));
+
+  delete name;
+  delete rhs;
 }
 
 void
-ParsingDriver::finish()
+ParsingDriver::hist_val(string *name, string *lag, ExpObj *rhs)
 {
-  string model_file_name(file);
+  check_symbol_existence(*name);
+  Type type = mod_file->symbol_table.getType(*name);
 
-  // Setting flags to compute what is necessary
-  if (mod_file->order == 1 || mod_file->linear == 1)
-    {
-      mod_file->model_tree.computeJacobianExo = true;
-      mod_file->model_tree.computeJacobian = false;
-    }
-  else if (mod_file->order != -1 && mod_file->linear != -1)
-    {
-      mod_file->model_tree.computeHessian = true;
-      mod_file->model_tree.computeJacobianExo = true;
-    }
-  // Removing extension chars
-  model_file_name.erase(model_file_name.size()-4,4);
-  mod_file->model_tree.ModelInitialization();
+  if (type != eEndogenous
+      && type != eExogenous
+      && type != eExogenousDet)
+    error("hist_val: " + *name + " should be an endogenous or exogenous variable");
 
-  if (mod_file->model_tree.computeHessian )
-    mod_file->model_tree.derive(2);
-  else
-    mod_file->model_tree.derive(1);
+  int ilag = atoi(lag->c_str());
+  pair<string, int> key(*name, ilag);
 
-  cout << "Processing outputs ..." << endl;
-  mod_file->model_tree.setStaticModel();
-  mod_file->model_tree.setDynamicModel();
+  if (hist_values.find(key) != hist_values.end())
+    error("hist_val: (" + *name + ", " + *lag + ") declared twice");
 
-  if (mod_file->model_tree.offset == 0)
-    {
-      mod_file->model_tree.OpenCFiles(model_file_name+"_static", model_file_name+"_dynamic");
-      mod_file->model_tree.SaveCFiles();
-    }
-  else
-    {
-      mod_file->model_tree.OpenMFiles(model_file_name+"_static", model_file_name+"_dynamic");
-      mod_file->model_tree.SaveMFiles();
-    }
-
-  *output << "save('" << model_file_name << "_results', 'oo_');\n";
-  *output << "diary off\n";
+  hist_values[key] = get_expression(rhs);
 
-  //  symbol_table.erase_local_parameters();
+  delete name;
+  delete lag;
+  delete rhs;
 }
 
 void
-ParsingDriver::begin_initval()
+ParsingDriver::use_dll()
 {
-  mod_file->numerical_initialization.BeginInitval();
+  // Seetting variable momber offset to use C outputs
+  mod_file->model_tree.offset = 0;
 }
 
 void
 ParsingDriver::end_initval()
 {
-  mod_file->numerical_initialization.EndInitval();
-}
-
-void
-ParsingDriver::begin_endval()
-{
-  mod_file->numerical_initialization.BeginEndval();
+  mod_file->addStatement(new InitValStatement(init_values, mod_file->symbol_table,
+                                              mod_file->model_parameters));
+  init_values.clear();
 }
 
 void
 ParsingDriver::end_endval()
 {
-  mod_file->numerical_initialization.EndEndval();
-}
-
-void
-ParsingDriver::begin_histval()
-{
-  mod_file->numerical_initialization.BeginHistval();
+  mod_file->addStatement(new EndValStatement(init_values, mod_file->symbol_table));
+  init_values.clear();
 }
 
 void
-ParsingDriver::begin_shocks()
+ParsingDriver::end_histval()
 {
-  mod_file->shocks.BeginShocks();
+  mod_file->addStatement(new HistValStatement(hist_values, mod_file->symbol_table));
+  hist_values.clear();
 }
 
 void
-ParsingDriver::begin_mshocks()
+ParsingDriver::end_shocks()
 {
-  mod_file->shocks.BeginMShocks();
+  mod_file->addStatement(new ShocksStatement(det_shocks, var_shocks, std_shocks,
+                                             covar_shocks, corr_shocks, mod_file->symbol_table));
+  det_shocks.clear();
+  var_shocks.clear();
+  std_shocks.clear();
+  covar_shocks.clear();
+  corr_shocks.clear();
 }
 
 void
-ParsingDriver::end_shocks()
+ParsingDriver::end_mshocks()
 {
-  mod_file->shocks.EndShocks();
+  mod_file->addStatement(new MShocksStatement(det_shocks, var_shocks, std_shocks,
+                                              covar_shocks, corr_shocks, mod_file->symbol_table));
+  det_shocks.clear();
+  var_shocks.clear();
+  std_shocks.clear();
+  covar_shocks.clear();
+  corr_shocks.clear();
 }
 
 void
 ParsingDriver::add_det_shock(string *var)
 {
   check_symbol_existence(*var);
-  int id = mod_file->symbol_table.getID(*var);
-  switch (mod_file->symbol_table.getType(*var))
+  Type type = mod_file->symbol_table.getType(*var);
+  if (type != eExogenous && type != eExogenousDet)
+    error("shocks: shocks can only be applied to exogenous variables");
+
+  if (det_shocks.find(*var) != det_shocks.end())
+    error("shocks: variable " + *var + " declared twice");
+
+  if (det_shocks_periods.size() != det_shocks_values.size())
+    error("shocks: variable " + *var + ": number of periods is different from number of shock values");
+
+  vector<ShocksStatement::DetShockElement> v;
+
+  for(unsigned int i = 0; i < det_shocks_periods.size(); i++)
     {
-    case eExogenous:
-      mod_file->shocks.AddDetShockExo(id);
-      break;
-    case eExogenousDet:
-      mod_file->shocks.AddDetShockExoDet(id);
-      break;
-    default:
-      error("Shocks can only be applied to exogenous variables");
+      ShocksStatement::DetShockElement dse;
+      dse.period1 = det_shocks_periods[i].first;
+      dse.period2 = det_shocks_periods[i].second;
+      dse.value = det_shocks_values[i];
+      v.push_back(dse);
     }
+
+  det_shocks[*var] = v;
+
+  det_shocks_periods.clear();
+  det_shocks_values.clear();
   delete var;
 }
 
@@ -362,8 +369,12 @@ void
 ParsingDriver::add_stderr_shock(string *var, ExpObj *value)
 {
   check_symbol_existence(*var);
-  int id = mod_file->symbol_table.getID(*var);
-  mod_file->shocks.AddSTDShock(id, get_expression(value));
+  if (var_shocks.find(*var) != var_shocks.end()
+      || std_shocks.find(*var) != std_shocks.end())
+    error("shocks: variance or stderr of shock on " + *var + " declared twice");
+
+  std_shocks[*var] = get_expression(value);
+
   delete var;
   delete value;
 }
@@ -372,8 +383,12 @@ void
 ParsingDriver::add_var_shock(string *var, ExpObj *value)
 {
   check_symbol_existence(*var);
-  int id = mod_file->symbol_table.getID(*var);
-  mod_file->shocks.AddVARShock(id, get_expression(value));
+  if (var_shocks.find(*var) != var_shocks.end()
+      || std_shocks.find(*var) != std_shocks.end())
+    error("shocks: variance or stderr of shock on " + *var + " declared twice");
+
+  var_shocks[*var] = get_expression(value);
+
   delete var;
   delete value;
 }
@@ -383,9 +398,18 @@ ParsingDriver::add_covar_shock(string *var1, string *var2, ExpObj *value)
 {
   check_symbol_existence(*var1);
   check_symbol_existence(*var2);
-  int id1 = mod_file->symbol_table.getID(*var1);
-  int id2 = mod_file->symbol_table.getID(*var2);
-  mod_file->shocks.AddCOVAShock(id1, id2, get_expression(value));
+
+  pair<string, string> key(*var1, *var2), key_inv(*var2, *var1);
+
+  if (covar_shocks.find(key) != covar_shocks.end()
+      || covar_shocks.find(key_inv) != covar_shocks.end()
+      || corr_shocks.find(key) != corr_shocks.end()
+      || corr_shocks.find(key_inv) != corr_shocks.end())
+    error("shocks: covariance or correlation shock on variable pair (" + *var1 + ", "
+          + *var2 + ") declared twice");
+  
+  covar_shocks[key] = get_expression(value);
+
   delete var1;
   delete var2;
   delete value;
@@ -396,9 +420,18 @@ ParsingDriver::add_correl_shock(string *var1, string *var2, ExpObj *value)
 {
   check_symbol_existence(*var1);
   check_symbol_existence(*var2);
-  int id1 = mod_file->symbol_table.getID(*var1);
-  int id2 = mod_file->symbol_table.getID(*var2);
-  mod_file->shocks.AddCORRShock(id1, id2, get_expression(value));
+
+  pair<string, string> key(*var1, *var2), key_inv(*var2, *var1);
+
+  if (covar_shocks.find(key) != covar_shocks.end()
+      || covar_shocks.find(key_inv) != covar_shocks.end()
+      || corr_shocks.find(key) != corr_shocks.end()
+      || corr_shocks.find(key_inv) != corr_shocks.end())
+    error("shocks: covariance or correlation shock on variable pair (" + *var1 + ", "
+          + *var2 + ") declared twice");
+  
+  corr_shocks[key] = get_expression(value);
+
   delete var1;
   delete var2;
   delete value;
@@ -407,7 +440,9 @@ ParsingDriver::add_correl_shock(string *var1, string *var2, ExpObj *value)
 void
 ParsingDriver::add_period(string *p1, string *p2)
 {
-  mod_file->shocks.AddPeriod(*p1, *p2);
+  int p1_val = atoi(p1->c_str());
+  int p2_val = atoi(p2->c_str());
+  det_shocks_periods.push_back(pair<int, int>(p1_val, p2_val));
   delete p1;
   delete p2;
 }
@@ -415,61 +450,76 @@ ParsingDriver::add_period(string *p1, string *p2)
 void
 ParsingDriver::add_period(string *p1)
 {
-  mod_file->shocks.AddPeriod(*p1, *p1);
+  int p1_val = atoi(p1->c_str());
+  det_shocks_periods.push_back(pair<int, int>(p1_val, p1_val));
   delete p1;
 }
 
 void
 ParsingDriver::add_value(string *value)
 {
-  mod_file->shocks.AddValue(*value);
+  det_shocks_values.push_back(*value);
   delete value;
 }
 
 void
 ParsingDriver::add_value(ExpObj *value)
 {
-  mod_file->shocks.AddValue(get_expression(value));
+  det_shocks_values.push_back(get_expression(value));
   delete value;
 }
 
 void
 ParsingDriver::do_sigma_e()
 {
-  mod_file->sigmae.set();
+  try
+    {
+      mod_file->addStatement(new SigmaeStatement(sigmae_matrix));
+    }
+  catch(SigmaeStatement::MatrixFormException &e)
+    {
+      error("Sigma_e: matrix is neither upper triangular nor lower triangular");
+    }
+  sigmae_matrix.clear();
 }
 
 void
 ParsingDriver::end_of_row()
 {
-  mod_file->sigmae.EndOfRow();
+  sigmae_matrix.push_back(sigmae_row);
+  sigmae_row.clear();
 }
 
 void
 ParsingDriver::add_to_row(string *s)
 {
-  mod_file->sigmae.AddExpression(*s);
+  sigmae_row.push_back(*s);
   delete s;
 }
 
 void
 ParsingDriver::add_to_row(ExpObj *v)
 {
-  mod_file->sigmae.AddExpression(get_expression(v));
+  sigmae_row.push_back(get_expression(v));
   delete v;
 }
 
 void
 ParsingDriver::steady()
 {
-  mod_file->computing_tasks.setSteady();
+  mod_file->addStatement(new SteadyStatement(options_list));
+  options_list.clear();
   mod_file->model_tree.computeJacobian = true;
 }
 
 void
 ParsingDriver::option_num(const string &name_option, string *opt1, string *opt2)
 {
-  mod_file->computing_tasks.setOption(name_option, *opt1, *opt2);
+  if (options_list.paired_num_options.find(name_option)
+      != options_list.paired_num_options.end())
+    error("option " + name_option + " declared twice");
+
+  options_list.paired_num_options[name_option] = pair<string, string>(*opt1, *opt2);
   delete opt1;
   delete opt2;
 }
@@ -484,11 +534,14 @@ ParsingDriver::option_num(const string &name_option, string *opt)
 void
 ParsingDriver::option_num(const string &name_option, const string &opt)
 {
-  mod_file->computing_tasks.setOption(name_option, opt);
+  if (options_list.num_options.find(name_option)
+      != options_list.num_options.end())
+    error("option " + name_option + " declared twice");
+
+  options_list.num_options[name_option] = opt;
+
   if (name_option == "order")
     mod_file->order = atoi(opt.c_str());
-  else if (name_option == "linear")
-    mod_file->linear = atoi(opt.c_str());
 }
 
 void
@@ -501,13 +554,25 @@ ParsingDriver::option_str(const string &name_option, string *opt)
 void
 ParsingDriver::option_str(const string &name_option, const string &opt)
 {
-  mod_file->computing_tasks.setOption(name_option, "'" + opt + "'");
+  if (options_list.string_options.find(name_option)
+      != options_list.string_options.end())
+    error("option " + name_option + " declared twice");
+
+  options_list.string_options[name_option] = opt;
+}
+
+void
+ParsingDriver::linear()
+{
+  mod_file->linear = 1;
 }
 
 void
 ParsingDriver::add_tmp_var(string *tmp_var1, string *tmp_var2)
 {
-  tmp_symbol_table.AddTempSymbol(*tmp_var1, *tmp_var2);
+  check_symbol_existence(*tmp_var1);
+  check_symbol_existence(*tmp_var2);
+  tmp_symbol_table->AddTempSymbol(*tmp_var1, *tmp_var2);
   delete tmp_var1;
   delete tmp_var2;
 }
@@ -515,15 +580,16 @@ ParsingDriver::add_tmp_var(string *tmp_var1, string *tmp_var2)
 void
 ParsingDriver::add_tmp_var(string *tmp_var)
 {
-  tmp_symbol_table.AddTempSymbol(*tmp_var);
+  check_symbol_existence(*tmp_var);
+  tmp_symbol_table->AddTempSymbol(*tmp_var);
   delete tmp_var;
 }
 
 void ParsingDriver::rplot()
 {
-  tmp_symbol_table.set("var_list_");
-  string tmp = tmp_symbol_table.get();
-  mod_file->computing_tasks.runRplot(tmp);
+  mod_file->addStatement(new RplotStatement(*tmp_symbol_table, options_list));
+  options_list.clear();
+  tmp_symbol_table->clear();
 }
 
 void ParsingDriver::stoch_simul()
@@ -535,187 +601,257 @@ void ParsingDriver::stoch_simul()
   if (mod_file->linear == -1)
     mod_file->linear = 0;
 
-  tmp_symbol_table.set("var_list_");
-  string tmp = tmp_symbol_table.get();
-  mod_file->computing_tasks.setStochSimul(tmp);
+  mod_file->addStatement(new StochSimulStatement(*tmp_symbol_table, options_list));
+  tmp_symbol_table->clear();
+  options_list.clear();
 }
 
 void
 ParsingDriver::simul()
 {
-  mod_file->computing_tasks.setSimul();
+  mod_file->addStatement(new SimulStatement(options_list));
+  options_list.clear();
   mod_file->model_tree.computeJacobian = true;
 }
 
 void
 ParsingDriver::check()
 {
-  mod_file->computing_tasks.setCheck();
+  mod_file->addStatement(new CheckStatement(options_list));
+  options_list.clear();
   mod_file->model_tree.computeJacobian = true;
 }
 
 void
-ParsingDriver::estimation_init()
+ParsingDriver::add_estimated_params_element()
 {
-  mod_file->computing_tasks.EstimParams = &estim_params;
-  mod_file->computing_tasks.setEstimationInit();
-  mod_file->model_tree.computeJacobianExo = true;
+  check_symbol_existence(estim_params.name);
+  if (estim_params.name2.size() > 0)
+    check_symbol_existence(estim_params.name2);
+
+  estim_params_list.push_back(estim_params);
+  estim_params.clear();
 }
 
 void
-ParsingDriver::set_estimated_elements()
+ParsingDriver::estimated_params()
 {
-  mod_file->computing_tasks.setEstimatedElements();
+  mod_file->addStatement(new EstimatedParamsStatement(estim_params_list, mod_file->symbol_table));
+  estim_params_list.clear();
+  mod_file->model_tree.computeJacobianExo = true;
 }
 
 void
-ParsingDriver::set_estimated_init_elements()
+ParsingDriver::estimated_params_init()
 {
-  mod_file->computing_tasks.setEstimatedInitElements();
+  mod_file->addStatement(new EstimatedParamsInitStatement(estim_params_list, mod_file->symbol_table));
+  estim_params_list.clear();
 }
 
 void
-ParsingDriver::set_estimated_bounds_elements()
+ParsingDriver::estimated_params_bounds()
 {
-  mod_file->computing_tasks.setEstimatedBoundsElements();
+  mod_file->addStatement(new EstimatedParamsBoundsStatement(estim_params_list, mod_file->symbol_table));
+  estim_params_list.clear();
 }
 
 void
 ParsingDriver::set_unit_root_vars()
 {
-  tmp_symbol_table.set("options_.unit_root_vars");
-  *output << tmp_symbol_table.get();
+  mod_file->addStatement(new UnitRootVarsStatement(*tmp_symbol_table));
+  tmp_symbol_table->clear();
 }
 
 void
 ParsingDriver::run_estimation()
 {
-  tmp_symbol_table.set("var_list_");
-  string tmp = tmp_symbol_table.get();
-  mod_file->computing_tasks.runEstimation(tmp);
+  mod_file->addStatement(new EstimationStatement(*tmp_symbol_table, options_list));
+  tmp_symbol_table->clear();
+  options_list.clear();
 }
 
 void
 ParsingDriver::run_prior_analysis()
 {
-  tmp_symbol_table.set("var_list_");
-  string tmp = tmp_symbol_table.get();
-  mod_file->computing_tasks.runPriorAnalysis(tmp);
+  mod_file->addStatement(new PriorAnalysisStatement(*tmp_symbol_table, options_list));
+  tmp_symbol_table->clear();
+  options_list.clear();
 }
 
 void
 ParsingDriver::run_posterior_analysis()
 {
-  tmp_symbol_table.set("var_list_");
-  string tmp = tmp_symbol_table.get();
-  mod_file->computing_tasks.runPosteriorAnalysis(tmp);
+  mod_file->addStatement(new PosteriorAnalysisStatement(*tmp_symbol_table, options_list));
+  tmp_symbol_table->clear();
+  options_list.clear();
+}
+
+void
+ParsingDriver::optim_options_helper(const string &name)
+{
+  if (options_list.string_options.find("optim_opt") == options_list.string_options.end())
+    options_list.string_options["optim_opt"] = "";
+  else
+    options_list.string_options["optim_opt"] += ",";
+  options_list.string_options["optim_opt"] += "''" + name + "'',";
+}
+
+void
+ParsingDriver::optim_options_string(string *name, string *value)
+{
+  optim_options_helper(*name);
+  options_list.string_options["optim_opt"] += "''" + *value + "''";
+  delete name;
+  delete value;
 }
 
 void
-ParsingDriver::optim_options(string *str1, string *str2, int task)
+ParsingDriver::optim_options_num(string *name, string *value)
 {
-  mod_file->computing_tasks.setOptimOptions(*str1, *str2, task);
-  delete str1;
-  delete str2;
+  optim_options_helper(*name);
+  options_list.string_options["optim_opt"] += *value;
+  delete name;
+  delete value;
 }
 
 void
 ParsingDriver::set_varobs()
 {
-  tmp_symbol_table.set("options_.varobs");
-  *output << tmp_symbol_table.get();
+  mod_file->addStatement(new VarobsStatement(*tmp_symbol_table));
+  tmp_symbol_table->clear();
 }
 
 void
-ParsingDriver::set_trend_init()
+ParsingDriver::set_trends()
 {
-  *output << "options_.trend_coeff_ = {};" << endl;
+  mod_file->addStatement(new ObservationTrendsStatement(trend_elements, mod_file->symbol_table));
+  trend_elements.clear();
 }
 
 void
 ParsingDriver::set_trend_element(string *arg1, ExpObj *arg2)
 {
-  mod_file->computing_tasks.set_trend_element(*arg1, get_expression(arg2));
+  check_symbol_existence(*arg1);
+  if (trend_elements.find(*arg1) != trend_elements.end())
+    error("observation_trends: " + *arg1 + " declared twice");
+  trend_elements[*arg1] = get_expression(arg2);
   delete arg1;
   delete arg2;
 }
 
 void
-ParsingDriver::begin_optim_weights()
+ParsingDriver::set_optim_weights(string *name, ExpObj *value)
 {
-  mod_file->computing_tasks.BeginOptimWeights();
+  check_symbol_existence(*name);
+  if (mod_file->symbol_table.getType(*name) != eEndogenous)
+    error("optim_weights: " + *name + " isn't an endogenous variable");
+  if (var_weights.find(*name) != var_weights.end())
+    error("optim_weights: " + *name + " declared twice");
+  var_weights[*name] = get_expression(value);
+  delete name;
+  delete value;
 }
 
 void
-ParsingDriver::set_optim_weights(string *arg1, ExpObj *arg2)
+ParsingDriver::set_optim_weights(string *name1, string *name2, ExpObj *value)
 {
-  mod_file->computing_tasks.setOptimWeights(*arg1, get_expression(arg2));
-  delete arg1;
-  delete arg2;
+  check_symbol_existence(*name1);
+  if (mod_file->symbol_table.getType(*name1) != eEndogenous)
+    error("optim_weights: " + *name1 + " isn't an endogenous variable");
+
+  check_symbol_existence(*name2);
+  if (mod_file->symbol_table.getType(*name2) != eEndogenous)
+    error("optim_weights: " + *name2 + " isn't an endogenous variable");
+
+  pair<string, string> covar_key(*name1, *name2);
+
+  if (covar_weights.find(covar_key) != covar_weights.end())
+    error("optim_weights: pair of variables (" + *name1 + ", " + *name2
+          + ") declared twice");
+
+  covar_weights[covar_key] = get_expression(value);
+  delete name1;
+  delete name2;
+  delete value;
 }
 
 void
-ParsingDriver::set_optim_weights(string *arg1, string *arg2, ExpObj *arg3)
+ParsingDriver::optim_weights()
 {
-  mod_file->computing_tasks.setOptimWeights(*arg1, *arg2, get_expression(arg3));
-  delete arg1;
-  delete arg2;
-  delete arg3;
+  mod_file->addStatement(new OptimWeightsStatement(var_weights, covar_weights, mod_file->symbol_table));
+  var_weights.clear();
+  covar_weights.clear();
 }
 
 void
 ParsingDriver::set_osr_params()
 {
-  tmp_symbol_table.set("osr_params_");
-  string tmp = tmp_symbol_table.get();
-  mod_file->computing_tasks.setOsrParams(tmp);
+  mod_file->addStatement(new OsrParamsStatement(*tmp_symbol_table));
+  tmp_symbol_table->clear();
 }
 
 void
 ParsingDriver::run_osr()
 {
-  tmp_symbol_table.set("var_list_");
-  string tmp = tmp_symbol_table.get();
-  mod_file->computing_tasks.runOsr(tmp);
+  mod_file->addStatement(new OsrStatement(*tmp_symbol_table, options_list));
+  tmp_symbol_table->clear();
+  options_list.clear();
   mod_file->model_tree.computeJacobianExo = true;
 }
 
 void
 ParsingDriver::set_olr_inst()
 {
-  tmp_symbol_table.set("options_.olr_inst");
-  string tmp = tmp_symbol_table.get();
-  mod_file->computing_tasks.setOlrInst(tmp);
+  mod_file->addStatement(new OlrInstStatement(*tmp_symbol_table));
+  tmp_symbol_table->clear();
 }
 
 void
 ParsingDriver::run_olr()
 {
-  tmp_symbol_table.set("var_list_");
-  string tmp = tmp_symbol_table.get();
-  mod_file->computing_tasks.runOlr(tmp);
-}
-
-void
-ParsingDriver::begin_calib_var()
-{
-  mod_file->computing_tasks.BeginCalibVar();
+  mod_file->addStatement(new OlrStatement(*tmp_symbol_table, options_list));
+  tmp_symbol_table->clear();
+  options_list.clear();
 }
 
 void
 ParsingDriver::set_calib_var(string *name, string *weight, ExpObj *expression)
 {
-  mod_file->computing_tasks.setCalibVar(*name, *weight, get_expression(expression));
+  check_symbol_existence(*name);
+  if (mod_file->symbol_table.getType(*name) != eEndogenous
+      && mod_file->symbol_table.getType(*name) != eExogenous)
+    error("calib_var: " + *name + " isn't an endogenous or exogenous variable");
+
+  if (calib_var.find(*name) != calib_var.end())
+    error("calib_var: " + *name + " declared twice");
+
+  calib_var[*name] = pair<string, string>(*weight, get_expression(expression));
+
   delete name;
   delete weight;
   delete expression;
 }
 
 void
-ParsingDriver::set_calib_var(string *name1, string *name2,
-                             string *weight, ExpObj *expression)
+ParsingDriver::set_calib_covar(string *name1, string *name2,
+                               string *weight, ExpObj *expression)
 {
-  mod_file->computing_tasks.setCalibVar(*name1, *name2, *weight, get_expression(expression));
+  check_symbol_existence(*name1);
+  check_symbol_existence(*name2);
+  if (mod_file->symbol_table.getType(*name1) != mod_file->symbol_table.getType(*name2))
+    error("calib_var: " + *name1 + " and " + *name2 + "dont't have the same type");
+  if (mod_file->symbol_table.getType(*name1) != eEndogenous
+      && mod_file->symbol_table.getType(*name1) != eExogenous)
+    error("calib_var: " + *name1 + " and " + *name2 + "aren't endogenous or exogenous variables");
+
+  pair<string, string> covar_key(*name1, *name2);
+
+  if (calib_covar.find(covar_key) != calib_covar.end())
+    error("calib_var: pair of variables (" + *name1 + ", " + *name2
+          + ") declared twice");
+
+  calib_covar[covar_key] = pair<string, string>(*weight, get_expression(expression));
+
   delete name1;
   delete name2;
   delete weight;
@@ -726,7 +862,18 @@ void
 ParsingDriver::set_calib_ac(string *name, string *ar,
                             string *weight, ExpObj *expression)
 {
-  mod_file->computing_tasks.setCalibAc(*name, *ar, *weight, get_expression(expression));
+  check_symbol_existence(*name);
+  if (mod_file->symbol_table.getType(*name) != eEndogenous)
+    error("calib_var: " + *name + "isn't an endogenous variable");
+
+  int iar = atoi(ar->c_str());
+  pair<string, int> ac_key(*name, iar);
+
+  if (calib_ac.find(ac_key) != calib_ac.end())
+    error("calib_var: autocorr " + *name + "(" + *ar + ") declared twice");
+
+  calib_ac[ac_key] = pair<string, string>(*weight, get_expression(expression));
+
   delete name;
   delete ar;
   delete weight;
@@ -734,41 +881,45 @@ ParsingDriver::set_calib_ac(string *name, string *ar,
 }
 
 void
-ParsingDriver::run_calib(int flag)
+ParsingDriver::run_calib_var()
 {
-  mod_file->computing_tasks.runCalib(flag);
+  mod_file->addStatement(new CalibVarStatement(calib_var, calib_covar, calib_ac,
+                                               mod_file->symbol_table));
+  calib_var.clear();
+  calib_covar.clear();
+  calib_ac.clear();
 }
 
 void
-ParsingDriver::run_dynatype(string *filename, string *ext)
+ParsingDriver::run_calib(int covar)
 {
-  tmp_symbol_table.set("var_list_");
-  string tmp = tmp_symbol_table.get();
-  mod_file->computing_tasks.runDynatype(*filename, *ext, tmp);
-  delete filename;
-  delete ext;
+  mod_file->addStatement(new CalibStatement(covar));
 }
 
 void
-ParsingDriver::run_dynasave(string *filename, string *ext)
+ParsingDriver::run_dynatype(string *filename, string *ext)
 {
-  tmp_symbol_table.set("var_list_");
-  string tmp = tmp_symbol_table.get();
-  mod_file->computing_tasks.runDynasave(*filename, *ext, tmp);
+  mod_file->addStatement(new DynaTypeStatement(*tmp_symbol_table, *filename, *ext));
+  tmp_symbol_table->clear();
   delete filename;
   delete ext;
 }
 
 void
-ParsingDriver::begin_model_comparison()
+ParsingDriver::run_dynasave(string *filename, string *ext)
 {
-  mod_file->computing_tasks.beginModelComparison();
+  mod_file->addStatement(new DynaSaveStatement(*tmp_symbol_table, *filename, *ext));
+  tmp_symbol_table->clear();
+  delete filename;
+  delete ext;
 }
 
 void
 ParsingDriver::add_mc_filename(string *filename, string *prior)
 {
-  mod_file->computing_tasks.addMcFilename(*filename, *prior);
+  if (filename_list.find(*filename) != filename_list.end())
+    error("model_comparison: filename " + *filename + " declared twice");
+  filename_list[*filename] = *prior;
   delete filename;
   delete prior;
 }
@@ -776,7 +927,9 @@ ParsingDriver::add_mc_filename(string *filename, string *prior)
 void
 ParsingDriver::run_model_comparison()
 {
-  mod_file->computing_tasks.runModelComparison();
+  mod_file->addStatement(new ModelComparisonStatement(filename_list, options_list));
+  filename_list.clear();
+  options_list.clear();
 }
 
 NodeID
@@ -937,5 +1090,5 @@ ParsingDriver::add_model_sqrt(NodeID arg1)
 void
 ParsingDriver::add_native(const char *s)
 {
-  *output << s;
+  mod_file->addStatement(new NativeStatement(s));
 }
diff --git a/parser.src/Shocks.cc b/parser.src/Shocks.cc
index fa12940ea53088ce09026ae961aac0a5b7e019d9..7b564fbd6201accd3e92811314d6c38c9935e4dc 100644
--- a/parser.src/Shocks.cc
+++ b/parser.src/Shocks.cc
@@ -11,185 +11,152 @@ using namespace std;
 #include "ModelParameters.hh"
 #include "Interface.hh"
 
-Shocks::Shocks() : mshock_flag(0), exo_det_length(0)
+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)
 {
-  // Empty
 }
 
-//------------------------------------------------------------------------------
-Shocks::~Shocks()
-{
-  // Empty
-}
-
-//------------------------------------------------------------------------------
-void Shocks::setOutput(ostringstream* iOutput)
+void
+AbstractShocksStatement::writeDetShocks(ostream &output) const
 {
-  output = iOutput;
-}
+  int exo_det_length = 0;
 
-//------------------------------------------------------------------------------
-void Shocks::BeginShocks(void)
-{
-  mshock_flag = 0;
-  // Writing a Matlab comment
-  *output << interfaces::comment() << "\n" << interfaces::comment() << "SHOCKS instructions \n"
-          << interfaces::comment() << "\n";
-  //Writing intstruction that initialize a shocks
-  *output << "make_ex_;\n";
-}
+  for(det_shocks_type::const_iterator it = det_shocks.begin();
+      it != det_shocks.end(); it++)
+    {
+      int id = symbol_table.getID(it->first) + 1;
+      bool exo_det = (symbol_table.getType(it->first) == eExogenousDet);
+      int set_shocks_index = ((int) mshocks) + 2 * ((int) exo_det);
 
-//------------------------------------------------------------------------------
-void Shocks::BeginMShocks(void)
-{
-  mshock_flag = 1;
-  // Writing a Matlab comment
-  *output << interfaces::comment() << "\n" << interfaces::comment() << "MSHOCKS instructions \n"
-          << interfaces::comment() << "\n";
-  //Writing intstruction that initialize a shocks
-  *output << "make_ex_;\n";
+      for (unsigned int i = 0; i < it->second.size(); i++)
+        {
+          const int &period1 = it->second[i].period1;
+          const int &period2 = it->second[i].period2;
+          const string &value = it->second[i].value;
+
+          if (period1 == period2)
+            output << "set_shocks(" << set_shocks_index << "," << period1
+                   << ", " << id << ", " << value << ");\n";
+          else
+            output << "set_shocks(" << set_shocks_index << "," << period1
+                   << ":" << period2 << ", " << id
+                   << ", " << value << ");\n";
+
+          if (exo_det && (period2 > exo_det_length))
+            exo_det_length = period2;
+        }
+    }
+  output << "M_.exo_det_length = " << exo_det_length << ";\n";
 }
 
-void Shocks::EndShocks(void)
+void
+AbstractShocksStatement::writeVarAndStdShocks(ostream &output) const
 {
-  *output << "M_.exo_det_length = " << exo_det_length << ";\n";
-}
+  var_and_std_shocks_type::const_iterator it;
 
-//------------------------------------------------------------------------------
-void Shocks::AddDetShockExo(int id1)
-{
-  if (mPeriod1.size() != mPeriod2.size() ||
-      mPeriod1.size() != mValues.size())
+  for(it = var_shocks.begin(); it != var_shocks.end(); it++)
     {
-      string msg = "in shocks statement, number of periods and values dont agree";
-      (* error) (msg.c_str());
+      int id = symbol_table.getID(it->first) + 1;
+      const string &value = it->second;
+      output << "M_.Sigma_e(" << id << ", " << id << ") = " << value << ";\n";
     }
-  for (unsigned int i = 0; i < mPeriod1.size(); i++)
+
+  for(it = std_shocks.begin(); it != std_shocks.end(); it++)
     {
-      string period1 = mPeriod1[i];
-      string period2 = mPeriod2[i];
-      if (period1 == period2)
-        {
-          *output << "set_shocks(" << mshock_flag << "," << period1
-                  << ", " << id1+1 << ", " << mValues[i]
-                  << ");\n";
-        }
-      else
-        {
-          *output << "set_shocks(" << mshock_flag << "," << period1
-                  << ":" << period2 << ", " << id1+1
-                  << ", " << mValues[i] << ");\n";
-        }
+      int id = symbol_table.getID(it->first) + 1;
+      const string &value = it->second;
+      output << "M_.Sigma_e(" << id << ", " << id << ") = " << value << "^2;\n";
     }
-  mPeriod1.clear();
-  mPeriod2.clear();
-  mValues.clear();
 }
 
-void Shocks::AddDetShockExoDet(int id1)
+void
+AbstractShocksStatement::writeCovarAndCorrShocks(ostream &output) const
 {
-  if (mPeriod1.size() != mPeriod2.size() ||
-      mPeriod1.size() != mValues.size())
+  covar_and_corr_shocks_type::const_iterator it;
+
+  for(it = covar_shocks.begin(); it != covar_shocks.end(); it++)
     {
-      string msg = "in shocks statement, number of periods and values dont agree";
-      (* error) (msg.c_str());
+      int id1 = symbol_table.getID(it->first.first) + 1;
+      int id2 = symbol_table.getID(it->first.second) + 1;
+      const string &value = it->second;
+      output << "M_.Sigma_e(" << id1 << ", " << id2 << ") = " << value
+             << "; M_.Sigma_e(" << id2 << ", " << id1 << ") = M_.Sigma_e("
+             << id1 << ", " << id2 << ");\n";
     }
-  for (unsigned int i = 0; i < mPeriod1.size(); i++)
+
+  for(it = corr_shocks.begin(); it != corr_shocks.end(); it++)
     {
-      string period1 = mPeriod1[i];
-      string period2 = mPeriod2[i];
-      if (period1 == period2)
-        {
-          *output << "set_shocks(" << mshock_flag + 2 << "," << period1
-                  << ", " << id1+1 << ", " << mValues[i]
-                  << ");\n";
-        }
-      else
-        {
-          *output << "set_shocks(" << mshock_flag+2 << "," << period1
-                  << ":" << period2 << ", " << id1+1
-                  << ", " << mValues[i] << ");\n";
-        }
-      int p2_int = atoi(period2.c_str());
-      if (p2_int > exo_det_length)
-        {
-          exo_det_length = p2_int;
-        }
+      int id1 = symbol_table.getID(it->first.first) + 1;
+      int id2 = symbol_table.getID(it->first.second) + 1;
+      const string &value = it->second;
+      output << "M_.Sigma_e(" << id1 << ", " << id2 << ") = " << value
+             << "*sqrt(M_.Sigma_e(" << id1 << ", " << id1 << ")*M_.Sigma_e("
+             << id2 << ", " << id2 << "); M_.Sigma_e(" << id2 << ", "
+             << id1 << ") = M_.Sigma_e(" << id1 << ", " << id2 << ");\n";
     }
-  mPeriod1.clear();
-  mPeriod2.clear();
-  mValues.clear();
 }
 
-void Shocks::AddSTDShock(int id1, std::string value)
-{
-  *output << "M_.Sigma_e(" << id1+1 << ", " << id1+1 << ") = " << value << "^2;\n";
-  mPeriod1.clear();
-  mPeriod2.clear();
-  mValues.clear();
-}
 
-void Shocks::AddVARShock(int id1, std::string value)
+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)
 {
-  *output << "M_.Sigma_e(" << id1+1 << ", " << id1+1 << ") = " << value << ";\n";
-  mPeriod1.clear();
-  mPeriod2.clear();
-  mValues.clear();
 }
 
-void Shocks::AddCOVAShock(int id1, int id2 , std::string value)
+void
+ShocksStatement::writeOutput(ostream &output) const
 {
-  *output << "M_.Sigma_e(" << id1+1 << ", " << id2+1 << ") = " << value <<
-    "; M_.Sigma_e(" << id2+1 << ", " << id1+1 << ") = M_.Sigma_e(" <<
-    id1+1 << ", " << id2+1 << ");\n";
-  mPeriod1.clear();
-  mPeriod2.clear();
-  mValues.clear();
-}
+  output << interfaces::comment() << endl << interfaces::comment();
+  output << "SHOCKS instructions \n";
+  output << interfaces::comment() << "\n";
 
-void Shocks::AddCORRShock(int id1, int id2 , std::string value)
-{
-  *output << "M_.Sigma_e(" << id1+1 << ", " << id2+1 << ") = " <<
-    value << "*sqrt(M_.Sigma_e(" << id1 << ", " << id1+1 << ")*M_.Sigma_e(" <<
-    id2 << ", " << id2+1 << "); M_.Sigma_e(" << id2+1 << ", " <<
-    id1 << ") = M_.Sigma_e(" << id1+1 << ", " << id2+1 << ");\n";
-  mPeriod1.clear();
-  mPeriod2.clear();
-  mValues.clear();
+  // Write instruction that initializes a shock
+  output << "make_ex_;\n";
+
+  writeDetShocks(output);
+  writeVarAndStdShocks(output);
+  writeCovarAndCorrShocks(output);
 }
 
-/*
-  sets a shock
-  for type == 0, set type, id1, shock_elem will be set by another method
-  for type == 1 or 2 set type, id1 and value
-  for type == 3 or 4 set type, id1, id2 and value
-*/
-//------------------------------------------------------------------------------
-/*
-  string Shocks::get(void)
-  {
-  return output.str();
-  }
-*/
-//------------------------------------------------------------------------------
-void Shocks::AddPeriod(string p1, string p2)
+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)
 {
-  mPeriod1.push_back(p1);
-  mPeriod2.push_back(p2);
 }
 
-//------------------------------------------------------------------------------
-void Shocks::AddPeriod(string p1)
+void
+MShocksStatement::writeOutput(ostream &output) const
 {
-  mPeriod1.push_back(p1);
-  mPeriod2.push_back(p1);
+  output << interfaces::comment() << endl << interfaces::comment();
+  output << "SHOCKS instructions \n";
+  output << interfaces::comment() << "\n";
 
-}
+  // Write instruction that initializes a shock
+  output << "make_ex_;\n";
 
-//------------------------------------------------------------------------------
-void Shocks::AddValue(string value)
-{
-  mValues.push_back(value);
+  writeDetShocks(output);
+  writeVarAndStdShocks(output);
+  writeCovarAndCorrShocks(output);
 }
-
-//------------------------------------------------------------------------------
diff --git a/parser.src/SigmaeInitialization.cc b/parser.src/SigmaeInitialization.cc
index 71e4599f57d963d57e18f7e0617646790dd41c7d..a87510751746799028b05d9fd35269696978b2a0 100644
--- a/parser.src/SigmaeInitialization.cc
+++ b/parser.src/SigmaeInitialization.cc
@@ -1,52 +1,17 @@
-/*! \file
-  \version 1.0
-  \date 04/13/2004
-  \par This file defines the SigmaeInitialization class.
-*/
-//------------------------------------------------------------------------------
-#include <iostream>
-
-using namespace std;
-
 #include "SigmaeInitialization.hh"
-//------------------------------------------------------------------------------
-SigmaeInitialization::SigmaeInitialization()
-{
-  // Empty
-}
-
-//------------------------------------------------------------------------------
-SigmaeInitialization::~SigmaeInitialization()
-{
-  // Empty
-}
-
-//------------------------------------------------------------------------------
-void SigmaeInitialization::setOutput(ostringstream* iOutput)
-{
-  output = iOutput;
-}
 
-//------------------------------------------------------------------------------
-void SigmaeInitialization::AddExpression(string expression)
+SigmaeStatement::SigmaeStatement(const matrix_type &matrix_arg) throw (MatrixFormException) :
+  matrix(matrix_arg),
+  matrix_form(determineMatrixForm(matrix))
 {
-  row.push_back(expression);
 }
 
-//------------------------------------------------------------------------------
-void SigmaeInitialization::EndOfRow()
+SigmaeStatement::matrix_form_type
+SigmaeStatement::determineMatrixForm(const matrix_type &matrix) throw (MatrixFormException)
 {
-
-  matrix.push_back(row);
-  row.clear();
-}
-
-//------------------------------------------------------------------------------
-void SigmaeInitialization::CheckMatrix(void)
-{
-  vector<vector<string> >::iterator ir;
   unsigned int nbe;
   int inc;
+  matrix_form_type type;
   // Checking if first or last row has one element.
   if (matrix.front().size() == 1)
     {
@@ -61,69 +26,55 @@ void SigmaeInitialization::CheckMatrix(void)
       type = eUpper;
     }
   else
-    {
-      string msg = "sigma_e isn't in triangular form";
-      (* error) (msg.c_str());
-    }
+    throw MatrixFormException();
+
   // Checking if matrix is triangular (upper or lower):
   // each row has one element more or less than the previous one
   // and first or last one has one element.
+  matrix_type::const_iterator ir;
   for (ir = matrix.begin(), ir++; ir != matrix.end(); ir++, nbe += inc )
     if (ir->size() != nbe)
-      {
-        string msg = "sigma_e isn't in triangular form!";
-        (* error) (msg.c_str());
-      }
+      throw MatrixFormException();
+
+  return type;
 }
 
-//------------------------------------------------------------------------------
-void SigmaeInitialization::SetMatrix(void)
+void
+SigmaeStatement::writeOutput(ostream &output) const
 {
   unsigned int ic, ic1;
   unsigned int ir, ir1;
 
-  *output << "M_.Sigma_e = [...\n";
+  output << "M_.Sigma_e = [...\n";
   for (ir = 0; ir < matrix.size(); ir++)
     {
-      *output << "\t";
+      output << "\t";
       for (ic = 0; ic < matrix.size(); ic++)
         {
-          if (ic >= ir && type == eUpper)
+          if (ic >= ir && matrix_form == eUpper)
             {
               ic1 = ic-ir;
               ir1 = ir;
             }
-          else if (ic < ir && type == eUpper)
+          else if (ic < ir && matrix_form == eUpper)
             {
               ic1 = ir-ic;
               ir1 = ic;
             }
-          else if (ic > ir && type == eLower)
+          else if (ic > ir && matrix_form == eLower)
             {
               ic1 = ir;
               ir1 = ic;
             }
-          else if (ic <= ir && type == eLower)
+          else if (ic <= ir && matrix_form == eLower)
             {
               ic1 = ic;
               ir1 = ir;
             }
 
-          *output << matrix[ir1][ic1] << " ";
+          output << matrix[ir1][ic1] << " ";
         }
-      *output << ";...\n";
+      output << ";...\n";
     }
-  *output << "\t];...\n";
+  output << "\t];...\n";
 }
-
-//------------------------------------------------------------------------------
-void SigmaeInitialization::set(void)
-{
-  CheckMatrix();
-
-  SetMatrix();
-  matrix.clear();
-  row.clear();
-}
-
-//------------------------------------------------------------------------------
diff --git a/parser.src/Statement.cc b/parser.src/Statement.cc
new file mode 100644
index 0000000000000000000000000000000000000000..540ca93d42e85d0633d72a947ea6bb6e5cdfc0e1
--- /dev/null
+++ b/parser.src/Statement.cc
@@ -0,0 +1,41 @@
+#include "Statement.hh"
+
+Statement::~Statement()
+{
+}
+
+NativeStatement::NativeStatement(const string &native_statement_arg) :
+  native_statement(native_statement_arg)
+{
+}
+
+void
+NativeStatement::writeOutput(ostream &output) const
+{
+  output << native_statement << endl;
+}
+
+void
+OptionsList::writeOutput(ostream &output) const
+{
+  for(num_options_type::const_iterator it = num_options.begin();
+      it != num_options.end(); it++)
+    output << "options_." << it->first << " = " << it->second << ";" << endl;
+
+  for(paired_num_options_type::const_iterator it = paired_num_options.begin();
+      it != paired_num_options.end(); it++)
+    output << "options_." << it->first << " = [" << it->second.first << "; "
+           << it->second.second << "];" << endl;
+
+  for(string_options_type::const_iterator it = string_options.begin();
+      it != string_options.end(); it++)
+    output << "options_." << it->first << " = '" << it->second << "';" << endl;
+}
+
+void
+OptionsList::clear()
+{
+  num_options.clear();
+  paired_num_options.clear();
+  string_options.clear();
+}
diff --git a/parser.src/SymbolTable.cc b/parser.src/SymbolTable.cc
index 66fcdedc5b3f7c01ff85bfd8151c13bb3fc3232c..306ee9c1131498d58f7efe817b37e05235bc48dc 100644
--- a/parser.src/SymbolTable.cc
+++ b/parser.src/SymbolTable.cc
@@ -189,10 +189,9 @@ void SymbolTable::clean()
 }
 #endif // Comment
 
-string SymbolTable::get()
+void
+SymbolTable::writeOutput(ostream &output)
 {
-  ostringstream output;
-
   if (mod_param.exo_nbr > 0)
     {
       output << "M_.exo_names = '" << getNameByID(eExogenous, 0) << "';\n";
@@ -251,5 +250,4 @@ string SymbolTable::get()
   output << "M_.endo_nbr = " << mod_param.endo_nbr << ";\n";
   output << "M_.recur_nbr = " << mod_param.recur_nbr << ";\n";
   output << "M_.param_nbr = " << mod_param.parameter_nbr << ";\n";
-  return output.str();
 }
diff --git a/parser.src/TmpSymbolTable.cc b/parser.src/TmpSymbolTable.cc
index 2d1ecc9dd4fbe16760484e3fe72730eb552d2603..449bb27f205021c53b7db060c53c8bad1ab4b555 100644
--- a/parser.src/TmpSymbolTable.cc
+++ b/parser.src/TmpSymbolTable.cc
@@ -10,7 +10,8 @@ using namespace std;
 #include "TmpSymbolTable.hh"
 #include "Interface.hh"
 //------------------------------------------------------------------------------
-TmpSymbolTable::TmpSymbolTable()
+TmpSymbolTable::TmpSymbolTable(const SymbolTable &symbol_table_arg) :
+  symbol_table(symbol_table_arg)
 {
   // Empty
 }
@@ -20,70 +21,39 @@ TmpSymbolTable::~TmpSymbolTable()
   // Empty
 }
 
-void TmpSymbolTable::setGlobalSymbolTable(SymbolTable *symbol_table_arg)
+void
+TmpSymbolTable::AddTempSymbol(const string &symbol)
 {
-  symbol_table = symbol_table_arg;
+  // FIXME: add check to verify that symbol exists in symbol_table
+  // FIXME: add check to verify that symbol doesn't yet exist in the present table
+  tmpsymboltable.push_back(symbol);
 }
 
-void TmpSymbolTable::AddTempSymbol(string symbol)
+void
+TmpSymbolTable::AddTempSymbol(const string &symbol1, const string &symbol2)
 {
-  if (symbol_table->Exist(symbol))
-    tmpsymboltable.push_back(symbol);
-  else
-    {
-      string msg = "Unknown symbol : "+symbol;
-      error(msg.c_str());
-    }
+  // FIXME: add checks to verify that symbol1 and symbol2 exist in symbol_table
+  // FIXME: add check to verify that symbol1 doesn't yet exist in the present table
+  tmpsymboltable.push_back(symbol1);
+  nameTable.push_back(symbol2);
 }
 
-//------------------------------------------------------------------------------
-void  TmpSymbolTable::AddTempSymbol(string symbol1, string symbol2)
+void
+TmpSymbolTable::writeOutput(const string &varname, ostream &output) const
 {
-  if (symbol_table->Exist(symbol1))
-    tmpsymboltable.push_back(symbol1);
-  else
-    {
-      string msg = "Unknown symbol : "+symbol1;
-      error(msg.c_str());
-    }
-
-  if (symbol_table->Exist(symbol2))
-    NameTable.push_back(symbol2);
-  else
-    {
-      string msg = "Unknown symbol : "+symbol2;
-      error(msg.c_str());
-    }
+  output << varname << "=[];" << endl;
+  for (vector<string>::const_iterator it = tmpsymboltable.begin();
+       it != tmpsymboltable.end(); it++)
+    if (symbol_table.isReferenced(*it) == eReferenced)
+      {
+        output << varname << " = ";
+        output << interfaces::strvcat(varname, "'" + *it + "'") << ";" << endl;
+      }
 }
 
-//------------------------------------------------------------------------------
-void TmpSymbolTable::set(string varname)
-{
-  list<string>::iterator it;
-  output.str("");
-  output << "\n" << varname << "=[];\n";
-  for (it = tmpsymboltable.begin(); it!= tmpsymboltable.end(); it++)
-    {
-      if (symbol_table->isReferenced(*it) == eReferenced)
-        {
-          output << varname << " = ";
-          output << interfaces::strvcat(varname,"'"+*it+"'")+";\n";
-        }
-    }
-}
-
-//------------------------------------------------------------------------------
-string TmpSymbolTable::get(void)
+void
+TmpSymbolTable::clear()
 {
   tmpsymboltable.clear();
-  NameTable.clear();
-  return output.str();
+  nameTable.clear();
 }
-
-//------------------------------------------------------------------------------
-int TmpSymbolTable::size(void)
-{
-  return tmpsymboltable.size();
-}
-
-//------------------------------------------------------------------------------
diff --git a/parser.src/include/ComputingTasks.hh b/parser.src/include/ComputingTasks.hh
index a37826ccda7afa114ea200dd23d125bc8b85caf9..25d5b8fcadd2cb14a5798b0642da1340a6d8d856 100644
--- a/parser.src/include/ComputingTasks.hh
+++ b/parser.src/include/ComputingTasks.hh
@@ -1,17 +1,221 @@
 #ifndef _COMPUTINGTASKS_HH
 #define _COMPUTINGTASKS_HH
-//------------------------------------------------------------------------------
-/** \file
- * \version 1.0
- * \date 12/16/2003
- * \par This file defines the ComputingTasks class .
- */
-//------------------------------------------------------------------------------
-#include <sstream>
-//------------------------------------------------------------------------------
+
+#include <ostream>
+
 #include "TmpSymbolTable.hh"
 #include "SymbolTable.hh"
-//------------------------------------------------------------------------------
+#include "Statement.hh"
+
+class SteadyStatement : public Statement
+{
+private:
+  const OptionsList options_list;
+public:
+  SteadyStatement(const OptionsList &options_list_arg);
+  virtual void writeOutput(ostream &output) const;
+};
+
+class CheckStatement : public Statement
+{
+private:
+  const OptionsList options_list;
+public:
+  CheckStatement(const OptionsList &options_list_arg);
+  virtual void writeOutput(ostream &output) const;
+};
+
+class SimulStatement : public Statement
+{
+private:
+  const OptionsList options_list;
+public:
+  SimulStatement(const OptionsList &options_list_arg);
+  virtual void writeOutput(ostream &output) const;
+};
+
+class StochSimulStatement : public Statement
+{
+private:
+  const TmpSymbolTable tmp_symbol_table;
+  const OptionsList options_list;
+public:
+  StochSimulStatement(const TmpSymbolTable &tmp_symbol_table_arg,
+                      const OptionsList &options_list_arg);
+  virtual void writeOutput(ostream &output) const;
+};
+
+class RplotStatement : public Statement
+{
+private:
+  const TmpSymbolTable tmp_symbol_table;
+  const OptionsList options_list;
+public:
+  RplotStatement(const TmpSymbolTable &tmp_symbol_table_arg,
+                 const OptionsList &options_list_arg);
+  virtual void writeOutput(ostream &output) const;
+};
+
+class UnitRootVarsStatement : public Statement
+{
+private:
+  const TmpSymbolTable tmp_symbol_table;
+public:
+  UnitRootVarsStatement(const TmpSymbolTable &tmp_symbol_table_arg);
+  virtual void writeOutput(ostream &output) const;
+};
+
+class PeriodsStatement : public Statement
+{
+private:
+  const int periods;
+public:
+  PeriodsStatement(int periods_arg);
+  virtual void writeOutput(ostream &output) const;
+};
+
+class DsampleStatement : public Statement
+{
+private:
+  const int val1, val2;
+public:
+  DsampleStatement(int val1_arg);
+  DsampleStatement(int val1_arg, int val2_arg);
+  virtual void writeOutput(ostream &output) const;
+};
+
+class EstimationStatement : public Statement
+{
+private:
+  const TmpSymbolTable tmp_symbol_table;
+  const OptionsList options_list;
+public:
+  EstimationStatement(const TmpSymbolTable &tmp_symbol_table_arg,
+                      const OptionsList &options_list_arg);
+  virtual void writeOutput(ostream &output) const;
+};
+
+class PriorAnalysisStatement : public Statement
+{
+private:
+  const TmpSymbolTable tmp_symbol_table;
+  const OptionsList options_list;
+public:
+  PriorAnalysisStatement(const TmpSymbolTable &tmp_symbol_table_arg,
+                         const OptionsList &options_list_arg);
+  virtual void writeOutput(ostream &output) const;
+};
+
+class PosteriorAnalysisStatement : public Statement
+{
+private:
+  const TmpSymbolTable tmp_symbol_table;
+  const OptionsList options_list;
+public:
+  PosteriorAnalysisStatement(const TmpSymbolTable &tmp_symbol_table_arg,
+                             const OptionsList &options_list_arg);
+  virtual void writeOutput(ostream &output) const;
+};
+
+class VarobsStatement : public Statement
+{
+private:
+  const TmpSymbolTable tmp_symbol_table;
+public:
+  VarobsStatement(const TmpSymbolTable &tmp_symbol_table_arg);
+  virtual void writeOutput(ostream &output) const;
+};
+
+class ObservationTrendsStatement : public Statement
+{
+public:
+  typedef map<string, string, less<string> > trend_elements_type;
+private:
+  const trend_elements_type trend_elements;
+  const SymbolTable &symbol_table;
+public:
+  ObservationTrendsStatement(const trend_elements_type &trend_elements_arg,
+                             const SymbolTable &symbol_table_arg);
+  virtual void writeOutput(ostream &output) const;
+};
+
+class OsrParamsStatement : public Statement
+{
+private:
+  const TmpSymbolTable tmp_symbol_table;
+public:
+  OsrParamsStatement(const TmpSymbolTable &tmp_symbol_table_arg);
+  virtual void writeOutput(ostream &output) const;
+};
+
+class OsrStatement : public Statement
+{
+private:
+  const TmpSymbolTable tmp_symbol_table;
+  const OptionsList options_list;
+public:
+  OsrStatement(const TmpSymbolTable &tmp_symbol_table_arg,
+               const OptionsList &options_list_arg);
+  virtual void writeOutput(ostream &output) const;
+};
+
+class OlrStatement : public Statement
+{
+private:
+  const TmpSymbolTable tmp_symbol_table;
+  const OptionsList options_list;
+public:
+  OlrStatement(const TmpSymbolTable &tmp_symbol_table_arg,
+               const OptionsList &options_list_arg);
+  virtual void writeOutput(ostream &output) const;
+};
+
+class OlrInstStatement : public Statement
+{
+private:
+  const TmpSymbolTable tmp_symbol_table;
+public:
+  OlrInstStatement(const TmpSymbolTable &tmp_symbol_table_arg);
+  virtual void writeOutput(ostream &output) const;
+};
+
+class DynaTypeStatement : public Statement
+{
+private:
+  const TmpSymbolTable tmp_symbol_table;
+  const string filename;
+  const string ext;
+public:
+  DynaTypeStatement(const TmpSymbolTable &tmp_symbol_table_arg,
+                    const string &filename_arg, const string &ext_arg);
+  virtual void writeOutput(ostream &output) const;
+};
+
+class DynaSaveStatement : public Statement
+{
+private:
+  const TmpSymbolTable tmp_symbol_table;
+  const string filename;
+  const string ext;
+public:
+  DynaSaveStatement(const TmpSymbolTable &tmp_symbol_table_arg,
+                    const string &filename_arg, const string &ext_arg);
+  virtual void writeOutput(ostream &output) const;
+};
+
+class ModelComparisonStatement : public Statement
+{
+public:
+  typedef map<string, string, less<string> > filename_list_type;
+private:
+  filename_list_type filename_list;
+  OptionsList options_list;
+public:
+  ModelComparisonStatement(const filename_list_type &filename_list_arg,
+                           const OptionsList &options_list_arg);
+  virtual void writeOutput(ostream &output) const;
+};
+
 /*!
   \class EstimationParams
   \brief EstimationParams
@@ -52,101 +256,84 @@ struct EstimationParams
   }
 };
 
-/*!
-  \class  ComputingTasks
-  \brief  This class concerns following Dynare commands :
-  steady, check, simul, stoch_simul, dynare_estimation,  calib, osr and olr
-*/
-class ComputingTasks
+class EstimatedParamsStatement : public Statement
 {
-private :
-  //! Output of this class
-  std::ostringstream  *output;
-  //! A reference to the symbol table
+private:
+  const vector<EstimationParams> estim_params_list;
   const SymbolTable &symbol_table;
+public:
+  EstimatedParamsStatement(const vector<EstimationParams> &estim_params_list_arg,
+                           const SymbolTable &symbol_table_arg);
+  virtual void writeOutput(ostream &output) const;
+};
 
-public :
-  //! Constructor
-  ComputingTasks(const SymbolTable &symbol_table_arg);
-  //! Destructor
-  ~ComputingTasks();
-  /*! Pointer to error function of parser class */
-  void (* error) (const char* m);
-  /*! Estimation parameters */
-  EstimationParams  *EstimParams;
-  /*!
-    Sets output reference
-    \param iOutput : reference to an ostringstream
-  */
-  void  setOutput(std::ostringstream* iOutput);
-  // Prints "steady;" to output
-  void  setSteady(void);
-  // Prints "check;" to output
-  void  setCheck(void);
-  // Prints "simul(oo_.dr);" to output
-  void  setSimul(void);
-  //! Prints tmp1 and "stoch_simul(var_list_);"  to output
-  void  setStochSimul(std::string tmp1);
-  /*! Sets an option by printing
-    option_.iName = iValue
-  */
-  void  setOption(std::string iName, std::string iValue);
-  void  setOption(std::string iName, std::string iValue1, std::string iValue2);
-  /*! Prints "dynare_estimation;" */
-  void  runEstimation(std::string);
-  void  runPriorAnalysis(std::string);
-  void  runPosteriorAnalysis(std::string);
-  void    runRplot(std::string);
-  /*! Prints some estimation initialisation */
-  void  setEstimationInit(void);
-  //! Prints optimization options
-  void  setOptimOptions(std::string str1, std::string str2, int task);
-  /*! Prints estimated elements */
-  void  setEstimatedElements(void);
-  void  setEstimatedInitElements(void);
-  void  setEstimatedBoundsElements(void);
-  void  setEstimationStandardError(void);
-  void    set_trend_element(std::string, std::string);
-
-  void  BeginCalibVar(void);
-
-  void  setCalibVar(std::string, std::string, std::string);
-
-  void  setCalibVar(std::string, std::string, std::string, std::string);
-
-  void  setCalibAc(std::string, std::string, std::string, std::string);
-
-  void  runCalib(int);
-
-  //! Prints "osr(var_list_,osr_params_,optim_weights_);"
-  void setOsrParams(std::string);
-
-  void  runOsr(std::string);
-
-  //! Prints "olr(var_list_,olr_inst_,obj_var_,optim_weights_);"
-  void setOlrInst(std::string);
-
-  void  runOlr(std::string);
-
-  void  BeginOptimWeights(void);
-
-  void  setOptimWeights(std::string, std::string);
-
-  void  setOptimWeights(std::string, std::string, std::string);
-
-  void    runDynasave(std::string,std::string,std::string);
-
-  void    runDynatype(std::string,std::string,std::string);
-
-  void    beginModelComparison(void);
+class EstimatedParamsInitStatement : public Statement
+{
+private:
+  const vector<EstimationParams> estim_params_list;
+  const SymbolTable &symbol_table;
+public:
+  EstimatedParamsInitStatement(const vector<EstimationParams> &estim_params_list_arg,
+                               const SymbolTable &symbol_table_arg);
+  virtual void writeOutput(ostream &output) const;
+};
 
-  void    addMcFilename(std::string,std::string);
+class EstimatedParamsBoundsStatement : public Statement
+{
+private:
+  const vector<EstimationParams> estim_params_list;
+  const SymbolTable &symbol_table;
+public:
+  EstimatedParamsBoundsStatement(const vector<EstimationParams> &estim_params_list_arg,
+                                 const SymbolTable &symbol_table_arg);
+  virtual void writeOutput(ostream &output) const;
+};
 
-  void    runModelComparison(void);
+class OptimWeightsStatement : public Statement
+{
+public:
+  typedef map<string, string, less<string> > var_weights_type;
+  typedef map<pair<string, string>, string, less<pair<string, string> > > covar_weights_type;
+private:
+  const var_weights_type var_weights;
+  const covar_weights_type covar_weights;
+  const SymbolTable &symbol_table;
+public:
+  OptimWeightsStatement(const var_weights_type &var_weights_arg,
+                        const covar_weights_type &covar_weights_arg,
+                        const SymbolTable &symbol_table_arg);
+  virtual void writeOutput(ostream &output) const;
+};
 
-  void  set(void);
+class CalibStatement : public Statement
+{
+private:
+  const int covar;
+public:
+  CalibStatement(int covar_arg);
+  virtual void writeOutput(ostream &output) const;
+};
 
-  std::string get(void);
+class CalibVarStatement : public Statement
+{
+public:
+  //! Maps a variable to a pair (weight, expression)
+  typedef map<string, pair<string, string>, less<string> > calib_var_type;
+  //! Maps a pair of variables to a pair (weight, expression)
+  typedef map<pair<string, string>, pair<string, string>, less<pair<string, string> > > calib_covar_type;
+  //! Maps a pair (variable, autocorr) to a pair (weight, expression)
+  typedef map<pair<string, int>, pair<string, string>, less<pair<string, int> > > calib_ac_type;
+private:
+  const calib_var_type calib_var;
+  const calib_covar_type calib_covar;
+  const calib_ac_type calib_ac;
+  const SymbolTable &symbol_table;
+public:
+  CalibVarStatement(const calib_var_type &calib_var_arg,
+                    const calib_covar_type &calib_covar_arg,
+                    const calib_ac_type &calib_ac_arg,
+                    const SymbolTable &symbol_table_arg);
+  virtual void writeOutput(ostream &output) const;
 };
-//------------------------------------------------------------------------------
+
 #endif
diff --git a/parser.src/include/ModFile.hh b/parser.src/include/ModFile.hh
index 1c989857e8d3b728dbcb45330384e892dd4fa939..20c255a58f81fc43a4a2fde964b3d1020321045f 100644
--- a/parser.src/include/ModFile.hh
+++ b/parser.src/include/ModFile.hh
@@ -1,44 +1,50 @@
 #ifndef _MOD_FILE_HH
 #define _MOD_FILE_HH
 
+using namespace std;
+
+#include <ostream>
+
 #include "ModelParameters.hh"
 #include "SymbolTable.hh"
 #include "NumericalConstants.hh"
-#include "NumericalInitialization.hh"
-#include "Shocks.hh"
-#include "SigmaeInitialization.hh"
-#include "ComputingTasks.hh"
 #include "ModelTree.hh"
 #include "VariableTable.hh"
+#include "Statement.hh"
 
 //! The abstract representation of a "mod" file
 class ModFile
 {
 public:
-  //! Constructor
   ModFile();
+  ~ModFile();
   //! Model parameters
   ModelParameters model_parameters;
   //! Symbol table
   SymbolTable symbol_table;
   //! Variable table
   VariableTable variable_table;
-  //! Numerical constants
+  //! Numerical constants table
   NumericalConstants num_constants;
-  //! Numerical initalisations
-  NumericalInitialization numerical_initialization;
-  //! Shocks instructions
-  Shocks shocks;
-  //! Sigma_e instructions
-  SigmaeInitialization sigmae;
-  //! Computing tasks instructions
-  ComputingTasks computing_tasks;
   //! Model equations and their derivatives
   ModelTree model_tree;
   //! Option order
   int order;
   //! Option linear
   int linear;
+private:
+  //! List of statements
+  vector<Statement *> statements;
+
+public:
+  //! Add a statement
+  void addStatement(Statement *st);
+  //! Writes Matlab/Scilab output files
+  /*!
+    \param basename The base name used for writing output files. Should be the name of the mod file without its extension
+    \param clear_all Should a "clear all" instruction be written to output ?
+  */
+  void writeOutputFiles(const string &basename, bool clear_all);
 };
 
 #endif // ! MOD_FILE_HH
diff --git a/parser.src/include/ModelTree.hh b/parser.src/include/ModelTree.hh
index ba5a40e279efe0b847798236604f827cbead9b01..16645ab9d2e76bb3ce5e2031c02ecce349f857ce 100644
--- a/parser.src/include/ModelTree.hh
+++ b/parser.src/include/ModelTree.hh
@@ -13,6 +13,7 @@
 #include <stack>
 #include <sstream>
 #include <fstream>
+#include <ostream>
 //------------------------------------------------------------------------------
 #include "SymbolTable.hh"
 #include "OperatorTable.hh"
@@ -34,8 +35,6 @@ private :
   std::ostringstream      StaticOutput;
   /*! Output for dynamic model */
   std::ostringstream      DynamicOutput;
-  /*! Output for main file */
-  std::ostringstream output;
   /*! Output file stream for static model */
   std::ofstream       mStaticModelFile;
   /*! Output file stream for dynamic model */
@@ -54,7 +53,6 @@ private :
   //! Reference to numerical constants table
   const NumericalConstants &num_constants;
 
-private :
   /*! Computes argument derivative */
   inline NodeID     DeriveArgument(NodeID iArg, Type iType, int iVarID);
   /*! Gets output argument of terminal token */
@@ -62,17 +60,8 @@ private :
   /*! Gets expression of part of model tree */
   inline std::string    getExpression(NodeID StartID, EquationType  iEquationType, int iEquationID = -1);
   inline int optimize(NodeID id);
-public :
-  /*! When Jacobian (vs endogenous) is writen this flag is set to true */
-  bool    computeJacobian;
-  /*! When Jacobian (vs endogenous and exogenous) is writen this flag is set to true */
-  bool    computeJacobianExo;
   /*! When Hessian  is writen this flag is set to true */
   bool    computeHessian;
-  /*! Constructor */
-  ModelTree(SymbolTable &symbol_table_arg, VariableTable &variable_table_arg, ModelParameters &mod_param_arg, const NumericalConstants &num_constants);
-  /*! Destructor */
-  ~ModelTree();
   /*! Opens output M files (1st and 2nd derivatives) */
   void OpenMFiles(std::string iModelFileName1, std::string iModelFileName2 = "");
   /*! Opens output C files (1st and 2nd derivatives) */
@@ -96,9 +85,19 @@ public :
   */
   std::string     setDynamicModel(void);
   /*! Writes initialization of various Matlab variables */
-  void    ModelInitialization(void);
-  /*! Returns string output for main file */
-  std::string get();
+  void ModelInitialization(std::ostream &output);
+
+public:
+  //! Constructor
+  ModelTree(SymbolTable &symbol_table_arg, VariableTable &variable_table_arg, ModelParameters &mod_param_arg, const NumericalConstants &num_constants);
+  //! Destructor
+  ~ModelTree();
+  //! When Jacobian (vs endogenous) is written this flag is set to true
+  bool computeJacobian;
+  //! When Jacobian (vs endogenous and exogenous) is written this flag is set to true
+  bool computeJacobianExo;
+  //! Writes model initialization to output and uses basename for dumping model static/dynamic files
+  void writeOutput(std::ostream &output, const std::string &basename, int order, int linear);
 };
 //------------------------------------------------------------------------------
 #endif
diff --git a/parser.src/include/NumericalInitialization.hh b/parser.src/include/NumericalInitialization.hh
index 241fa1484ca6d67b74c1d475eac1b2949aa560bc..32ecd7b5bbb9480d5688b5210643f81084aca19b 100644
--- a/parser.src/include/NumericalInitialization.hh
+++ b/parser.src/include/NumericalInitialization.hh
@@ -1,94 +1,78 @@
 #ifndef _NUMERICALINITIALIZATION_HH
 #define _NUMERICALINITIALIZATION_HH
-//------------------------------------------------------------------------------
-/*! \file
-  \version 1.0
-  \date 04/16/2004
-  \par This file defines the NumericalInitialization class.
-*/
-//------------------------------------------------------------------------------
-#include <iostream>
+
+using namespace std;
+
 #include <string>
-#include <sstream>
-//------------------------------------------------------------------------------
+#include <map>
+
 #include "SymbolTable.hh"
-//------------------------------------------------------------------------------
-/*! \class  NumericalInitialization
-  \brief  Handles numerical initialization of Endogenousous and Exogenousous variables.
-*/
-class NumericalInitialization
+#include "Statement.hh"
+
+class InitParamStatement : public Statement
 {
-private :
-  /*! Output of this class */
-  std::ostringstream  *output;
-  //! Reference to symbol table
+private:
+  const string param_name;
+  const string param_value;
   const SymbolTable &symbol_table;
-  //! Reference to model parameters
-  const ModelParameters &mod_param;
-public :
-  /*! Constructor */
-  NumericalInitialization(const SymbolTable &symbol_table_arg, const ModelParameters &mod_param_arg);
-  /*! Destrcutor */
-  ~NumericalInitialization();
-  /*! Pointer to error function of parser class */
-  void (* error) (const char* m);
-  /*! Sets reference to external output string */
-  void setOutput(std::ostringstream* iOutput);
-  /*!
-    \param name a string.
-    \param expression an Expression type.
-    \par Description
-    Set constant value
-    - in interpretative languages (Matlab, Gauss, Scilab)\n
-    - in C++, evaluate expression and set value for Name in Parameters_Table\n
-  */
-  void SetConstant(std::string name,std::string expression);
-  /*!
-    \par Description
-    Initializes an initval block to set initial values for variables
-  */
-  void BeginInitval(void);
-  /*!
-    \param name a string.
-    \param expression an Expression type.
-    \par Description
-    Set initial or terminal value for a variable :
-    - name must be searched for in SymbolTable
-    - error if symbol isn't a variable (Endogenousous or Exogenousous or Exogenousous deterministic or recursive)
-    - write "oo_.steady_state(symbol.id) = expression;" (for Endogenousous variables) or "oo_.exo_steady_state(symbol.id) = expression;" (for Exogenousous variables)
-  */
-  void SetInit(std::string name,std::string expression);
-  /*!
-    \par Description
-    Closes an initval block.
-  */
-  void EndInitval(void);
-  /*!
-    \par Description
-    Initializes an endval block to set terminal values for variables.
-  */
-  void BeginEndval(void);
-  /*!
-    \par Description
-    Closes an endval block.
-  */
-  void EndEndval(void);
+public:
+  InitParamStatement(const string &param_name_arg, const string &param_value_arg,
+                     const SymbolTable &symbol_table_arg);
+  virtual void writeOutput(ostream &output) const;
+};
+
+class InitOrEndValStatement : public Statement
+{
+public:
   /*!
-    \par Description
-    Initializes an histval block to set historical values for variables.
+    We use a vector instead of a map, since the order of declaration matters:
+    an initialization can depend on a previously initialized variable inside the block
   */
-  void BeginHistval(void);
+  typedef vector<pair<string, string> > init_values_type;
+protected:
+  const init_values_type init_values;
+  const SymbolTable &symbol_table;
+public:
+  InitOrEndValStatement(const init_values_type &init_values_arg,
+                        const SymbolTable &symbol_table_arg);
+protected:
+  void writeInitValues(ostream &output) const;
+};
+
+class InitValStatement : public InitOrEndValStatement
+{
+private:
+  const ModelParameters &mod_param;
+public:
+  InitValStatement(const init_values_type &init_values_arg,
+                   const SymbolTable &symbol_table_arg,
+                   const ModelParameters &mod_param_arg);
+  virtual void writeOutput(ostream &output) const;
+};
+
+class EndValStatement : public InitOrEndValStatement
+{
+public:
+  EndValStatement(const init_values_type &init_values_arg,
+                  const SymbolTable &symbol_table_arg);
+  virtual void writeOutput(ostream &output) const;
+};
+
+class HistValStatement : public Statement
+{
+public:
   /*!
-    \param name a string.
-    \param lag a string.
-    \param expression an Expression type.
-    \par Description
-    Set historical for a variable :
-    - name must be searched for in SymbolTable
-    - error if symbol isn't a variable (Endogenousous or Exogenousous or Exogenousous deterministic or recursive)
-    - write "oo_.y_simul(symbol.id,MP.max_lag+lag+1) = expression;" (for Endogenousous variables) or "oo_.exo_simul(MP.max_lag+lag+1,symbol.id) = expression;" (for Exogenousous variables) (lag is <= 0, MP.max_lag will be set in ModelTree)
+    Contrary to Initval and Endval, we use a map, since it is impossible to reuse
+    a given initialization value in a second initialization inside the block.
   */
-  void SetHist(std::string name,int lag, std::string expression);
+  typedef map<pair<string, int>, string> hist_values_type;
+private:
+  const hist_values_type hist_values;
+  const SymbolTable &symbol_table;
+public:
+  HistValStatement(const hist_values_type &hist_values_arg,
+                   const SymbolTable &symbol_table_arg);
+  virtual void writeOutput(ostream &output) const;
 };
-//------------------------------------------------------------------------------
+
 #endif
diff --git a/parser.src/include/OutputFile.hh b/parser.src/include/OutputFile.hh
deleted file mode 100644
index a837e531b2276550801c4bcfa98f5e8a9c9d35d2..0000000000000000000000000000000000000000
--- a/parser.src/include/OutputFile.hh
+++ /dev/null
@@ -1,38 +0,0 @@
-#ifndef _OUTPUTFILE_HH
-#define _OUTPUTFILE_HH
-//------------------------------------------------------------------------------
-/*! \file
-  \version 1.0
-  \date 04/26/2004
-  \par This file defines the OutputFile class.
-*/
-//------------------------------------------------------------------------------
-#include <fstream>
-#include <string>
-
-#include "ModFile.hh"
-//------------------------------------------------------------------------------
-/*!
-  \class  OutputFile
-  \brief  Handles opening, writing and closing of output file
-*/
-class OutputFile
-{
-private :
-  /*! Output file stream */
-  ofstream  mOutputFile;
-public :
-  /*! Flag if set to true, execute "clear all"*/
-  bool    clear_all;
-
-  /*! Constructor */
-  OutputFile();
-  /*! Destructor */
-  ~OutputFile();
-  /*! Opens a given file and writes some initialization */
-  void Open(std::string iFileName, ModFile *mod_file);
-  /*! Writes output data from SymbolTable and passed strings to output file */
-  void Save(std::ostringstream& iOutput, ModFile *mod_file);
-};
-//------------------------------------------------------------------------------
-#endif
diff --git a/parser.src/include/ParsingDriver.hh b/parser.src/include/ParsingDriver.hh
index 8b3e93346a537f9200b17af5ee996fc23ffbf04d..968aaec16e1566a537c9eed7e57e954cd52dc058 100644
--- a/parser.src/include/ParsingDriver.hh
+++ b/parser.src/include/ParsingDriver.hh
@@ -8,6 +8,9 @@
 #include "TmpSymbolTable.hh"
 #include "DynareBison.hh"
 #include "ComputingTasks.hh"
+#include "Shocks.hh"
+#include "SigmaeInitialization.hh"
+#include "NumericalInitialization.hh"
 
 using namespace std;
 
@@ -39,19 +42,59 @@ private:
   //! Checks that a given symbol exists, and stops with an error message if it doesn't
   void check_symbol_existence(const string &name);
 
+  //! Creates option "optim_opt" in OptionsList if it doesn't exist, else add a comma, and adds the option name
+  void optim_options_helper(const string &name);
+
   //! Stores expressions
   Expression expression;
   //! Stores temporary symbol table
-  TmpSymbolTable tmp_symbol_table;
+  TmpSymbolTable *tmp_symbol_table;
+  //! Stores options lists
+  OptionsList options_list;
   //! Stores operator table
   OperatorTable op_table;
+  //! Temporary storage for trend elements
+  ObservationTrendsStatement::trend_elements_type trend_elements;
+  //! Temporary storage for filename list of ModelComparison
+  ModelComparisonStatement::filename_list_type filename_list;
+  //! Temporary storage for list of EstimationParams (from estimated_params* statements)
+  vector<EstimationParams> estim_params_list;
+  //! Temporary storage of variances from optim_weights
+  OptimWeightsStatement::var_weights_type var_weights;
+  //! Temporary storage of covariances from optim_weights
+  OptimWeightsStatement::covar_weights_type covar_weights;
+  //! Temporary storage of variances from calib_var
+  CalibVarStatement::calib_var_type calib_var;
+  //! Temporary storage of covariances from calib_var
+  CalibVarStatement::calib_covar_type calib_covar;
+  //! Temporary storage of autocorrelations from calib_var
+  CalibVarStatement::calib_ac_type calib_ac;
+  //! Temporary storage for deterministic shocks
+  ShocksStatement::det_shocks_type det_shocks;
+  //! Temporary storage for periods of deterministic shocks
+  vector<pair<int, int> > det_shocks_periods;
+  //! Temporary storage for values of deterministic shocks
+  vector<string> det_shocks_values;
+  //! Temporary storage for variances of shocks
+  ShocksStatement::var_and_std_shocks_type var_shocks;
+  //! Temporary storage for standard errors of shocks
+  ShocksStatement::var_and_std_shocks_type std_shocks;
+  //! Temporary storage for covariances of shocks
+  ShocksStatement::covar_and_corr_shocks_type covar_shocks;
+  //! Temporary storage for correlations of shocks
+  ShocksStatement::covar_and_corr_shocks_type corr_shocks;
+  //! Temporary storage for Sigma_e rows
+  SigmaeStatement::row_type sigmae_row;
+  //! Temporary storage for Sigma_e matrix
+  SigmaeStatement::matrix_type sigmae_matrix;
+  //! Temporary storage for initval/endval blocks
+  InitOrEndValStatement::init_values_type init_values;
+  //! Temporary storage for histval blocks
+  HistValStatement::hist_values_type hist_values;
 
   //! The mod file representation constructed by this ParsingDriver
   ModFile *mod_file;
 
-  //! Reference to output string
-  ostringstream *output;
-
 public:
   //! Constructor
   ParsingDriver();
@@ -94,10 +137,6 @@ public:
     exit(-1);
   }
 
-  //! Sets reference to output string
-  void setoutput(ostringstream *ostr);
-  //! Executes final instructions
-  void finish();
   //! Check if a given symbol exists in the parsing context
   bool exists_symbol(const char *s);
   //! Sets variable offset of ModelTree class to use C output
@@ -128,26 +167,28 @@ public:
   ExpObj *add_expression_token(ExpObj *arg1, int op);
   //! Adds a unary token to an expression, with function name unknown
   ExpObj *add_expression_token(ExpObj *arg1, string *op_name);
+  //! Adds a "periods" statement
+  void periods(string *periods);
+  //! Adds a "dsample" statement
+  void dsample(string *arg1);
+  //! Adds a "dsample" statement
+  void dsample(string *arg1, string *arg2);
   //! Writes parameter intitialisation expression
   void init_param(string *name, ExpObj *rhs);
   //! Writes an initval block
   void init_val(string *name, ExpObj *rhs);
   //! Writes an histval block
   void hist_val(string *name, string *lag, ExpObj *rhs);
-  //! Writes begining of an initval block
-  void begin_initval();
   //! Writes end of an initval block
   void end_initval();
-  //! Writes begining of an endval block
-  void begin_endval();
   //! Writes end of an endval block
   void end_endval();
-  //! Writes begining of an histval block
-  void begin_histval();
-  //! Write begining of a shock block
-  void begin_shocks();
-  void begin_mshocks();
+  //! Writes end of an histval block
+  void end_histval();
+  //! Writes a shocks statement
   void end_shocks();
+  //! Writes a mshocks statement
+  void end_mshocks();
   //! Adds a deterministic chock
   void add_det_shock(string *var);
   //! Adds a std error chock
@@ -186,6 +227,8 @@ public:
   void option_str(const string &name_option, string *opt);
   //! Sets an option to a string value
   void option_str(const string &name_option, const string &opt);
+  //! Indicates that the model is linear
+  void linear();
   //! Adds a variable to temp symbol table and sets its value
   void add_tmp_var(string *tmp_var1, string *tmp_var2);
   //! Adds a variable to temp symbol table
@@ -198,40 +241,43 @@ public:
   void simul();
   //! Writes check command
   void check();
-  //! Writes instructions for estimation initialization
-  void estimation_init();
-  //! Writes instructions for estimated elements
-  void set_estimated_elements();
-  void set_estimated_init_elements();
-  void set_estimated_bounds_elements();
+  //! Writes estimated params command
+  void estimated_params();
+  //! Writes estimated params init command
+  void estimated_params_init();
+  //! Writes estimated params bound command
+  void estimated_params_bounds();
+  //! Add a line in an estimated params block
+  void add_estimated_params_element();
   //! Runs estimation process
   void run_estimation();
   //! Runs prior_analysis();
   void run_prior_analysis();
   //! Runs posterior_analysis();
   void run_posterior_analysis();
-  //! Prints optimization options
-  void optim_options(string *str1, string *str2, int task);
+  //! Adds an optimization option (string value)
+  void optim_options_string(string *name, string *value);
+  //! Adds an optimization option (numeric value)
+  void optim_options_num(string *name, string *value);
   //! Prints varops instructions
   void set_varobs();
-  void set_trend_init();
+  void set_trends();
   void set_trend_element(string *arg1, ExpObj *arg2);
   void set_unit_root_vars();
-  void begin_optim_weights();
-  void set_optim_weights(string *arg1, ExpObj *arg2);
-  void set_optim_weights(string *arg1, string *arg2, ExpObj *arg3);
+  void optim_weights();
+  void set_optim_weights(string *name, ExpObj *value);
+  void set_optim_weights(string *name1, string *name2, ExpObj *value);
   void set_osr_params();
   void run_osr();
   void set_olr_inst();
   void run_olr();
-  void begin_calib_var();
+  void run_calib_var();
   void set_calib_var(string *name, string *weight, ExpObj *expression);
-  void set_calib_var(string *name1, string *name2, string *weight, ExpObj *expression);
+  void set_calib_covar(string *name1, string *name2, string *weight, ExpObj *expression);
   void set_calib_ac(string *name, string *ar, string *weight, ExpObj *expression);
-  void run_calib(int);
+  void run_calib(int covar);
   void run_dynasave(string *arg1, string *arg2 = new string);
   void run_dynatype(string *arg1, string *arg2 = new string);
-  void begin_model_comparison();
   void add_mc_filename(string *filename, string *prior = new string("1"));
   void run_model_comparison();
   //! Writes token "arg1=arg2" to model tree
diff --git a/parser.src/include/Shocks.hh b/parser.src/include/Shocks.hh
index 1f760004e58cfd705e01cf7f5b51ebfde0d365b6..57494d7da75b82f31ba21912f30e4cb57f26ab3f 100644
--- a/parser.src/include/Shocks.hh
+++ b/parser.src/include/Shocks.hh
@@ -1,85 +1,69 @@
 #ifndef _SHOCKS_HH
 #define _SHOCKS_HH
-//------------------------------------------------------------------------------
-/** \file
- * \version 1.0
- * \date 04/26/2004
- * \par This file defines the Shocks class.
- */
-//------------------------------------------------------------------------------
+
+using namespace std;
+
 #include <string>
-#include <sstream>
-#include <list>
 #include <vector>
+#include <map>
 
-//------------------------------------------------------------------------------
-/*!
-  \class  Shocks
-  \brief  Handles Shocks command
-*/
-class Shocks
+#include "Statement.hh"
+#include "SymbolTable.hh"
+
+class AbstractShocksStatement : public Statement
 {
-private :
-  int mshock_flag;
-  int exo_det_length;
-  /*!
-    \class ShockElement
-    \brief Shock element strcuture
-  */
-  struct ShockElement
+public:
+  struct DetShockElement
   {
-    std::string period1;
-    std::string period2;
-    std::string value;
+    int period1;
+    int period2;
+    string value;
   };
-  /*!
-    \class Shock
-    \brief Shock Structure
-  */
-  struct Shock
-  {
-    int     id1;
-    int     id2;
-    std::list<ShockElement> shock_elems;
-    std::string value;
-  };
-  /*! Output string of this class */
-  std::ostringstream    *output;
-  /*! Vector of begin period range */
-  std::vector<std::string>    mPeriod1;
-  /*! vector of end period range */
-  std::vector<std::string>    mPeriod2;
-  /*! vector of shock values */
-  std::vector<std::string>    mValues;
-public :
-  /*! Constructor */
-  Shocks();
-  /*! Destructor */
-  ~Shocks();
-  /*! Pointer to error function of parser class */
-  void (* error) (const char* m);
-  /*!
-    Set output reference
-    \param iOutput : reference to an ostringstream
-  */
-  void  setOutput(std::ostringstream* iOutput);
-  /*! Initialize shocks (or mshocks, for multiplicative shocks) block */
-  void  BeginShocks(void);
-  void  BeginMShocks(void);
-  void  EndShocks(void);
-  /*! Sets a shock */
-  void  AddDetShockExo(int id1);
-  void  AddDetShockExoDet(int id1);
-  void  AddSTDShock(int id1, std::string value);
-  void  AddVARShock(int id1, std::string value);
-  void  AddCOVAShock(int id1, int id2 , std::string value);
-  void  AddCORRShock(int id1, int id2 , std::string value);
-  /*! Adds a period rage */
-  void  AddPeriod(std::string p1, std::string p2);
-  /*! Adds a period */
-  void  AddPeriod(std::string p1);
-  /*! Adds a value */
-  void  AddValue(std::string value);
+  typedef map<string, vector<DetShockElement> > det_shocks_type;
+  typedef map<string, string> var_and_std_shocks_type;
+  typedef map<pair<string, string>, string> covar_and_corr_shocks_type;
+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:
+  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);
+  virtual void writeOutput(ostream &output) const;
+};
+
+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;
+};
+
 #endif
diff --git a/parser.src/include/SigmaeInitialization.hh b/parser.src/include/SigmaeInitialization.hh
index 79a6fa47c8a40e3725f236cd670357d6e615d5cb..45ee65cc3f640de794e6f0614bb0a1d538e10e52 100644
--- a/parser.src/include/SigmaeInitialization.hh
+++ b/parser.src/include/SigmaeInitialization.hh
@@ -1,65 +1,45 @@
 #ifndef _SIGMAEINITIALIZATION_HH
 #define _SIGMAEINITIALIZATION_HH
-//------------------------------------------------------------------------------
-/*! \file
-  \version 1.0
-  \date 04/13/2004
-  \par This file defines the SigmaeInitialization class.
-*/
-//------------------------------------------------------------------------------
+
+using namespace std;
+
 #include <string>
-#include <sstream>
 #include <vector>
-//------------------------------------------------------------------------------
-/*!
-  \class  SigmaeInitialization
-  \brief  Handles Sigma_e command
-*/
-class SigmaeInitialization
-{
 
-  /*! Matrix type enum */
-  enum MatrixType
+#include "Statement.hh"
+
+//! Stores a Sigma_e statement
+class SigmaeStatement : public Statement
+{
+public:
+  //! Matrix form (lower or upper triangular) enum
+  enum matrix_form_type
     {
-      eLower = 0,                  //!< Lower matrix
-      eUpper = 1                   //!< Upper matrix
+      eLower = 0,              //!< Lower triangular matrix
+      eUpper = 1               //!< Upper triangular matrix
     };
-private :
-  /*!  A row of Sigma_e */
-  std::vector<std::string>        row;
-  //!  The hole matrix Sigma_e */
-  std::vector<std::vector<std::string> >    matrix;
-  /*! Output of this class */
-  std::ostringstream        *output;
-  /*! Matrix type (eLower(lower) or eUpper) */
-  MatrixType          type;
+  //! Type of a matrix row
+  typedef vector<string> row_type;
+  //! Type of a complete matrix
+  typedef vector<row_type> matrix_type;
+
+  //! An exception indicating that a matrix is neither upper triangular nor lower triangular
+  class MatrixFormException
+  {
+  };
+private:
+  //! The matrix
+  const matrix_type matrix;
+  //! Matrix form (lower or upper)
+  const matrix_form_type matrix_form;
 
-  /*! Check that the matrix is triangular or square */
-  void  CheckMatrix(void);
-  /*! Print matrix to output */
-  void  SetMatrix(void);
+  //! Returns the type (upper or lower triangular) of a given matrix
+  /*! Throws an exception if it is neither upper triangular nor lower triangular */
+  static matrix_form_type determineMatrixForm(const matrix_type &matrix) throw (MatrixFormException);
 
 public :
-  /*! Constructor */
-  SigmaeInitialization();
-  /*! Destructor */
-  ~SigmaeInitialization();
-  /*! Pointer to error function of parser class */
-  void (* error) (const char* m);
-  /*!
-    Set output reference
-    \param iOutput : reference to an ostringstream
-  */
-  void  setOutput(std::ostringstream* iOutput);
-  /*!
-    Add an expression to current row
-    \param expression : a string expression
-  */
-  void  AddExpression(std::string expression);
-  /*! Add current row to matrix and clears the row */
-  void  EndOfRow();
-  /*! Check matrix and print it to output */
-  void  set(void);
+  SigmaeStatement(const matrix_type &matrix_arg) throw (MatrixFormException);
+  virtual void writeOutput(ostream &output) const;
 };
-//------------------------------------------------------------------------------
+
 #endif
diff --git a/parser.src/include/Statement.hh b/parser.src/include/Statement.hh
new file mode 100644
index 0000000000000000000000000000000000000000..fb170ba49849a3457c47108d410e58ac31c0f82c
--- /dev/null
+++ b/parser.src/include/Statement.hh
@@ -0,0 +1,39 @@
+#ifndef _STATEMENT_HH
+#define _STATEMENT_HH
+
+using namespace std;
+
+#include <ostream>
+#include <string>
+#include <map>
+
+class Statement
+{
+public:
+  virtual ~Statement();
+  virtual void writeOutput(ostream &output) const = 0;
+};
+
+class NativeStatement : public Statement
+{
+private:
+  const string native_statement;
+public:
+  NativeStatement(const string &native_statement_arg);
+  virtual void writeOutput(ostream &output) const;
+};
+
+class OptionsList
+{
+public:
+  typedef map<string, string> num_options_type;
+  typedef map<string, pair<string, string> > paired_num_options_type;
+  typedef map<string, string> string_options_type;
+  num_options_type num_options;
+  paired_num_options_type paired_num_options;
+  string_options_type string_options;
+  void writeOutput(ostream &output) const;
+  void clear();
+};
+
+#endif // ! _STATEMENT_HH
diff --git a/parser.src/include/SymbolTable.hh b/parser.src/include/SymbolTable.hh
index 2a158c3aea1c40c992f97e0756dcfcd25d88f6b3..1338d356a11799f06755d3fc6a2b66bdf4252e28 100644
--- a/parser.src/include/SymbolTable.hh
+++ b/parser.src/include/SymbolTable.hh
@@ -10,6 +10,7 @@
 #include <map>
 #include <string>
 #include <vector>
+#include <ostream>
 //------------------------------------------------------------------------------
 #include "ModelParameters.hh"
 #include "SymbolTableTypes.hh"
@@ -75,8 +76,8 @@ public :
   inline Type getType(const std::string &name) const;
   /*! Gets ID by name */
   inline int getID(const std::string &name) const;
-  /*! Gets output string of this class */
-  std::string get();
+  //! Write output of this class
+  void writeOutput(std::ostream &output);
 #if 0 // Commented out on 27/11/2006, SV
   /*! Checks if symbols are used in model equations, removes unused symbol */
   void clean();
diff --git a/parser.src/include/TmpSymbolTable.hh b/parser.src/include/TmpSymbolTable.hh
index 8b2fbaec1eef2b8603644624beedf7839df8789d..b0d2dc864a8e5d16fc195999852f78b7b7253601 100644
--- a/parser.src/include/TmpSymbolTable.hh
+++ b/parser.src/include/TmpSymbolTable.hh
@@ -7,8 +7,9 @@
  * \par This file defines the TmpSymbolTable class.
  */
 //------------------------------------------------------------------------------
-#include <list>
-#include <sstream>
+#include <string>
+#include <vector>
+#include <ostream>
 
 #include "SymbolTable.hh"
 
@@ -20,32 +21,24 @@ class TmpSymbolTable
 {
 private :
   /*! list of string TempSymbolTable */
-  std::list<std::string>  tmpsymboltable;
+  std::vector<std::string> tmpsymboltable;
   /*! List of symbol Values */
-  std::list<std::string>  NameTable;
-  /*! Output of this class */
-  std::ostringstream  output;
-  //! Pointer to global symbol table
-  SymbolTable *symbol_table;
+  std::vector<std::string> nameTable;
+  //! A reference to enclosing symbol table
+  const SymbolTable &symbol_table;
 public :
   /*! Constrcutor */
-  TmpSymbolTable();
+  TmpSymbolTable(const SymbolTable &symbol_table_arg);
   /*! Destructor*/
   ~TmpSymbolTable();
-  //! Sets global symbol table pointer
-  void setGlobalSymbolTable(SymbolTable *symbol_table_arg);
-  /*! Pointer to error function of parser class */
-  void (* error) (const char* m);
   /*! Adds a temp symbol */
-  void  AddTempSymbol(std::string symbol);
+  void AddTempSymbol(const std::string &symbol);
   /*! Adds a temp symbol and its value */
-  void  AddTempSymbol(std::string symbol1, std::string symbol2);
+  void AddTempSymbol(const std::string &symbol1, const std::string &symbol2);
   /*! Write TempSymbolTable to output string */
-  void  set(std::string varname);
-  /*! Gets size of TempSymbolTable */
-  int   size(void);
-  /*! Gets output of this class*/
-  std::string  get(void);
+  void writeOutput(const std::string &varname, std::ostream &output) const;
+  //! Clears all content
+  void clear();
 };
 //------------------------------------------------------------------------------
 #endif