diff --git a/doc/manual.xml b/doc/manual.xml
index 7c1a32b307c8aee0f6a42588ea9cbefd04614e75..bdce693bf4080a60c08749f3a3bc90bc48c3c74e 100644
--- a/doc/manual.xml
+++ b/doc/manual.xml
@@ -98,7 +98,7 @@ Currently the development team of Dynare is composed of S. Adjemian, H. Bastani,
 <chapter><title>Installation and configuration</title>
 
 
-<sect1><title>Software requirements</title>
+<sect1 id="software-requirements"><title>Software requirements</title>
 <para>
 Packaged versions of Dynare are available for <trademark class="registered">Windows</trademark> XP/Vista, <ulink url="http://www.debian.org">Debian GNU/Linux</ulink> and <ulink url="http://www.ubuntu.com/">Ubuntu</ulink>.
 Dynare should work on other systems, but some compilation steps are necessary in that case.
@@ -109,13 +109,10 @@ Dynare should work on other systems, but some compilation steps are necessary in
 <listitem><para>GNU Octave version 3.0.0 or above.</para></listitem>
 </itemizedlist>
 </para>
-</sect1>
-
-<sect1><title>Installation of GNU Octave</title>
 
-<para>You can skip this step if you are planning to use only <trademark class="registered">MATLAB</trademark> with Dynare.</para>
+<para>Some installation instructions for GNU Octave can be found on <ulink url="http://www.dynare.org/DynareWiki/DynareOctave">Dynare Wiki</ulink>.</para>
 
-<para>Please refer to <ulink url="http://www.dynare.org/DynareWiki/DynareOctave">Dynare Wiki</ulink> for detailed instructions.</para>
+<para>If you are using MATLAB for Windows, and if you plan to use options <xref linkend="use_dll"/> or <xref linkend="k_order_solver"/>, you will need to install a C++ compiler on your machine. The easiest solution is to install <ulink url="http://www.microsoft.com/Express/VC/"><trademark class="registered">Microsoft</trademark> Visual C++ 2008 Express Edition</ulink>, and then to type <literal>mex -setup</literal> on the MATLAB prompt (it should autodetect the compiler). For users of MATLAB 64-bit, please refer to <ulink url="http://www.mathworks.fr/support/solutions/en/data/1-6IJJ3L/index.html?solution=1-6IJJ3L">these instructions</ulink>. Users of MATLAB under Linux and MacOS, and users of GNU Octave normally need to do nothing, since a working compilation environment is available by default.</para>
 
 </sect1>
 
@@ -840,9 +837,9 @@ Inside the model block, Dynare allows the creation of <emphasis>model-local vari
       <term><option>linear</option></term>
       <listitem><para>Declares the model as being linear. It spares oneself from having to declare initial values for computing the steady state, and it sets automatically <option>order</option><literal>=1</literal> in <xref linkend="stoch_simul" />.</para></listitem>
     </varlistentry>
-    <varlistentry>
+    <varlistentry id="use_dll">
       <term><option>use_dll</option></term>
-      <listitem><para>Instructs the preprocessor to create dynamic loadable libraries (DLL) containing the model equations and derivatives, instead of writing those in <filename class="extension">M</filename>-files. You need to having a working compilation environment (<foreignphrase>i.e.</foreignphrase> the <command>mex</command> command of <trademark class="registered">MATLAB</trademark> or Octave must be operational). Using this option can result in faster simulations or estimations, at the expense of some initial compilation time.<footnote><para>In particular, for big models, the compilation step can be very time-consuming, and use of this option may be counter-productive in those cases.</para></footnote></para></listitem>
+      <listitem><para>Instructs the preprocessor to create dynamic loadable libraries (DLL) containing the model equations and derivatives, instead of writing those in <filename class="extension">M</filename>-files. You need a working compilation environment, <foreignphrase>i.e.</foreignphrase> a working <literal>mex</literal> command (see <xref linkend="software-requirements"/> for more details). Using this option can result in faster simulations or estimations, at the expense of some initial compilation time.<footnote><para>In particular, for big models, the compilation step can be very time-consuming, and use of this option may be counter-productive in those cases.</para></footnote></para></listitem>
     </varlistentry>
   </variablelist>
 </refsect1>
@@ -1032,7 +1029,7 @@ For models with lags on more than one period, the command <xref linkend='histval
 
 <para>This steady state will be used as the initial condition at all the periods preceeding the first simulation period for the two possible types of simulations in stochastic mode:</para>
 <itemizedlist>
-  <listitem><para>in <xref linkend="stoch_simul"/>, if the <option>simul</option> or <option>periods</option> options are specified</para></listitem>
+  <listitem><para>in <xref linkend="stoch_simul"/>, if the <option>periods</option> options is specified</para></listitem>
   <listitem><para>in <xref linkend="forecast"/> (in this case, note that it is still possible to modify some of these initial values with <xref linkend="histval"/>)</para></listitem>
 </itemizedlist>
 
@@ -1789,7 +1786,7 @@ The simulated endogenous variables are available in global matrix <varname>oo_.e
   </varlistentry>
   <varlistentry>
     <term><option>hp_filter</option> = <replaceable>INTEGER</replaceable></term>
-    <listitem><para>Uses HP filter with &lambda; = <replaceable>INTEGER</replaceable> before computing moments. Default: no filter</para></listitem>
+    <listitem><para>Uses HP filter with &lambda; = <replaceable>INTEGER</replaceable> before computing moments. Note that this option is currently not available when computing empirical moments (see <option>periods</option> option). Default: no filter</para></listitem>
   </varlistentry>
   <varlistentry>
     <term><option>hp_ngrid</option> = <replaceable>INTEGER</replaceable></term>
@@ -1833,11 +1830,15 @@ The simulated endogenous variables are available in global matrix <varname>oo_.e
   </varlistentry>
   <varlistentry id="order" xreflabel="order">
     <term><option>order = <replaceable>INTEGER</replaceable></option></term>
-    <listitem><para>Order of Taylor approximation. Acceptable values are <literal>1</literal> and <literal>2</literal>. Default: <literal>2</literal></para></listitem>
+    <listitem><para>Order of Taylor approximation. Acceptable values are <literal>1</literal>, <literal>2</literal> and <literal>3</literal>. Note that for third order, <option>k_order_solver</option> option is implied, and only empirical moments are available (you must provide a value for <option>periods</option> option). Default: <literal>2</literal></para></listitem>
   </varlistentry>
+    <varlistentry id="k_order_solver">
+      <term><option>k_order_solver</option></term>
+      <listitem><para>Use a k-order solver, implemented in C++, instead of the default Dynare solver. You need a working compilation environment, <foreignphrase>i.e.</foreignphrase> a working <literal>mex</literal> command (see <xref linkend="software-requirements"/> for more details). Default: disabled for order 1 and 2, enabled otherwise</para></listitem>
+    </varlistentry>
   <varlistentry>
     <term><option>periods</option> = <replaceable>INTEGER</replaceable></term>
-    <listitem><para>Specifies the number of periods to use in simulations. If <option>order</option>=<literal>1</literal>, no simulation is necessary to compute theoretical moments and IRFs. A number of periods larger than one triggers automatically option <option>simul</option>. Default: <literal>0</literal></para></listitem>
+    <listitem><para>If different from zero, empirical moments will be computed instead of theoretical moments. The value of the option specifies the number of periods to use in the simulations. Values of the <xref linkend='initval'/> block, possibly recomputed by <xref linkend='steady'/>, will be used as starting point for the simulation. The simulated endogenous variables are made available to the user in a vector for each variable and in the global matrix <varname>oo_.endo_simul</varname>. The variables in the <varname>oo_.endo_simul</varname> matrix, in their order of declaration (as in <varname>M_.endo_names</varname>) Default: <literal>0</literal></para></listitem>
   </varlistentry>
   <varlistentry>
     <term><option>qz_criterium</option> = <replaceable>DOUBLE</replaceable></term>
@@ -1847,10 +1848,6 @@ The simulated endogenous variables are available in global matrix <varname>oo_.e
     <term><option>replic</option> = <replaceable>INTEGER</replaceable></term>
     <listitem><para>Number of simulated series used to compute the IRFs. Default: <literal>1</literal> if <option>order</option>=<literal>1</literal>, and <literal>50</literal> otherwise</para></listitem>
   </varlistentry>
-  <varlistentry>
-    <term><option>simul</option></term>
-    <listitem><para>Computes a stochastic simulation of the model for the number of periods specified in the <option>periods</option> option. Uses <xref linkend='initval'/> values, possibly recomputed by <xref linkend='steady'/>, as initial values for the simulation. The simulated endogenous variables are made available to the user in a vector for each variable and in the global matrix <varname>oo_.endo_simul</varname>. The variables in the <varname>oo_.endo_simul</varname> matrix, in their order of declaration (as in <varname>M_.endo_names</varname>). Default: no simulation</para></listitem>
-  </varlistentry>
   <varlistentry>
     <term><option>simul_seed</option> = <replaceable>INTEGER</replaceable></term>
     <listitem><para>Specifies a seed for the random generator so as to obtain the same random sample at each run of the program. Otherwise a different sample is used for each run. Default: seed not specified</para></listitem>
diff --git a/matlab/disp_moments.m b/matlab/disp_moments.m
index a4dbf1f5d46891c284a3058ce11cf8dd94b9df04..734130861723bfb065b59dd308bfd8517159baaf 100644
--- a/matlab/disp_moments.m
+++ b/matlab/disp_moments.m
@@ -23,6 +23,10 @@ function disp_moments(y,var_list)
   warning_old_state = warning;
   warning off
 
+  if options_.hp_filter
+      error('STOCH_SIMUL: HP filter is not yet implemented for empirical moments')
+  end
+  
   if size(var_list,1) == 0
       var_list = M_.endo_names(1:M_.orig_endo_nbr, :);
   end
diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index 467c6927eb82256dd577b1e32e26c692095cd555..3656583749a4798653f68fbc381bd7dead0d5e53 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -90,7 +90,6 @@ function global_initialization()
   options_.nocorr = 0;
   options_.periods = 0;
   options_.noprint = 0;
-  options_.simul = 0;
   options_.SpectralDensity = 0;
 
   % TeX output
diff --git a/matlab/osr.m b/matlab/osr.m
index 2e9a31e0341518642eed1c0df3777eaee3b7ad39..275763fc50ce93fe64072864d1231527577fb37d 100644
--- a/matlab/osr.m
+++ b/matlab/osr.m
@@ -32,12 +32,6 @@ function osr(var_list,params,i_var,W)
   options_ = set_default_option(options_,'hp_ngrid',512);
   options_ = set_default_option(options_,'simul',0);
   options_ = set_default_option(options_,'periods',1);
-  if options_.simul & ~isempty(options_.periods) & options_.periods == 0
-    options_.periods = options_.periods;
-  end
-
-  options_.periods = max(options_.periods,1);
-  options_.periods = options_.periods;
   
   make_ex_;
 
diff --git a/matlab/simul.m b/matlab/simul.m
index fd7b5dc360981004135c53a360ceb4e0d46d758e..f2c3df2c8912e60c391523ab7362ccd53b03fc3d 100644
--- a/matlab/simul.m
+++ b/matlab/simul.m
@@ -37,14 +37,10 @@ if size(M_.lead_lag_incidence,2)-nnz(M_.lead_lag_incidence(M_.maximum_endo_lag+1
   error (mess) ;
 end
 
-if ~isfield(options_,'periods') & ~isempty(options_.periods)
-  options_.periods = options_.periods
-end
 options_ = set_default_option(options_,'periods',0);
 if options_.periods == 0
   error('SIMUL: number of periods for the simulation isn''t specified')
 end
-options_.periods = options_.periods;
 
   if ~ options_.initval_file
       if ~isfield(options_,'datafile')
diff --git a/matlab/simult.m b/matlab/simult.m
index 728af37012e0e5109e15b98c8e5ef3dd13b5bb62..fba16268df2784543996e35dfa3335aafd434e68 100644
--- a/matlab/simult.m
+++ b/matlab/simult.m
@@ -39,7 +39,6 @@ if replic == 0
   replic = 1;
 end
 seed = options_.simul_seed;
-options_.periods = options_.periods;
 
 it_ = M_.maximum_lag + 1 ;
 
diff --git a/matlab/stoch_simul.m b/matlab/stoch_simul.m
index 6e068a2b1c9bdf2a5e18f7f3442b9abbc125a05a..f7c1321cb3c4e4c5c8bb8587ca24d804ae9efb72 100644
--- a/matlab/stoch_simul.m
+++ b/matlab/stoch_simul.m
@@ -26,7 +26,6 @@ function info=stoch_simul(var_list)
   if options_.order == 1
       options_.replic = 1;
   elseif options_.order == 3
-      options_.simul = 1;
       options_.k_order_solver = 1;
   end
   
@@ -70,12 +69,9 @@ function info=stoch_simul(var_list)
     end
   end
 
-  if options_.simul == 0 & options_.nomoments == 0
+  if options_.periods == 0 && options_.nomoments == 0
     disp_th_moments(oo_.dr,var_list); 
-  elseif options_.simul == 1
-    if options_.periods == 0
-      error('STOCH_SIMUL error: number of periods for the simulation isn''t specified')
-    end
+  elseif options_.periods ~= 0
     if options_.periods < options_.drop
       disp(['STOCH_SIMUL error: The horizon of simulation is shorter' ...
 	    ' than the number of observations to be DROPed'])
diff --git a/matlab/stoch_simul_sparse.m b/matlab/stoch_simul_sparse.m
index 04741fbd926876a62c41d850a1a4602427f5dd40..f5f47bec3cf63426cc0523e3300ba1ff2c8ed371 100644
--- a/matlab/stoch_simul_sparse.m
+++ b/matlab/stoch_simul_sparse.m
@@ -73,12 +73,9 @@ function info=stoch_simul_sparse(var_list)
     disp_dr_sparse(oo_.dr,options_.order,var_list);
   end
 
-  if options_.simul == 0 & options_.nomoments == 0
+  if options_.periods == 0 && options_.nomoments == 0
     disp_th_moments(oo_.dr,var_list); 
-  elseif options_.simul == 1
-    if options_.periods == 0
-      error('STOCH_SIMUL error: number of periods for the simulation isn''t specified')
-    end
+  elseif options_.periods ~= 0
     if options_.periods < options_.drop
       disp(['STOCH_SIMUL error: The horizon of simulation is shorter' ...
 	    ' than the number of observations to be DROPed'])
diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc
index 90f7e9a8749c65ee0bc3a153fa4fe5e91c07fa3b..f802f917d00816d443beb2b1d7b5c588ed9e510b 100644
--- a/preprocessor/ComputingTasks.cc
+++ b/preprocessor/ComputingTasks.cc
@@ -118,14 +118,6 @@ StochSimulStatement::checkPass(ModFileStructure &mod_file_struct)
   it = options_list.num_options.find("partial_information");
   if (it != options_list.num_options.end() && it->second == "1")
     mod_file_struct.partial_information = true;
-
-  // This (temporary) check is present in stoch_simul, osr and ramsey_policy
-  if (options_list.num_options.find("simul") != options_list.num_options.end()
-      && options_list.num_options.find("hp_filter") != options_list.num_options.end())
-    {
-      cerr << "ERROR: stoch_simul: HP filter is not yet implemented when computing empirical simulations" << endl;
-      exit(EXIT_FAILURE);
-    }
 }
 
 void
@@ -182,14 +174,6 @@ RamseyPolicyStatement::checkPass(ModFileStructure &mod_file_struct)
   it = options_list.num_options.find("partial_information");
   if (it != options_list.num_options.end() && it->second == "1")
     mod_file_struct.partial_information = true;
-
-  // This (temporary) check is present in stoch_simul, osr and ramsey_policy
-  if (options_list.num_options.find("simul") != options_list.num_options.end()
-      && options_list.num_options.find("hp_filter") != options_list.num_options.end())
-    {
-      cerr << "ERROR: ramsey_policy: HP filter is not yet implemented when computing empirical simulations" << endl;
-      exit(EXIT_FAILURE);
-    }
 }
 
 void
@@ -316,7 +300,6 @@ void
 PeriodsStatement::writeOutput(ostream &output, const string &basename) const
 {
   output << "options_.periods = " << periods << ";" << endl;
-  output << "options_.simul = 1;" << endl;
 }
 
 DsampleStatement::DsampleStatement(int val1_arg) : val1(val1_arg), val2(-1)
@@ -747,14 +730,6 @@ OsrStatement::checkPass(ModFileStructure &mod_file_struct)
   it = options_list.num_options.find("partial_information");
   if (it != options_list.num_options.end() && it->second == "1")
     mod_file_struct.partial_information = true;
-
-  // This (temporary) check is present in stoch_simul, osr and ramsey_policy
-  if (options_list.num_options.find("simul") != options_list.num_options.end()
-      && options_list.num_options.find("hp_filter") != options_list.num_options.end())
-    {
-      cerr << "ERROR: osr: HP filter is not yet implemented when computing empirical simulations" << endl;
-      exit(EXIT_FAILURE);
-    }
 }
 
 void
diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy
index 925dc3e0171e67db72b687b5846c54764fbfaa94..61ca875f42f8b64596dcbaf59413c352af85f940 100644
--- a/preprocessor/DynareBison.yy
+++ b/preprocessor/DynareBison.yy
@@ -1655,13 +1655,12 @@ o_nomoments : NOMOMENTS { driver.option_num("nomoments", "1"); };
 o_irf : IRF EQUAL INT_NUMBER { driver.option_num("irf", $3); };
 o_hp_filter : HP_FILTER EQUAL INT_NUMBER { driver.option_num("hp_filter", $3); };
 o_hp_ngrid : HP_NGRID EQUAL INT_NUMBER { driver.option_num("hp_ngrid", $3); };
-o_periods : PERIODS EQUAL INT_NUMBER
-            { driver.option_num("periods", $3); driver.option_num("simul", "1"); };
+o_periods : PERIODS EQUAL INT_NUMBER { driver.option_num("periods", $3); };
 o_cutoff : CUTOFF EQUAL number { driver.cutoff($3); }
 o_markowitz : MARKOWITZ EQUAL number { driver.option_num("markowitz", $3); };
 o_minimal_solving_periods : MINIMAL_SOLVING_PERIODS EQUAL number { driver.option_num("minimal_solving_periods", $3); };
 o_mfs : MFS EQUAL INT_NUMBER { driver.mfs($3); };
-o_simul : SIMUL { driver.option_num("simul", "1"); };
+o_simul : SIMUL; // Do nothing, only here for backward compatibility
 o_simul_seed : SIMUL_SEED EQUAL INT_NUMBER { driver.option_num("simul_seed", $3); } ;
 o_qz_criterium : QZ_CRITERIUM EQUAL number { driver.option_num("qz_criterium", $3); };
 o_datafile : DATAFILE EQUAL filename { driver.option_str("datafile", $3); };
diff --git a/preprocessor/ParsingDriver.cc b/preprocessor/ParsingDriver.cc
index e4b000827b5aa8144b7c0c221c5fd575afdf801d..890f5741492bad899aae50efacdf609516f8df46 100644
--- a/preprocessor/ParsingDriver.cc
+++ b/preprocessor/ParsingDriver.cc
@@ -648,9 +648,7 @@ ParsingDriver::option_num(const string &name_option, string *opt)
 void
 ParsingDriver::option_num(const string &name_option, const string &opt)
 {
-  // Since "periods" option automatically sets "simul" option, we don't want to fail if user explicitly sets both "simul" and "periods"
-  if (name_option != "simul"
-      && (options_list.num_options.find(name_option) != options_list.num_options.end()))
+  if (options_list.num_options.find(name_option) != options_list.num_options.end())
     error("option " + name_option + " declared twice");
 
   if ((name_option == "periods") && mod_file->block)