diff --git a/doc/manual/source/running-dynare.rst b/doc/manual/source/running-dynare.rst
index b7ca7d7dd7203454ae8753702b964c4428e0c5e5..d36c3fc88aad1a1fac6988cfa72645b062ece405 100644
--- a/doc/manual/source/running-dynare.rst
+++ b/doc/manual/source/running-dynare.rst
@@ -357,6 +357,13 @@ by the ``dynare`` command.
         without executing the ``.mod`` file. See :ref:`conf-file`, for
         more information about the configuration file.
 
+    .. option:: parallel_use_psexec=true|false
+
+        For local execution under Windows operating system,
+        set ``parallel_use_psexec=false`` to use ``start``
+        instead of ``psexec``, to properly allocate affinity when there are
+        more than 32 cores in the local machine. [default=true]
+
     .. option:: -DMACRO_VARIABLE=MACRO_EXPRESSION
 
         Defines a macro-variable from the command line (the same effect as
diff --git a/doc/manual/source/the-model-file.rst b/doc/manual/source/the-model-file.rst
index 48f5be41c42df9093d8907b99de580db2ded06fc..c9f432c5fe23c06783096af132fb4a6c344499bf 100644
--- a/doc/manual/source/the-model-file.rst
+++ b/doc/manual/source/the-model-file.rst
@@ -4416,10 +4416,15 @@ The coefficients of the decision rules are stored as follows:
   to all endogenous in the declaration order.
 * :math:`A` is stored in ``oo_.dr.ghx``. The matrix rows correspond to
   all endogenous in DR-order. The matrix columns correspond to state
-  variables in DR-order.
+  variables in DR-order, as given by ``oo_.dr.state_var``. (N.B.: if the
+  ``block`` option to the ``model`` block has been specified, then rows
+  are in declaration order, and columns are ordered
+  according to ``oo_.dr.state_var`` which may differ from DR-order.)
 * :math:`B` is stored ``oo_.dr.ghu``. The matrix rows correspond to
   all endogenous in DR-order. The matrix columns correspond to
-  exogenous variables in declaration order.
+  exogenous variables in declaration order.  (N.B.: if the
+  ``block`` option to the ``model`` block has been specified, then rows
+  are in declaration order.)
 
 Of course, the shown form of the approximation is only stylized,
 because it neglects the required different ordering in :math:`y^s` and
@@ -5111,6 +5116,16 @@ block decomposition of the model (see :opt:`block`).
 
        See :opt:`graph_format <graph_format = ( FORMAT, FORMAT... )>`.
 
+    .. option:: no_init_estimation_check_first_obs
+
+       Do not check for stochastic singularity in first period. If used, `ESTIMATION CHECKS`
+       does not return an error if the check fails only in first observation.
+       This should only be used when observing stock variables (e.g. capital) in first period, on top of their associated flow (e.g. investment).
+       Using this option may lead to a crash or provide undesired/wrong results for badly specified problems 
+       (e.g. the additional variable observed in first period is not predetermined).
+
+       For advanced use only.
+
     .. option:: lik_init = INTEGER
 
        Type of initialization of Kalman filter:
@@ -5542,6 +5557,48 @@ block decomposition of the model (see :opt:`block`).
        density, and posterior statistics from an existing ``_results``
        file instead of recomputing them.
 
+    .. option:: mh_initialize_from_previous_mcmc
+
+       This option allows to pick initial values for new MCMC from a previous one, 
+       where the model specification, the number of estimated parameters, 
+       (some) prior might have changed (so a situation where ``load_mh_file`` would not work).
+       If an additional parameter is estimated, it is automatically initialized from prior_draw. 
+       Note that, if this option is used to skip the optimization step, you should use a sampling method which does not require
+       a proposal density, like slice. Otherwise, optimization should always be done beforehand or a mode file with
+       an appropriate posterior covariance matrix should be used.
+
+    .. option:: mh_initialize_from_previous_mcmc_directory = FILENAME
+
+       If ``mh_initialize_from_previous_mcmc`` is set, users must provide here 
+       the path to the standard FNAME folder from where to load prior definitions and
+       last MCMC values to be used to initialize the new MCMC. 
+
+       Example: if previous project directory is ``/my_previous_dir`` and FNAME is ``mymodel``, 
+       users should set the option as 
+
+       ``mh_initialize_from_previous_mcmc_directory = '/my_previous_dir/mymodel'``
+
+       Dynare will then look for the last record file into 
+
+       ``/my_previous_dir/mymodel/metropolis/mymodel_mh_history_<LAST>.mat``
+
+       and for the prior definition file into 
+
+       ``/my_previous_dir/mymodel/prior/definition.mat``
+
+    .. option:: mh_initialize_from_previous_mcmc_record = FILENAME
+
+       If ``mh_initialize_from_previous_mcmc`` is set, and whenever the standard file or directory tree
+       is not applicable to load initial values, users may directly provide here 
+       the path to the record file from which to load
+       values to be used to initialize the new MCMC.
+
+    .. option:: mh_initialize_from_previous_mcmc_prior = FILENAME
+
+       If ``mh_initialize_from_previous_mcmc`` is set, and whenever the standard file or directory tree
+       is not applicable to load initial values, users may directly provide here 
+       the path to the prior definition file, to get info in the priors used in previous MCMC.
+
     .. option:: optim = (NAME, VALUE, ...)
 
        A list of NAME and VALUE pairs. Can be used to set options for
@@ -7821,11 +7878,11 @@ Shock Decomposition
         decomposition. Default: ``0``.
 
     .. option:: fast_realtime = INTEGER
+                fast_realtime = [INTEGER1:INTEGER2]
+                fast_realtime = [INTEGER1 INTEGER2 ...]
 
-        Runs the smoother only twice: once for the last in-sample and
-        once for the last out-of-sample data point, where the provided
-        integer defines the last observation (equivalent to
-        :opt:`nobs`). Default: not enabled.
+        Runs the smoother only for the data vintages provided
+        by the specified integer (vector).
 
     .. option:: with_epilogue
 
diff --git a/license.txt b/license.txt
index 43f6280eb41fef21376223c3e61b027bea4e75c8..c45679977adeedf6d9db90138f1fbf3bb1c99750 100644
--- a/license.txt
+++ b/license.txt
@@ -246,6 +246,7 @@ License: GPL-3+
 
 Files: m4/ax_blas.m4 m4/ax_lapack.m4
 Copyright: 2008-2009 Steven G. Johnson <stevenj@alum.mit.edu>
+           2019 Geoffrey M. Oxberry <goxberry@gmail.com>
 License: GPL-3+ with Autoconf exception
 
 Files: m4/ax_compare_version.m4
@@ -262,6 +263,7 @@ Copyright: 2008 Benjamin Kosnik <bkoz@redhat.com>
            2015 Moritz Klammler <moritz@klammler.eu>
            2016, 2018 Krzesimir Nowak <qdlacz@gmail.com>
            2019 Enji Cooper <yaneurabeya@gmail.com>
+           2020 Jason Merrill <jason@redhat.com>
 License: FSFAP
 
 Files: m4/ax_latex_class.m4 m4/ax_latex_test.m4
@@ -311,6 +313,7 @@ Copyright: 2008 Benjamin Kosnik <bkoz@redhat.com>
            2015 Moritz Klammler <moritz@klammler.eu>
            2016, 2018 Krzesimir Nowak <qdlacz@gmail.com>
            2019 Enji Cooper <yaneurabeya@gmail.com>
+           2020 Jason Merrill <jason@redhat.com>
 License: FSFAP
 
 Files: preprocessor/m4/ax_latex_class.m4
diff --git a/m4/ax_blas.m4 b/m4/ax_blas.m4
index bf9f29c628823cba855f919a9090edad2ee7cd92..95719050a97cfda2421f66fb8065f560e167d3e0 100644
--- a/m4/ax_blas.m4
+++ b/m4/ax_blas.m4
@@ -36,6 +36,7 @@
 # LICENSE
 #
 #   Copyright (c) 2008 Steven G. Johnson <stevenj@alum.mit.edu>
+#   Copyright (c) 2019 Geoffrey M. Oxberry <goxberry@gmail.com>
 #
 #   This program is free software: you can redistribute it and/or modify it
 #   under the terms of the GNU General Public License as published by the
@@ -63,11 +64,11 @@
 #   modified version of the Autoconf Macro, you may extend this special
 #   exception to the GPL to apply to your modified version as well.
 
-#serial 15
+#serial 17
 
 AU_ALIAS([ACX_BLAS], [AX_BLAS])
 AC_DEFUN([AX_BLAS], [
-AC_PREREQ(2.50)
+AC_PREREQ([2.55])
 AC_REQUIRE([AC_F77_LIBRARY_LDFLAGS])
 AC_REQUIRE([AC_CANONICAL_HOST])
 ax_blas_ok=no
@@ -77,7 +78,9 @@ AC_ARG_WITH(blas,
 case $with_blas in
 	yes | "") ;;
 	no) ax_blas_ok=disable ;;
-	-* | */* | *.a | *.so | *.so.* | *.o) BLAS_LIBS="$with_blas" ;;
+	-* | */* | *.a | *.so | *.so.* | *.dylib | *.dylib.* | *.o)
+		BLAS_LIBS="$with_blas"
+	;;
 	*) BLAS_LIBS="-l$with_blas" ;;
 esac
 
@@ -93,7 +96,7 @@ if test $ax_blas_ok = no; then
 if test "x$BLAS_LIBS" != x; then
 	save_LIBS="$LIBS"; LIBS="$BLAS_LIBS $LIBS"
 	AC_MSG_CHECKING([for $sgemm in $BLAS_LIBS])
-	AC_TRY_LINK_FUNC($sgemm, [ax_blas_ok=yes], [BLAS_LIBS=""])
+	AC_LINK_IFELSE([AC_LANG_CALL([], [$sgemm])], [ax_blas_ok=yes], [BLAS_LIBS=""])
 	AC_MSG_RESULT($ax_blas_ok)
 	LIBS="$save_LIBS"
 fi
@@ -103,7 +106,7 @@ fi
 if test $ax_blas_ok = no; then
 	save_LIBS="$LIBS"; LIBS="$LIBS"
 	AC_MSG_CHECKING([if $sgemm is being linked in already])
-	AC_TRY_LINK_FUNC($sgemm, [ax_blas_ok=yes])
+	AC_LINK_IFELSE([AC_LANG_CALL([], [$sgemm])], [ax_blas_ok=yes])
 	AC_MSG_RESULT($ax_blas_ok)
 	LIBS="$save_LIBS"
 fi
@@ -174,7 +177,7 @@ fi
 if test $ax_blas_ok = no; then
 	save_LIBS="$LIBS"; LIBS="-framework vecLib $LIBS"
 	AC_MSG_CHECKING([for $sgemm in -framework vecLib])
-	AC_TRY_LINK_FUNC($sgemm, [ax_blas_ok=yes;BLAS_LIBS="-framework vecLib"])
+	AC_LINK_IFELSE([AC_LANG_CALL([], [$sgemm])], [ax_blas_ok=yes;BLAS_LIBS="-framework vecLib"])
 	AC_MSG_RESULT($ax_blas_ok)
 	LIBS="$save_LIBS"
 fi
diff --git a/m4/ax_cxx_compile_stdcxx.m4 b/m4/ax_cxx_compile_stdcxx.m4
index 43087b2e6889ec6f8ebd2f8ba77f4a9a716f8ac2..9413da624d2545123501b7788b7ac6d96fd322e8 100644
--- a/m4/ax_cxx_compile_stdcxx.m4
+++ b/m4/ax_cxx_compile_stdcxx.m4
@@ -16,7 +16,7 @@
 #   The second argument, if specified, indicates whether you insist on an
 #   extended mode (e.g. -std=gnu++11) or a strict conformance mode (e.g.
 #   -std=c++11).  If neither is specified, you get whatever works, with
-#   preference for an extended mode.
+#   preference for no added switch, and then for an extended mode.
 #
 #   The third argument, if specified 'mandatory' or if left unspecified,
 #   indicates that baseline support for the specified C++ standard is
@@ -35,13 +35,14 @@
 #   Copyright (c) 2015 Moritz Klammler <moritz@klammler.eu>
 #   Copyright (c) 2016, 2018 Krzesimir Nowak <qdlacz@gmail.com>
 #   Copyright (c) 2019 Enji Cooper <yaneurabeya@gmail.com>
+#   Copyright (c) 2020 Jason Merrill <jason@redhat.com>
 #
 #   Copying and distribution of this file, with or without modification, are
 #   permitted in any medium without royalty provided the copyright notice
 #   and this notice are preserved.  This file is offered as-is, without any
 #   warranty.
 
-#serial 11
+#serial 12
 
 dnl  This macro is based on the code from the AX_CXX_COMPILE_STDCXX_11 macro
 dnl  (serial version number 13).
@@ -62,6 +63,16 @@ AC_DEFUN([AX_CXX_COMPILE_STDCXX], [dnl
   AC_LANG_PUSH([C++])dnl
   ac_success=no
 
+  m4_if([$2], [], [dnl
+    AC_CACHE_CHECK(whether $CXX supports C++$1 features by default,
+		   ax_cv_cxx_compile_cxx$1,
+      [AC_COMPILE_IFELSE([AC_LANG_SOURCE([_AX_CXX_COMPILE_STDCXX_testbody_$1])],
+        [ax_cv_cxx_compile_cxx$1=yes],
+        [ax_cv_cxx_compile_cxx$1=no])])
+    if test x$ax_cv_cxx_compile_cxx$1 = xyes; then
+      ac_success=yes
+    fi])
+
   m4_if([$2], [noext], [], [dnl
   if test x$ac_success = xno; then
     for alternative in ${ax_cxx_compile_alternatives}; do
diff --git a/m4/ax_lapack.m4 b/m4/ax_lapack.m4
index e7eadd17b7d2de0eca19f921b6f6f3304a07f4a0..abaff9d1714ab3158d6235188eb8c01b6e69a0fa 100644
--- a/m4/ax_lapack.m4
+++ b/m4/ax_lapack.m4
@@ -37,6 +37,7 @@
 # LICENSE
 #
 #   Copyright (c) 2009 Steven G. Johnson <stevenj@alum.mit.edu>
+#   Copyright (c) 2019 Geoffrey M. Oxberry <goxberry@gmail.com>
 #
 #   This program is free software: you can redistribute it and/or modify it
 #   under the terms of the GNU General Public License as published by the
@@ -64,7 +65,7 @@
 #   modified version of the Autoconf Macro, you may extend this special
 #   exception to the GPL to apply to your modified version as well.
 
-#serial 8
+#serial 10
 
 AU_ALIAS([ACX_LAPACK], [AX_LAPACK])
 AC_DEFUN([AX_LAPACK], [
@@ -76,7 +77,9 @@ AC_ARG_WITH(lapack,
 case $with_lapack in
         yes | "") ;;
         no) ax_lapack_ok=disable ;;
-        -* | */* | *.a | *.so | *.so.* | *.o) LAPACK_LIBS="$with_lapack" ;;
+        -* | */* | *.a | *.so | *.so.* | *.dylib | *.dylib.* | *.o)
+                 LAPACK_LIBS="$with_lapack"
+        ;;
         *) LAPACK_LIBS="-l$with_lapack" ;;
 esac
 
@@ -93,7 +96,7 @@ fi
 if test "x$LAPACK_LIBS" != x; then
         save_LIBS="$LIBS"; LIBS="$LAPACK_LIBS $BLAS_LIBS $LIBS $FLIBS"
         AC_MSG_CHECKING([for $cheev in $LAPACK_LIBS])
-        AC_TRY_LINK_FUNC($cheev, [ax_lapack_ok=yes], [LAPACK_LIBS=""])
+        AC_LINK_IFELSE([AC_LANG_CALL([], [$cheev])], [ax_lapack_ok=yes], [LAPACK_LIBS=""])
         AC_MSG_RESULT($ax_lapack_ok)
         LIBS="$save_LIBS"
         if test $ax_lapack_ok = no; then
diff --git a/m4/ax_mexopts.m4 b/m4/ax_mexopts.m4
index eff7dc0c905bca5d65a738def4883795cc5bc7be..51fbfe816c177836bc4d257724c1b3141aab4cad 100644
--- a/m4/ax_mexopts.m4
+++ b/m4/ax_mexopts.m4
@@ -29,35 +29,27 @@ AC_MSG_CHECKING([for options to compile MEX for MATLAB])
 MATLAB_CPPFLAGS="-I$MATLAB/extern/include"
 
 case ${MATLAB_ARCH} in
-  glnx86 | glnxa64)
-    MATLAB_DEFS="$MATLAB_DEFS -D_GNU_SOURCE -DNDEBUG"
-    MATLAB_CFLAGS="-fexceptions -fPIC -pthread -g -O2"
-    MATLAB_CXXFLAGS="-fPIC -pthread -g -O2"
-    MATLAB_FCFLAGS="-fPIC -g -O2 -fexceptions"
-    MATLAB_LDFLAGS_NOMAP="-shared -Wl,--no-undefined -Wl,-rpath-link,$MATLAB/bin/${MATLAB_ARCH} -L$MATLAB/bin/${MATLAB_ARCH}"
-    MATLAB_LDFLAGS="$MATLAB_LDFLAGS_NOMAP -Wl,--version-script,$MATLAB/extern/lib/${MATLAB_ARCH}/mexFunction.map"
-    MATLAB_LIBS="-lmx -lmex -lmat -lm -lstdc++ -lmwlapack -lmwblas"
-    if test "${MATLAB_ARCH}" = "glnx86"; then
-      MATLAB_DEFS="$MATLAB_DEFS -D_FILE_OFFSET_BITS=64"
-      MATLAB_CFLAGS="$MATLAB_CFLAGS -m32"
-      MATLAB_CXXFLAGS="$MATLAB_CXXFLAGS -m32"
-    else # glnxa64
-      MATLAB_CFLAGS="$MATLAB_CFLAGS -fno-omit-frame-pointer"
-      MATLAB_CXXFLAGS="$MATLAB_CXXFLAGS -fno-omit-frame-pointer"
-    fi
+  glnxa64)
+    MATLAB_DEFS="-D_GNU_SOURCE -DNDEBUG"
+    MATLAB_CFLAGS="-fexceptions -fPIC -pthread"
+    MATLAB_CXXFLAGS="-fPIC -pthread"
+    MATLAB_FCFLAGS="-fPIC -fexceptions"
+    MATLAB_LDFLAGS_NOMAP="-shared -Wl,--no-undefined -Wl,-rpath-link,$MATLAB/bin/glnxa64 -L$MATLAB/bin/glnxa64"
+    MATLAB_LDFLAGS="$MATLAB_LDFLAGS_NOMAP -Wl,--version-script,$MATLAB/extern/lib/glnxa64/mexFunction.map"
+    MATLAB_LIBS="-lmx -lmex -lmat -lm -lmwlapack -lmwblas"
     ax_mexopts_ok="yes"
     ;;
-  win32 | win64)
-    MATLAB_CFLAGS="-fexceptions -g -O2"
-    MATLAB_CXXFLAGS="-g -O2"
-    MATLAB_FCFLAGS="-g -O2 -fexceptions -fno-underscoring"
-    MATLAB_DEFS="$MATLAB_DEFS -DNDEBUG"
+  win64)
+    MATLAB_CFLAGS="-fexceptions"
+    MATLAB_CXXFLAGS=""
+    MATLAB_FCFLAGS="-fexceptions -fno-underscoring"
+    MATLAB_DEFS="-DNDEBUG"
     # The hack for libquadmath is needed because -static-libgfortran
     # unfortunately does not imply the static linking of the former.
     # The last part about winpthread is a hack to avoid dynamically linking
     # against libwinpthread DLL (which is pulled in by libstdc++, even without
     # using threads, since we are using the POSIX threads version of MinGW)
-    MATLAB_LDFLAGS_NOMAP="-static-libgcc -static-libstdc++ -static-libgfortran -Wl,-Bstatic,--whole-archive -lquadmath -Wl,-Bdynamic,--no-whole-archive -shared -L$MATLAB/bin/${MATLAB_ARCH} -Wl,-Bstatic,--whole-archive -lwinpthread -Wl,-Bdynamic,--no-whole-archive"
+    MATLAB_LDFLAGS_NOMAP="-static-libgcc -static-libstdc++ -static-libgfortran -Wl,-Bstatic,--whole-archive -lquadmath -Wl,-Bdynamic,--no-whole-archive -shared -L$MATLAB/bin/win64 -Wl,-Bstatic,--whole-archive -lwinpthread -Wl,-Bdynamic,--no-whole-archive"
     MATLAB_LDFLAGS="$MATLAB_LDFLAGS_NOMAP \$(abs_top_srcdir)/mex.def"
     MATLAB_LIBS="-lmex -lmx -lmat -lmwlapack -lmwblas"
     # Hack for static linking of libgomp, needed for OpenMP
@@ -65,10 +57,10 @@ case ${MATLAB_ARCH} in
     ax_mexopts_ok="yes"
     ;;
   maci64)
-    MATLAB_DEFS="$MATLAB_DEFS -DNDEBUG"
+    MATLAB_DEFS="-DNDEBUG"
     MATLAB_CFLAGS="-fno-common -fexceptions"
-    MATLAB_CXXFLAGS="$MATLAB_CFLAGS"
-    MATLAB_FCFLAGS="-g -O2 -fexceptions -fbackslash"
+    MATLAB_CXXFLAGS="-fno-common -fexceptions"
+    MATLAB_FCFLAGS="-fexceptions -fbackslash"
     MATLAB_LDFLAGS_NOMAP="-Wl,-twolevel_namespace -undefined error -bundle"
     MATLAB_LDFLAGS="$MATLAB_LDFLAGS_NOMAP -Wl,-exported_symbols_list,\$(abs_top_srcdir)/mexFunction-MacOSX.map"
     # This -L flag is put here, hence later on the linker command line, so as
diff --git a/m4/ax_slicot.m4 b/m4/ax_slicot.m4
index 37bfaa843f4afec2e09901d582d198738551ad01..c3bc018e883014c45ca0827106ef8a3d19cc92c1 100644
--- a/m4/ax_slicot.m4
+++ b/m4/ax_slicot.m4
@@ -5,7 +5,7 @@ dnl
 dnl AX_SLICOT([matlab])
 dnl AX_SLICOT([octave])
 dnl
-dnl Copyright © 2012-2019 Dynare Team
+dnl Copyright © 2012-2021 Dynare Team
 dnl
 dnl This file is part of Dynare.
 dnl
@@ -37,28 +37,17 @@ AC_DEFUN([AX_SLICOT],
   else
     LDFLAGS_SLICOT=""
   fi
-  ac_save_LDFLAGS="$LDFLAGS"
-  LDFLAGS_SAVED="$LDFLAGS"
+  ac_save_LDFLAGS=$LDFLAGS
 
+  # At this point we should add MATLAB_FCFLAGS to FCFLAGS for Windows (which has -fno-underscoring),
+  # but that does not work. The actual underscore test seems to happen at the very beginning of the
+  # macro. Hence the modification of FCFLAGS was moved higher (in mex/build/matlab/configure.ac).
   AC_FC_FUNC(sb02od)
 
   if test "$1" = matlab; then
     LDFLAGS="$LDFLAGS $MATLAB_LDFLAGS_NOMAP $LDFLAGS_SLICOT"
 
-    case ${MATLAB_ARCH} in
-       glnxa64 | win64 | maci64)
-         AX_COMPARE_VERSION([$MATLAB_VERSION], [ge], [7.8], [use_64_bit_indexing=yes], [use_64_bit_indexing=no])
-         ;;
-       *)
-         use_64_bit_indexing=no
-         ;;
-    esac
-
-    if test "$use_64_bit_indexing" = yes; then
-       AC_CHECK_LIB([slicot64_pic], [$sb02od], [LIBADD_SLICOT="-lslicot64_pic"], [has_slicot=no], [$MATLAB_LIBS])
-    else
-       AC_CHECK_LIB([slicot_pic], [$sb02od], [LIBADD_SLICOT="-lslicot_pic"], [has_slicot=no], [$MATLAB_LIBS])
-    fi
+    AC_CHECK_LIB([slicot64_pic], [$sb02od], [LIBADD_SLICOT="-lslicot64_pic"], [has_slicot=no], [$MATLAB_LIBS])
   else
     LDFLAGS="$LDFLAGS $($MKOCTFILE -p LFLAGS) $LDFLAGS_SLICOT"
     # Fallback on libslicot_pic if dynamic libslicot not found
@@ -69,7 +58,7 @@ AC_DEFUN([AX_SLICOT],
              [$($MKOCTFILE -p BLAS_LIBS) $($MKOCTFILE -p LAPACK_LIBS)])
   fi
 
-  LDFLAGS="$ac_save_LDFLAGS"
+  LDFLAGS=$ac_save_LDFLAGS
   AC_SUBST(LDFLAGS_SLICOT)
   AC_SUBST(LIBADD_SLICOT)
 ])
diff --git a/matlab/WriteShockDecomp2Excel.m b/matlab/WriteShockDecomp2Excel.m
index f339f294e223ee90f3743389cb6c0f7c4565980c..8be1b6f8816c6f8bce4e5e9322d1ac25834e3fb6 100644
--- a/matlab/WriteShockDecomp2Excel.m
+++ b/matlab/WriteShockDecomp2Excel.m
@@ -11,7 +11,7 @@ function WriteShockDecomp2Excel(z,shock_names,endo_names,i_var,initial_date,Dyna
 %   DynareModel     [structure]                     Dynare model structure
 %   DynareOptions   [structure]                     Dynare options structure
 
-% Copyright (C) 2016-2018 Dynare Team
+% Copyright (C) 2016-2021 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -120,7 +120,7 @@ for j=1:nvar
     fig_name1 = strrep(fig_name1,'.','');
     
     if ~ismac
-        [STATUS,MESSAGE] = xlswrite([OutputDirectoryName,filesep,DynareModel.fname,'_shock_decomposition',fig_mode,fig_name1],d0,endo_names{i_var(j)});
+        STATUS = xlswrite([OutputDirectoryName,filesep,DynareModel.fname,'_shock_decomposition',fig_mode,fig_name1],d0,endo_names{i_var(j)});
     else
         writetable(cell2table(d0), [OutputDirectoryName,filesep,DynareModel.fname,'_shock_decomposition',fig_mode,fig_name1 '.xls'], 'Sheet', endo_names{i_var(j)},'WriteVariableNames',false);
     end
diff --git a/matlab/default_option_values.m b/matlab/default_option_values.m
index 1a6b3883edf9e1b518b91633dad6170a9cbcf3a4..dca273c2f9362d27abf97b6decbd7e8966a03fee 100644
--- a/matlab/default_option_values.m
+++ b/matlab/default_option_values.m
@@ -346,6 +346,7 @@ options_.ramsey.maxit = 500;
 
 % estimation
 options_.initial_period = NaN; %dates(1,1);
+options_.no_init_estimation_check_first_obs=false;
 options_.dataset.file = [];
 options_.dataset.series = [];
 options_.dataset.firstobs = dates();
@@ -395,6 +396,10 @@ options_.mh_tune_jscale.c1 = .02;
 options_.mh_tune_jscale.c2 = .05;
 options_.mh_tune_jscale.c3 = 4;
 options_.mh_init_scale = 2*options_.mh_jscale;
+options_.mh_initialize_from_previous_mcmc.status = false;
+options_.mh_initialize_from_previous_mcmc.directory = '';
+options_.mh_initialize_from_previous_mcmc.record = '';
+options_.mh_initialize_from_previous_mcmc.prior = '';
 options_.mh_mode = 1;
 options_.mh_nblck = 2;
 options_.mh_recover = false;
@@ -465,6 +470,7 @@ options_.parallel = 0;
 options_.parallel_info.isHybridMatlabOctave = false;
 options_.parallel_info.leaveSlaveOpen = 0;
 options_.parallel_info.RemoteTmpFolder = '';
+options_.parallel_info.use_psexec = true;
 options_.number_of_grid_points_for_kde = 2^9;
 quarter = 1;
 years = [1 2 3 4 5 10 20 30 40 50];
diff --git a/matlab/dr_block.m b/matlab/dr_block.m
index bb51b72610bbb94f35760680f8d37d2a65c17669..0f6069bb53adf635f61544bc3b88d186adb4dbb5 100644
--- a/matlab/dr_block.m
+++ b/matlab/dr_block.m
@@ -1,8 +1,11 @@
 function [dr,info,M_,options_,oo_] = dr_block(dr,task,M_,options_,oo_,varargin)
 % function [dr,info,M_,options_,oo_] = dr_block(dr,task,M_,options_,oo_,varargin)
-% computes the reduced form solution of a rational expectations model
+% computes the reduced form solution of a rational expectations block-decomposed model
 % (first order approximation of the stochastic model around the deterministic steady state).
 %
+% NB: This code does not work with option mfs > 0. The preprocessor has a check to avoid this
+%     configuration. See also #1726.
+%
 % INPUTS
 %   dr         [matlab structure] Decision rules for stochastic simulations.
 %   task       [integer]          if task = 0 then dr_block computes decision rules.
@@ -34,7 +37,7 @@ function [dr,info,M_,options_,oo_] = dr_block(dr,task,M_,options_,oo_,varargin)
 %   none.
 %
 
-% Copyright (C) 2010-2020 Dynare Team
+% Copyright (C) 2010-2021 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -63,13 +66,12 @@ end
 
 z = repmat(dr.ys,1,M_.maximum_lead + M_.maximum_lag + 1);
 zx = repmat([oo_.exo_simul oo_.exo_det_simul],M_.maximum_lead + M_.maximum_lag + 1, 1);
-if (isfield(M_,'block_structure'))
-    data = M_.block_structure.block;
-    Size = length(M_.block_structure.block);
-else
-    data = M_;
-    Size = 1;
+
+if ~isfield(M_,'block_structure')
+    error('Option ''block'' has not been specified')
 end
+data = M_.block_structure.block;
+
 if options_.bytecode
     [~, data]= bytecode('dynamic','evaluate', z, zx, M_.params, dr.ys, 1, data);
 else
@@ -93,7 +95,7 @@ n_sv = size(dr.state_var, 2);
 dr.ghx = zeros(M_.endo_nbr, length(dr.state_var));
 dr.exo_var = 1:M_.exo_nbr;
 dr.ghu = zeros(M_.endo_nbr, M_.exo_nbr);
-for i = 1:Size
+for i = 1:length(data)
     ghx = [];
     indexi_0 = 0;
     if (verbose)
diff --git a/matlab/dynare.m b/matlab/dynare.m
index 57bff4d1841ede9b31682260037fd3ed8fa9dd4a..65be458847206bf2ff0edc7181c47a82906fb914 100644
--- a/matlab/dynare.m
+++ b/matlab/dynare.m
@@ -280,8 +280,13 @@ clear(['+' fname '/driver'])
 try
     evalin('base',[fname '.driver']) ;
 catch ME
+    W = evalin('base','whos');
     diary off
-    rethrow(ME)
+    if ismember(fname,{W(:).name})
+        error('Your base workspace already contains a variable with the same name as the mod-file. You need to delete it or rename the mod-file.')        
+    else
+        rethrow(ME)    
+    end    
 end
 
 diary off
diff --git a/matlab/forcst.m b/matlab/forcst.m
index 4aa28230b6420cfb10a037f6c19d633ffd23fbab..ffa992cb350374c730a5aba69065d8d1c66d86bd 100644
--- a/matlab/forcst.m
+++ b/matlab/forcst.m
@@ -85,7 +85,7 @@ for i = 1:horizon
 end
 if nargout==3
     var_yf_ME=var_yf;
-    [loc_H, loc_varlist] = ismember(options_.varobs', options_.varlist);
+    [loc_H, loc_varlist] = ismember(options_.varobs, var_list);
     loc_varlist(loc_varlist==0) = [];
     if ~isempty(loc_varlist)
         var_yf_ME(:,loc_varlist) = var_yf(:,loc_varlist)+repmat(diag(M_.H(loc_H,loc_H))', horizon, 1);
diff --git a/matlab/initial_estimation_checks.m b/matlab/initial_estimation_checks.m
index e9100e4644b94687d0bea27f81d966ceff15ec77..9f6618282f2f67ba779fb24f49bb3962b37ba13b 100644
--- a/matlab/initial_estimation_checks.m
+++ b/matlab/initial_estimation_checks.m
@@ -39,7 +39,8 @@ function DynareResults = initial_estimation_checks(objective_function,xparam1,Dy
 
 %get maximum number of simultaneously observed variables for stochastic
 %singularity check
-maximum_number_non_missing_observations=max(sum(~isnan(DynareDataset.data),2));
+maximum_number_non_missing_observations=max(sum(~isnan(DynareDataset.data(2:end,:)),2));
+init_number_non_missing_observations=sum(~isnan(DynareDataset.data(1,:)),2);
 
 if DynareOptions.order>1
     if any(any(isnan(DynareDataset.data)))
@@ -82,13 +83,33 @@ end
 
 non_zero_ME=length(EstimatedParameters.H_entries_to_check_for_positive_definiteness);
 
+print_init_check_warning=false;
 if maximum_number_non_missing_observations>Model.exo_nbr+non_zero_ME
     error(['initial_estimation_checks:: Estimation can''t take place because there are less declared shocks than observed variables!'])
 end
+if init_number_non_missing_observations>Model.exo_nbr+non_zero_ME
+    if DynareOptions.no_init_estimation_check_first_obs
+        print_init_check_warning=true;
+    else
+        error(['initial_estimation_checks:: Estimation can''t take place because there are less declared shocks than observed variables in first period!'])
+    end
+end
 
 if maximum_number_non_missing_observations>length(find(diag(Model.Sigma_e)))+non_zero_ME
     error(['initial_estimation_checks:: Estimation can''t take place because too many shocks have been calibrated with a zero variance!'])
 end
+if init_number_non_missing_observations>length(find(diag(Model.Sigma_e)))+non_zero_ME
+    if DynareOptions.no_init_estimation_check_first_obs
+        print_init_check_warning=true;
+    else
+        error(['initial_estimation_checks:: Estimation can''t take place because too many shocks have been calibrated with a zero variance in first period!'])
+    end
+end
+if print_init_check_warning
+    fprintf('ESTIMATION_CHECKS: You decided to ignore test of stochastic singularity in first_obs.\n');
+    fprintf('ESTIMATION_CHECKS: If this was not done on purpose (typically when observing a stock variable [capital] in first period, on top of its flow [investment]),\n');
+    fprintf('ESTIMATION_CHECKS: it may lead to a crash or provide undesired/wrong results later on!\n');
+end
 
 if (any(BayesInfo.pshape  >0 ) && DynareOptions.mh_replic) && DynareOptions.mh_nblck<1
     error(['initial_estimation_checks:: Bayesian estimation cannot be conducted with mh_nblocks=0.'])
diff --git a/matlab/missing_DiffuseKalmanSmootherH1_Z.m b/matlab/missing_DiffuseKalmanSmootherH1_Z.m
index 2a4ed39c4ad778f9eb8a8812c2a37ed989cc5aa5..fbdf7e21173e4a70d4796c84c578cee04bdbe658 100644
--- a/matlab/missing_DiffuseKalmanSmootherH1_Z.m
+++ b/matlab/missing_DiffuseKalmanSmootherH1_Z.m
@@ -211,10 +211,12 @@ while notsteady && t<smpl
         P(:,:,t+1)  = T*P(:,:,t)*L(:,:,t)' + QQ;
     end
     a(:,t+1)    = T*atilde(:,t);
-    Pf          = P(:,:,t);
+    Pf          = P(:,:,t+1);
     aK(1,:,t+1) = a(:,t+1);
     for jnk=1:nk
-        Pf = T*Pf*T' + QQ;
+        if jnk>1
+            Pf = T*Pf*T' + QQ;
+        end
         PK(jnk,:,:,t+jnk) = Pf;
         if jnk>1
             aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
diff --git a/matlab/missing_DiffuseKalmanSmootherH3_Z.m b/matlab/missing_DiffuseKalmanSmootherH3_Z.m
index cfdbdd08c35c5a23d4ee39ebb83522e0695ffc78..ff30820ef16124c25608de342abc86573fcbc1e5 100644
--- a/matlab/missing_DiffuseKalmanSmootherH3_Z.m
+++ b/matlab/missing_DiffuseKalmanSmootherH3_Z.m
@@ -1,5 +1,5 @@
-function [alphahat,epsilonhat,etahat,a,P,aK,PK,decomp,V] = missing_DiffuseKalmanSmootherH3_Z(a_initial,T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag,state_uncertainty_flag)
-% function [alphahat,epsilonhat,etahat,a1,P,aK,PK,d,decomp] = missing_DiffuseKalmanSmootherH3_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,decomp_flag,state_uncertainty_flag)
+function [alphahat,epsilonhat,etahat,a,P1,aK,PK,decomp,V] = missing_DiffuseKalmanSmootherH3_Z(a_initial,T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,diffuse_kalman_tol,decomp_flag,state_uncertainty_flag)
+% function [alphahat,epsilonhat,etahat,a1,P1,aK,PK,d,decomp] = missing_DiffuseKalmanSmootherH3_Z(T,Z,R,Q,H,Pinf1,Pstar1,Y,pp,mm,smpl,data_index,nk,kalman_tol,decomp_flag,state_uncertainty_flag)
 % Computes the diffuse kalman smoother in the case of a singular var-cov matrix.
 % Univariate treatment of multivariate time series.
 %
@@ -31,7 +31,7 @@ function [alphahat,epsilonhat,etahat,a,P,aK,PK,decomp,V] = missing_DiffuseKalman
 %    a:        matrix of updated variables (a_{t|t})
 %    aK:       3D array of k step ahead filtered state variables (a_{t+k|t})
 %              (meaningless for periods 1:d)
-%    P:        3D array of one-step ahead forecast error variance
+%    P1:        3D array of one-step ahead forecast error variance
 %              matrices
 %    PK:       4D array of k-step ahead forecast error variance
 %              matrices (meaningless for periods 1:d)
@@ -221,18 +221,21 @@ while notsteady && t<smpl
         end
     end
     a1(:,t+1) = T*a(:,t);                                                   %transition according to (6.14) in DK (2012)
-    Pf          = P(:,:,t);
+    P(:,:,t+1) = T*P(:,:,t)*T' + QQ;                                        %transition according to (6.14) in DK (2012)
+    Pf          = P(:,:,t+1);
     aK(1,:,t+1) = a1(:,t+1);
     for jnk=1:nk
-        Pf = T*Pf*T' + QQ;
+        if jnk>1
+            Pf = T*Pf*T' + QQ;
+        end
         PK(jnk,:,:,t+jnk) = Pf;
         if jnk>1
             aK(jnk,:,t+jnk) = T*dynare_squeeze(aK(jnk-1,:,t+jnk-1));
         end
     end
-    P(:,:,t+1) = T*P(:,:,t)*T' + QQ;                                        %transition according to (6.14) in DK (2012)
                                                                             %  notsteady   = ~(max(max(abs(P(:,:,t+1)-P(:,:,t))))<kalman_tol);
 end
+P1(:,:,t+1) = P(:,:,t+1);
 % $$$ P_s=tril(P(:,:,t))+tril(P(:,:,t),-1)';
 % $$$ P1_s=tril(P1(:,:,t))+tril(P1(:,:,t),-1)';
 % $$$ Fi_s = Fi(:,t);
diff --git a/matlab/modules/dseries b/matlab/modules/dseries
index a02a4ce6745c1e23ff9aaedd4db651c1acf3cd47..bd31aec3782b379d7197848643cd41e28257c429 160000
--- a/matlab/modules/dseries
+++ b/matlab/modules/dseries
@@ -1 +1 @@
-Subproject commit a02a4ce6745c1e23ff9aaedd4db651c1acf3cd47
+Subproject commit bd31aec3782b379d7197848643cd41e28257c429
diff --git a/matlab/ols/pooled_ols.m b/matlab/ols/pooled_ols.m
index e47d77646265136e9b30b55de4dd1b63f17189e3..df51efd2577ab681ead1dcf92aba1c5c3462a025 100644
--- a/matlab/ols/pooled_ols.m
+++ b/matlab/ols/pooled_ols.m
@@ -76,7 +76,11 @@ else
 end
 
 st = dbstack(1);
-if strcmp(st(1).name, 'pooled_fgls')
+if isoctave
+    % Workaround for bug in Octave 6, see https://savannah.gnu.org/bugs/?60531
+    st = st(2:end);
+end
+if ~isempty(st) && strcmp(st(1).name, 'pooled_fgls')
     save_structure_name = 'pooled_fgls';
 else
     save_structure_name = 'pooled_ols';
diff --git a/matlab/ols/sur.m b/matlab/ols/sur.m
index 618e4cbe69b4c14318571b5beaafad07d30b6ae4..8002d13bf5b2ba0586ea9264f928ac8b725cfbd3 100644
--- a/matlab/ols/sur.m
+++ b/matlab/ols/sur.m
@@ -121,6 +121,10 @@ end
 %
 
 st = dbstack(1);
+if isoctave
+    % Workaround for bug in Octave 6, see https://savannah.gnu.org/bugs/?60531
+    st = st(2:end);
+end
 if ~isempty(st) && strcmp(st(1).name, 'surgibbs')
     varargout{1} = nobs;
     varargout{2} = X{param_names{:}}.data;
diff --git a/matlab/optimization/gmhmaxlik.m b/matlab/optimization/gmhmaxlik.m
index 4f1db628d06cf6070678076d47a8d36daa8dbe08..339526efa6458bc5884d0325ed793f6221468276 100644
--- a/matlab/optimization/gmhmaxlik.m
+++ b/matlab/optimization/gmhmaxlik.m
@@ -105,6 +105,7 @@ for i=1:gmhmaxlikOptions.iterations
     printline(58,'=')
     disp(['   Change in the posterior covariance matrix = ' num2str(dVariance) '.'])
     disp(['   Change in the posterior mean = ' num2str(dMean) '.'])
+    disp(['   Current mode = ' num2str(ModeValue)])
     disp(['   Mode improvement = ' num2str(abs(OldModeValue-ModeValue))])
     disp(['   New value of jscale = ' num2str(Scale)])
     printline(58,'=')
diff --git a/matlab/parallel/masterParallel.m b/matlab/parallel/masterParallel.m
index 1c102ecf47860a82e3f53de5ee6b3c6ea89fe1ce..89581db984a0e549782db7d643503cf139b9cdf1 100644
--- a/matlab/parallel/masterParallel.m
+++ b/matlab/parallel/masterParallel.m
@@ -353,10 +353,21 @@ if parallel_recover ==0
                         command1=[Parallel(indPC).MatlabOctavePath,' -nosplash -nodesktop -minimize ',compThread,' -r "addpath(''',Parallel(indPC).DynarePath,'''), dynareroot = dynare_config(); fParallel(',int2str(offset+1),',',int2str(sum(nBlockPerCPU(1:j))),',',int2str(j),',',int2str(indPC),',''',fname,''')" &'];
                     end
                 else    % Hybrid computing Matlab(Master)->Octave(Slaves) and Vice Versa!
+                    if Parallel_info.use_psexec
+                        token1 = ['psexec -accepteula -d -W "',DyMo, '" -a ',my_affinity,' -low  '];
+                    else
+                        itmp = intmin('uint64');
+                        cpus = eval([ '['  my_affinity ']' ])+1;
+                        for icpu=1:length(cpus)
+                            itmp = bitset(itmp,cpus(icpu));
+                        end
+                        hex_affinity = dec2hex(itmp);
+                        token1 = ['start /B /D "',DyMo, '" /affinity ',hex_affinity,' /LOW  '];
+                    end
                     if  regexpi([Parallel(indPC).MatlabOctavePath], 'octave')
-                        command1=['psexec -accepteula -d -W "',DyMo, '" -a ',my_affinity,' -low  ',Parallel(indPC).MatlabOctavePath,' -f --eval "default_save_options(''-v7''); addpath(''',Parallel(indPC).DynarePath,'''), dynareroot = dynare_config(); fParallel(',int2str(offset+1),',',int2str(sum(nBlockPerCPU(1:j))),',',int2str(j),',',int2str(indPC),',''',fname,''')"'];
+                        command1=[token1,Parallel(indPC).MatlabOctavePath,' -f --eval "default_save_options(''-v7''); addpath(''',Parallel(indPC).DynarePath,'''), dynareroot = dynare_config(); fParallel(',int2str(offset+1),',',int2str(sum(nBlockPerCPU(1:j))),',',int2str(j),',',int2str(indPC),',''',fname,''')"'];
                     else
-                        command1=['psexec -accepteula -d -W "',DyMo, '" -a ',my_affinity,' -low  ',Parallel(indPC).MatlabOctavePath,' -nosplash -nodesktop -minimize ',compThread,' -r "addpath(''',Parallel(indPC).DynarePath,'''), dynareroot = dynare_config(); fParallel(',int2str(offset+1),',',int2str(sum(nBlockPerCPU(1:j))),',',int2str(j),',',int2str(indPC),',''',fname,''')"'];
+                        command1=[token1,Parallel(indPC).MatlabOctavePath,' -nosplash -nodesktop -minimize ',compThread,' -r "addpath(''',Parallel(indPC).DynarePath,'''), dynareroot = dynare_config(); fParallel(',int2str(offset+1),',',int2str(sum(nBlockPerCPU(1:j))),',',int2str(j),',',int2str(indPC),',''',fname,''')"'];
                     end
                 end
             else                                                            % 0.2 Parallel(indPC).Local==0: Run using network on remote machine or also on local machine.
@@ -424,10 +435,21 @@ if parallel_recover ==0
                         command1=[Parallel(indPC).MatlabOctavePath,' -nosplash -nodesktop -minimize ',compThread,' -r "addpath(''',Parallel(indPC).DynarePath,'''), dynareroot = dynare_config(); slaveParallel(',int2str(j),',',int2str(indPC),')" &'];
                     end
                 else    % Hybrid computing Matlab(Master)->Octave(Slaves) and Vice Versa!
+                    if Parallel_info.use_psexec
+                        token1 = ['psexec -accepteula -d -W "',DyMo, '" -a ',my_affinity,' -low  '];
+                    else
+                        itmp = intmin('uint64');
+                        cpus = eval([ '['  my_affinity ']' ])+1;
+                        for icpu=1:length(cpus)
+                            itmp = bitset(itmp,cpus(icpu));
+                        end
+                        hex_affinity = dec2hex(itmp);
+                        token1 = ['start /B /D "',DyMo, '" /affinity ',hex_affinity,' /LOW  '];
+                    end
                     if  regexpi([Parallel(indPC).MatlabOctavePath], 'octave')
-                        command1=['psexec -accepteula -d -W "',DyMo, '" -a ',my_affinity,' -low  ',Parallel(indPC).MatlabOctavePath,' -f --eval "default_save_options(''-v7'');addpath(''',Parallel(indPC).DynarePath,'''), dynareroot = dynare_config(); slaveParallel(',int2str(j),',',int2str(indPC),')"'];
+                        command1=[token1,Parallel(indPC).MatlabOctavePath,' -f --eval "default_save_options(''-v7'');addpath(''',Parallel(indPC).DynarePath,'''), dynareroot = dynare_config(); slaveParallel(',int2str(j),',',int2str(indPC),')"'];
                     else
-                        command1=['psexec -accepteula -d -W "',DyMo, '" -a ',my_affinity,' -low  ',Parallel(indPC).MatlabOctavePath,' -nosplash -nodesktop -minimize ',compThread,' -r "addpath(''',Parallel(indPC).DynarePath,'''), dynareroot = dynare_config(); slaveParallel(',int2str(j),',',int2str(indPC),')"'];
+                        command1=[token1,Parallel(indPC).MatlabOctavePath,' -nosplash -nodesktop -minimize ',compThread,' -r "addpath(''',Parallel(indPC).DynarePath,'''), dynareroot = dynare_config(); slaveParallel(',int2str(j),',',int2str(indPC),')"'];
                     end
                 end
             elseif Parallel(indPC).Local==0                                % 1.2 Run using network on remote machine or also on local machine.
diff --git a/matlab/posterior_sampler_initialization.m b/matlab/posterior_sampler_initialization.m
index b950d3861d30faabb5077417b0e48d59d7f8ee67..89e32fddc60699973eccfd158fbfb3e47f6e48ee 100644
--- a/matlab/posterior_sampler_initialization.m
+++ b/matlab/posterior_sampler_initialization.m
@@ -115,13 +115,80 @@ if ~options_.load_mh_file && ~options_.mh_recover
     if isempty(d)
         prior_draw(bayestopt_,options_.prior_trunc);
     end
+    if options_.mh_initialize_from_previous_mcmc.status
+        PreviousFolder0 = options_.mh_initialize_from_previous_mcmc.directory;
+        RecordFile0 = options_.mh_initialize_from_previous_mcmc.record;
+        PriorFile0 = options_.mh_initialize_from_previous_mcmc.prior;
+        if ~isempty(RecordFile0)
+            %% check for proper filesep char in user defined paths
+            RecordFile0=strrep(RecordFile0,'\',filesep);
+            if isempty(dir(RecordFile0))
+                disp('Estimation::mcmc: wrong value for mh_initialize_from_previous_mcmc_record option')
+                error('Estimation::mcmc: path to record file is not found')
+            else
+                record0=load(RecordFile0);
+            end
+            record0=record0.record;
+            MetropolisFolder0 = fileparts(RecordFile0);
+            PreviousFolder0=fileparts(MetropolisFolder0);
+            [~, ModelName0]=fileparts(PreviousFolder0);
+        else            
+            %% check for proper filesep char in user defined paths
+            PreviousFolder0=strrep(PreviousFolder0,'\',filesep);
+            MetropolisFolder0 = [PreviousFolder0 filesep 'metropolis'];
+            [~, ModelName0]=fileparts(PreviousFolder0);
+            record0=load_last_mh_history_file(MetropolisFolder0, ModelName0);
+        end
+        if ~isnan(record0.MCMCConcludedSuccessfully) && ~record0.MCMCConcludedSuccessfully
+            error('Estimation::mcmc: You are trying to load an MCMC that did not finish successfully. Please use mh_recover.')
+        end
+%         mh_files = dir([ MetropolisFolder0 filesep ModelName0 '_mh*.mat']);
+%         if ~length(mh_files)
+%             error('Estimation::mcmc: I cannot find any MH file to load here!')
+%         end
+        disp('Estimation::mcmc: Initializing from past Metropolis-Hastings simulations...')
+        disp(['Estimation::mcmc: Past MH path ' MetropolisFolder0 ])
+        disp(['Estimation::mcmc: Past model name ' ModelName0 ])
+        fprintf(fidlog,'  Loading initial values from previous MH\n');
+        fprintf(fidlog,'  Past MH path: %s\n', MetropolisFolder0 );
+        fprintf(fidlog,'  Past model name: %s\n', ModelName0);
+        fprintf(fidlog,' \n');
+        past_number_of_blocks = record0.Nblck;
+        if past_number_of_blocks ~= NumberOfBlocks
+            disp('Estimation::mcmc: The specified number of blocks doesn''t match with the previous number of blocks!')
+            disp(['Estimation::mcmc: You declared ' int2str(NumberOfBlocks) ' blocks, but the previous number of blocks was ' int2str(past_number_of_blocks) '.'])
+            disp(['Estimation::mcmc: I will run the Metropolis-Hastings with ' int2str(past_number_of_blocks) ' blocks.' ])
+            NumberOfBlocks = past_number_of_blocks;
+            options_.mh_nblck = NumberOfBlocks;
+        end
+        if ~isempty(PriorFile0)
+            if isempty(dir(PriorFile0))
+                disp('Estimation::mcmc: wrong value for mh_initialize_from_previous_mcmc_prior option')
+                error('Estimation::mcmc: path to prior file is not found')
+            else
+                bayestopt0 = load(PriorFile0);
+            end
+        else
+            bayestopt0 = load([PreviousFolder0 filesep 'prior' filesep 'definition.mat']);
+        end
+        [common_parameters,IA,IB] = intersect(bayestopt_.name,bayestopt0.bayestopt_.name);
+        new_estimated_parameters = ~ismember(bayestopt_.name,bayestopt0.bayestopt_.name);
+        ix2 = zeros(NumberOfBlocks,npar);
+        ilogpo2 = zeros(NumberOfBlocks,1);
+        ilogpo2(:,1) = record0.LastLogPost;
+        ix2(:,IA) = record0.LastParameters(:,IB);
+    else
+        new_estimated_parameters = true(1,npar);
+    end
     % Find initial values for the NumberOfBlocks chains...
-    if NumberOfBlocks > 1% Case 1: multiple chains
+    if NumberOfBlocks > 1 || options_.mh_initialize_from_previous_mcmc.status% Case 1: multiple chains
         set_dynare_seed('default');
         fprintf(fidlog,['  Initial values of the parameters:\n']);
         disp('Estimation::mcmc: Searching for initial values...')
-        ix2 = zeros(NumberOfBlocks,npar);
-        ilogpo2 = zeros(NumberOfBlocks,1);
+        if ~options_.mh_initialize_from_previous_mcmc.status
+            ix2 = zeros(NumberOfBlocks,npar);
+            ilogpo2 = zeros(NumberOfBlocks,1);
+        end
         for j=1:NumberOfBlocks
             validate    = 0;
             init_iter   = 0;
@@ -133,7 +200,7 @@ if ~options_.load_mh_file && ~options_.mh_recover
                     candidate = rand_multivariate_normal( transpose(xparam1), d * options_.mh_init_scale, npar);
                 end
                 if all(candidate(:) >= mh_bounds.lb) && all(candidate(:) <= mh_bounds.ub)
-                    ix2(j,:) = candidate;
+                    ix2(j,new_estimated_parameters) = candidate(new_estimated_parameters);
                     ilogpo2(j) = - feval(TargetFun,ix2(j,:)',dataset_,dataset_info,options_,M_,estim_params_,bayestopt_,mh_bounds,oo_);
                     if ~isfinite(ilogpo2(j)) % if returned log-density is
                                              % Inf or Nan (penalized value)
@@ -212,11 +279,16 @@ if ~options_.load_mh_file && ~options_.mh_recover
     record.MAX_nruns=MAX_nruns;
     record.AcceptanceRatio = zeros(1,NumberOfBlocks);
     record.FunctionEvalPerIteration = NaN(1,NumberOfBlocks);
-    for j=1:NumberOfBlocks
-        % we set a different seed for the random generator for each block then we record the corresponding random generator state (vector)
-        set_dynare_seed(options_.DynareRandomStreams.seed+j);
-        % record.Seeds keeps a vector of the random generator state and not the scalar seed despite its name
-        [record.InitialSeeds(j).Unifor,record.InitialSeeds(j).Normal] = get_dynare_random_generator_state();
+    if options_.mh_initialize_from_previous_mcmc.status
+        record.InitialSeeds = record0.LastSeeds;
+    else
+        
+        for j=1:NumberOfBlocks
+            % we set a different seed for the random generator for each block then we record the corresponding random generator state (vector)
+            set_dynare_seed(options_.DynareRandomStreams.seed+j);
+            % record.Seeds keeps a vector of the random generator state and not the scalar seed despite its name
+            [record.InitialSeeds(j).Unifor,record.InitialSeeds(j).Normal] = get_dynare_random_generator_state();
+        end
     end
     record.InitialParameters = ix2;
     record.InitialLogPost = ilogpo2;
diff --git a/matlab/realtime_shock_decomposition.m b/matlab/realtime_shock_decomposition.m
index 6686d5c8e0bb7f39236fb4596b9b951f36ef1ec3..ab29c861fe7c153ace50584103c9a5fce14e8b67 100644
--- a/matlab/realtime_shock_decomposition.m
+++ b/matlab/realtime_shock_decomposition.m
@@ -116,56 +116,55 @@ if forecast_ && any(forecast_params)
     [~,~,~,~,~,~,oo1] = dynare_resolve(M1,options_,oo_);
 end
 
-if fast_realtime
-    skipline()
-    skipline()
-    running_text = 'Fast realtime shock decomposition ';
-    newString=sprintf(running_text);
-    fprintf(['%s'],newString);
-    options_.nobs=fast_realtime;
-    [oo0,M_,~,~,Smoothed_Variables_deviation_from_mean0] = evaluate_smoother(parameter_set,varlist,M_,oo_,options_,bayestopt_,estim_params_);
-    gend0 = size(oo0.SmoothedShocks.(M_.exo_names{1}),1);
-    prctdone=0.5;
-    if isoctave
-        printf([running_text,' %3.f%% done\r'], prctdone*100);
-    else
-        s0=repmat('\b',1,length(newString)+1);
-        newString=sprintf([running_text,' %3.1f%% done'], prctdone*100);
-        fprintf([s0,'%s'],newString);
-    end
-    options_.nobs=nobs;
-    [oo2,M_,~,~,Smoothed_Variables_deviation_from_mean2] = evaluate_smoother(parameter_set,varlist,M_,oo_,options_,bayestopt_,estim_params_);
-    prctdone=1;
-    if isoctave
-        printf([running_text,' %3.f%% done\r'], prctdone*100);
-    else
-        s0=repmat('\b',1,length(newString)+1);
-        newString=sprintf([running_text,' %3.1f%% done'], prctdone*100);
-        fprintf([s0,'%s'],newString);
-    end
-end
+gend0=0;
 
 skipline()
 skipline()
-running_text = 'Realtime shock decomposition ';
+if isequal(fast_realtime,0)
+    running_text = 'Realtime shock decomposition ';
+else
+    running_text = 'Fast realtime shock decomposition ';
+end
 newString=sprintf(running_text);
 fprintf(['%s'],newString);
 
 for j=presample+1:nobs
     %    evalin('base',['options_.nobs=' int2str(j) ';'])
     options_.nobs=j;
-    if ~fast_realtime
+    if isequal(fast_realtime,0)
         [oo,M_,~,~,Smoothed_Variables_deviation_from_mean] = evaluate_smoother(parameter_set,varlist,M_,oo_,options_,bayestopt_,estim_params_);
         gend = size(oo.SmoothedShocks.(M_.exo_names{1}),1);
     else
-        gend = gend0+j-fast_realtime;
-        if j>fast_realtime
-            oo=oo2;
-            Smoothed_Variables_deviation_from_mean = Smoothed_Variables_deviation_from_mean2(:,1:gend);
+        if j<min(fast_realtime) && gend0<j
+            options_.nobs=min(fast_realtime);
+            [oo0,M_,~,~,Smoothed_Variables_deviation_from_mean0] = evaluate_smoother(parameter_set,varlist,M_,oo_,options_,bayestopt_,estim_params_);
+            gend0 = size(oo0.SmoothedShocks.(M_.exo_names{1}),1);
+            options_.nobs=j;
+        end
+        
+        if ismember(j,fast_realtime) && gend0<j
+            [oo,M_,~,~,Smoothed_Variables_deviation_from_mean] = evaluate_smoother(parameter_set,varlist,M_,oo_,options_,bayestopt_,estim_params_);
+            gend = size(oo.SmoothedShocks.(M_.exo_names{1}),1);
+            gend0 = gend;
+            oo0=oo;
+            Smoothed_Variables_deviation_from_mean0=Smoothed_Variables_deviation_from_mean;
         else
+            if j>gend0
+                if j>max(fast_realtime)
+                    options_.nobs = nobs;
+                else
+                    options_.nobs=min(fast_realtime(fast_realtime>j));
+                end
+                [oo0,M_,~,~,Smoothed_Variables_deviation_from_mean0] = evaluate_smoother(parameter_set,varlist,M_,oo_,options_,bayestopt_,estim_params_);
+                gend0 = size(oo0.SmoothedShocks.(M_.exo_names{1}),1);
+                options_.nobs=j;
+            end
+            
+            gend = j;
             oo=oo0;
             Smoothed_Variables_deviation_from_mean = Smoothed_Variables_deviation_from_mean0(:,1:gend);
         end
+        
     end
     % reduced form
     dr = oo.dr;
diff --git a/matlab/simult_.m b/matlab/simult_.m
index 508af26e435fe890db154dd1a1f9733c29b46911..1ab1729cfebcc9780e6c2dde88ec6f36a431a4a6 100644
--- a/matlab/simult_.m
+++ b/matlab/simult_.m
@@ -42,7 +42,8 @@ y_ = zeros(size(y0,1),iter+M_.maximum_lag);
 y_(:,1) = y0;
 
 if options_.loglinear && ~options_.logged_steady_state
-    dr.ys=log(dr.ys);
+    k = get_all_variables_but_lagged_leaded_exogenous(M_);
+    dr.ys(k)=log(dr.ys(k));
 end
 
 if ~options_.k_order_solver || (options_.k_order_solver && options_.pruning) %if k_order_pert is not used or if we do not use Dynare++ with k_order_pert
diff --git a/matlab/store_smoother_results.m b/matlab/store_smoother_results.m
index cc529fab24c183b9a3d991c8fb92b62459d28bf6..b53283412888242bdead1eb75d649d068bd2c433 100644
--- a/matlab/store_smoother_results.m
+++ b/matlab/store_smoother_results.m
@@ -34,8 +34,9 @@ function [oo_, yf]=store_smoother_results(M_,oo_,options_,bayestopt_,dataset_,da
 %                   oo_.Smoother.Variance: one-step ahead forecast error variance (declaration order)
 %                   oo_.Smoother.Constant: structure storing the constant term of the smoother
 %                   oo_.Smoother.Trend: structure storing the trend term of the smoother
-%                   oo_.FilteredVariablesKStepAhead: k-step ahead forecast error variance matrices (decision-rule order)
-%                   oo_.FilteredVariablesShockDecomposition: shock decomposition of k-step ahead filtered variables (decision-rule order)
+%                   oo_.FilteredVariablesKStepAhead: k-step ahead filtered variables (declaration order)
+%                   oo_.FilteredVariablesKStepAheadVariances: k-step ahead forecast error variance matrices (declaration order)
+%                   oo_.FilteredVariablesShockDecomposition: shock decomposition of k-step ahead filtered variables (declaration order)
 %                   oo_.FilteredVariables: structure storing the filtered variables
 %                   oo_.UpdatedVariables: structure storing the updated variables
 %                   oo_.SmoothedShocks: structure storing the smoothed shocks
@@ -78,6 +79,12 @@ end
 
 if options_.loglinear
     oo_.Smoother.loglinear = true;
+    if ~isempty(M_.aux_vars) % deal with lead/lag of exogenous variables
+        exo_lead_lag_index = M_.orig_endo_nbr + find(([M_.aux_vars.type] == 2) | ([M_.aux_vars.type]  == 3));
+    else
+        exo_lead_lag_index =[];
+    end
+
 else
     oo_.Smoother.loglinear = false;
 end
@@ -136,6 +143,9 @@ if ~isempty(options_.nk) && options_.nk ~= 0
     i_endo_declaration_order = oo_.dr.order_var(i_endo_in_dr_matrices); %get indices of smoothed variables in name vector
     if  options_.loglinear %logged steady state must be used
         constant_all_variables=repmat(log(ys(i_endo_declaration_order))',[length(options_.filter_step_ahead),1,gend+max(options_.filter_step_ahead)]);
+        if ~isempty(exo_lead_lag_index) && any(ismember(i_endo_declaration_order,exo_lead_lag_index)) % deal with lead/lag of exogenous variables
+            constant_all_variables(:,ismember(i_endo_declaration_order,exo_lead_lag_index),:)=repmat(ys(i_endo_declaration_order(ismember(i_endo_declaration_order,exo_lead_lag_index)))',[length(options_.filter_step_ahead),1,gend+max(options_.filter_step_ahead)]);
+        end
     elseif ~options_.loglinear %unlogged steady state must be used
         constant_all_variables=repmat((ys(i_endo_declaration_order))',[length(options_.filter_step_ahead),1,gend+max(options_.filter_step_ahead)]);
     end
@@ -156,7 +166,11 @@ for i_endo_in_bayestopt_smoother_varlist=bayestopt_.smoother_saved_var_list'
     i_endo_declaration_order = oo_.dr.order_var(i_endo_in_dr); %get indices of smoothed variables in name vector
     %% Compute constant
     if  options_.loglinear == 1 %logged steady state must be used
-        constant_current_variable=repmat(log(ys(i_endo_declaration_order)),gend,1);
+        if ~isempty(exo_lead_lag_index) && ismember(i_endo_declaration_order,exo_lead_lag_index)% deal with lead/lag of exogenous variables
+            constant_current_variable=repmat(ys(i_endo_declaration_order),gend,1);
+        else
+            constant_current_variable=repmat(log(ys(i_endo_declaration_order)),gend,1);
+        end
     elseif options_.loglinear == 0 %unlogged steady state must be used
         constant_current_variable=repmat((ys(i_endo_declaration_order)),gend,1);
     end
diff --git a/mex/build/matlab/configure.ac b/mex/build/matlab/configure.ac
index c9460676a4172feab47f6c316d34f1c8f1b22cb1..63ddf61a248bc148386998481fbc171d3b553e71 100644
--- a/mex/build/matlab/configure.ac
+++ b/mex/build/matlab/configure.ac
@@ -45,10 +45,6 @@ case ${host_os} in
     ;;
 esac
 
-CFLAGS="$MATLAB_CFLAGS -Wall -Wno-parentheses"
-FCFLAGS="$MATLAB_FCFLAGS -Wall"
-CXXFLAGS="$MATLAB_CXXFLAGS -Wall -Wno-parentheses -Wold-style-cast"
-
 AC_PROG_FC
 AC_PROG_CC
 AC_PROG_CC_C99
@@ -89,7 +85,13 @@ fi
 
 # Check for libslicot, needed by kalman_steady_state
 if test "$enable_mex_kalman_steady_state" = yes; then
+  # FCFLAGS must be temporarily modified, because otherwise -fno-underscoring is not
+  # taken into account by AC_FC_FUNC in the AX_SLICOT macro.
+  # For some obscure reason, it is necessary to do it at this level and not within the macro.
+  ac_save_FCFLAGS=$FCFLAGS
+  FCFLAGS="$FCFLAGS $MATLAB_FCFLAGS"
   AX_SLICOT([matlab])
+  FCFLAGS=$ac_save_FCFLAGS
   test "$has_slicot" != yes && AC_MSG_ERROR([slicot cannot be found. If you want to skip the compilation of the kalman_steady_state MEX, pass the --disable-mex-kalman-steady-state flag.])
 fi
 
diff --git a/mex/build/matlab/mex.am b/mex/build/matlab/mex.am
index 5ba01b72306f8c083666f8f115302fe7d4c23249..0bd51ad0bf1f129903328bd25ff76a87fb7b732e 100644
--- a/mex/build/matlab/mex.am
+++ b/mex/build/matlab/mex.am
@@ -7,6 +7,9 @@ DEFS += $(MATLAB_DEFS)
 DEFS += -DMATLAB_MEX_FILE
 DEFS += -DMEXEXT=\"$(MEXEXT)\"
 
+AM_CFLAGS = $(MATLAB_CFLAGS) -Wall -Wno-parentheses
+AM_FCFLAGS = $(MATLAB_FCFLAGS) -Wall
+AM_CXXFLAGS = $(MATLAB_CXXFLAGS) -Wall -Wno-parentheses -Wold-style-cast
 AM_LDFLAGS = $(MATLAB_LDFLAGS)
 LIBS += $(MATLAB_LIBS)
 
diff --git a/mex/build/octave/configure.ac b/mex/build/octave/configure.ac
index b9ffc14ce25e084e92fe9515dd74b2c311b96d4f..1ce99dbb646ccb22ff9d82ca341b134dd5854312 100644
--- a/mex/build/octave/configure.ac
+++ b/mex/build/octave/configure.ac
@@ -26,14 +26,14 @@ AX_OCTAVE
 
 test "$ax_enable_octave" != yes && AC_MSG_ERROR([Octave cannot be found])
 
+# Let mkoctfile set the default compilers and flags (except for FC and FCFLAGS which are not supported)
+# NB: mkoctfile honors overrides via environment variables
 CC=$($MKOCTFILE -p CC)
 CXX=$($MKOCTFILE -p CXX)
 AR=$($MKOCTFILE -p AR)
 RANLIB=$($MKOCTFILE -p RANLIB)
-CFLAGS="$($MKOCTFILE -p CFLAGS) -Wall -Wno-parentheses"
-FCFLAGS="$($MKOCTFILE -p FFLAGS) -Wall -std=gnu" # Override -std=legacy that is in FFLAGS in Octave 5
-FFLAGS="$($MKOCTFILE -p FFLAGS) -Wall"
-CXXFLAGS="$($MKOCTFILE -p CXXFLAGS) -Wall -Wno-parentheses -Wold-style-cast"
+CFLAGS=$($MKOCTFILE -p CFLAGS)
+CXXFLAGS=$($MKOCTFILE -p CXXFLAGS)
 LDFLAGS="$($MKOCTFILE -p LFLAGS) $($MKOCTFILE -p LDFLAGS)"
 
 AC_CANONICAL_HOST
diff --git a/mex/build/octave/mex.am b/mex/build/octave/mex.am
index 77eccb3e7df62bc644d46f4e619ef933dc10b922..4766d5d25f7344af4e008a729e6c91b5188be052 100644
--- a/mex/build/octave/mex.am
+++ b/mex/build/octave/mex.am
@@ -5,10 +5,9 @@ AM_CPPFLAGS += -I$(top_srcdir)/../../sources
 DEFS += -DOCTAVE_MEX_FILE
 DEFS += -DMEXEXT=\".mex\"
 
-AM_CFLAGS = $(shell $(MKOCTFILE) -p CPICFLAG)
-AM_FFLAGS = $(shell $(MKOCTFILE) -p FPICFLAG)
-AM_FCFLAGS = $(shell $(MKOCTFILE) -p FPICFLAG)
-AM_CXXFLAGS = $(shell $(MKOCTFILE) -p CXXPICFLAG)
+AM_CFLAGS = $(shell $(MKOCTFILE) -p CPICFLAG) -Wall -Wno-parentheses
+AM_FCFLAGS = $(shell $(MKOCTFILE) -p FPICFLAG) -Wall
+AM_CXXFLAGS = $(shell $(MKOCTFILE) -p CXXPICFLAG) -Wall -Wno-parentheses -Wold-style-cast
 AM_LDFLAGS = $(shell $(MKOCTFILE) -p DL_LDFLAGS) -L"$(shell $(MKOCTFILE) -p OCTLIBDIR)"
 
 LIBS += $(shell $(MKOCTFILE) -p OCTAVE_LIBS)
@@ -16,7 +15,6 @@ LIBS += $(shell $(MKOCTFILE) -p BLAS_LIBS)
 LIBS += $(shell $(MKOCTFILE) -p LAPACK_LIBS)
 LIBS += $(shell $(MKOCTFILE) -p FFTW_LIBS)
 LIBS += $(shell $(MKOCTFILE) -p LIBS)
-LIBS += $(shell $(MKOCTFILE) -p FLIBS)
 
 mexdir = $(libdir)/dynare/mex/octave
 
diff --git a/preprocessor b/preprocessor
index 910649dd28ea3f78d6cf4cb245a2747207828bf5..a2c96a9d7927a47ad6017c64df55f32595ccb1e5 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit 910649dd28ea3f78d6cf4cb245a2747207828bf5
+Subproject commit a2c96a9d7927a47ad6017c64df55f32595ccb1e5
diff --git a/scripts/dynare.el b/scripts/dynare.el
index b6d102b49daec856fd3a3d9227b77f346d5064ce..59153dd4ccb131d1a0755452ac583ec5bcc1c4c8 100644
--- a/scripts/dynare.el
+++ b/scripts/dynare.el
@@ -1,7 +1,7 @@
 ;;; dynare.el --- major mode for editing Dynare mod files
 
 ;; Copyright © 2010 Yannick Kalantzis
-;; Copyright © 2019-2020 Dynare Team
+;; Copyright © 2019-2021 Dynare Team
 ;;
 ;; This program is free software; you can redistribute it and/or modify
 ;; it under the terms of the GNU General Public License as published by
@@ -94,7 +94,7 @@
     '("model" "steady_state_model" "initval" "endval" "histval" "shocks"
       "shock_groups" "init2shocks" "mshocks" "estimated_params" "epilogue" "priors"
       "estimated_params_init" "estimated_params_bounds" "osr_params_bounds"
-      "observation_trends" "optim_weights" "homotopy_setup"
+      "observation_trends" "deterministic_trends" "optim_weights" "homotopy_setup"
       "conditional_forecast_paths" "svar_identification" "moment_calibration"
       "irf_calibration" "ramsey_constraints" "generate_irfs"
       "matched_moments" "verbatim")
diff --git a/tests/.gitignore b/tests/.gitignore
index b1d306b68fd5eedc797c3ba5929eafd51aeb3719..355c011cf96a51aae73316261c0dd3868ba50815 100644
--- a/tests/.gitignore
+++ b/tests/.gitignore
@@ -54,6 +54,9 @@ wsOct
 !/estimation/method_of_moments/RBC/RBC_Andreasen_Data_2.mat
 !/estimation/method_of_moments/AFVRR/AFVRR_data.mat
 !/estimation/method_of_moments/AFVRR/AFVRR_steady_helper.m
+!estimation/no_init_estimation_check_first_obs/fsdat_mat.m
+!estimation/no_init_estimation_check_first_obs/fsdat_mat_XFAIL.m
+!estimation/no_init_estimation_check_first_obs/fsdat_simul.m
 !/expectations/expectation_ss_old_steadystate.m
 !/external_function/extFunDeriv.m
 !/external_function/extFunNoDerivs.m
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 86d3e1a42b4ffbc97fc24594cda0d5bbcdbd08da..7befb095f5e100c32abb355f864132397632d8da 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -44,6 +44,7 @@ MODFILES = \
 	estimation/fs2000_calibrated_covariance.mod \
 	estimation/fs2000_model_comparison.mod \
 	estimation/fs2000_fast.mod \
+	estimation/fs2000_init_from_previous.mod \
 	estimation/ls2003_endog_prior_restrict_estimation.mod \
 	estimation/independent_mh/fs2000_independent_mh.mod \
 	estimation/MH_recover/fs2000_recover.mod \
@@ -60,6 +61,7 @@ MODFILES = \
 	estimation/method_of_moments/AFVRR/AFVRR_M0.mod \
 	estimation/method_of_moments/AFVRR/AFVRR_MFB.mod \
 	estimation/method_of_moments/AFVRR/AFVRR_MFB_RRA.mod \
+	estimation/no_init_estimation_check_first_obs/fs2000_init_check.mod \
 	moments/example1_var_decomp.mod \
 	moments/example1_bp_test.mod \
 	moments/test_AR1_spectral_density.mod \
@@ -158,6 +160,12 @@ MODFILES = \
 	steady_state_operator/bytecode_test.mod \
 	block_bytecode/ireland.mod \
 	block_bytecode/ramst_normcdf_and_friends.mod \
+	block_bytecode/lola_solve_one_boundary.mod \
+	block_bytecode/lola_solve_one_boundary_mfs1.mod \
+	block_bytecode/lola_solve_one_boundary_mfs2.mod \
+	block_bytecode/lola_solve_one_boundary_mfs3.mod \
+	block_bytecode/lola_stochastic.mod \
+	block_bytecode/lola_stochastic_block.mod \
 	k_order_perturbation/fs2000k2a.mod \
 	k_order_perturbation/fs2000k2_use_dll.mod \
 	k_order_perturbation/fs2000k_1_use_dll.mod \
@@ -338,10 +346,6 @@ MODFILES = \
 	deterministic_simulations/multiple_lead_lags/sim_endo_lead_lag.mod \
 	deterministic_simulations/multiple_lead_lags/sim_lead_lag_aux_vars.mod \
 	deterministic_simulations/multiple_lead_lags/sim_lead_lag.mod \
-	deterministic_simulations/lola_solve_one_boundary.mod \
-	deterministic_simulations/lola_solve_one_boundary_mfs1.mod \
-	deterministic_simulations/lola_solve_one_boundary_mfs2.mod \
-	deterministic_simulations/lola_solve_one_boundary_mfs3.mod \
 	deterministic_simulations/ramst_block_mfs1.mod \
 	deterministic_simulations/linear_approximation/sw.mod \
 	deterministic_simulations/multiple_lead_lags/AR2.mod \
@@ -571,7 +575,9 @@ XFAIL_MODFILES = ramst_xfail.mod \
 	kalman_initial_state/fs2000_kalman_initial_xfail.mod \
 	example1_extra_exo_xfail.mod \
 	estimation/tune_mh_jscale/fs2000_1_xfail.mod \
-	estimation/tune_mh_jscale/fs2000_2_xfail.mod
+	estimation/tune_mh_jscale/fs2000_2_xfail.mod \
+	estimation/no_init_estimation_check_first_obs/fs2000_init_check_XFAIL.mod \
+	estimation/no_init_estimation_check_first_obs/fs2000_init_check_XFAIL2.mod
 
 MFILES = histval_initval_file/ramst_initval_file_data.m
 
@@ -642,6 +648,8 @@ AIM/ls2003_2L2L_AIM.o.trs: AIM/ls2003_2L2L.o.trs
 
 estimation/fs2000_model_comparison.m.trs: estimation/fs2000.m.trs estimation/fs2000_initialize_from_calib.m.trs estimation/fs2000_calibrated_covariance.m.trs
 estimation/fs2000_model_comparison.o.trs: estimation/fs2000.o.trs estimation/fs2000_initialize_from_calib.o.trs estimation/fs2000_calibrated_covariance.o.trs
+estimation/fs2000_init_from_previous.m.trs: estimation/fs2000.m.trs
+estimation/fs2000_init_from_previous.o.trs: estimation/fs2000.o.trs
 
 optimizers/fs2000_102.m.trs: estimation/fs2000.m.trs
 optimizers/fs2000_102.o.trs: estimation/fs2000.o.trs
@@ -696,13 +704,15 @@ deterministic_simulations/rbc_det_stack_solve_algo_7_exo_lead.m.trs: determinist
 deterministic_simulations/rbc_det_stack_solve_algo_7_exo_lead.o.trs: deterministic_simulations/rbc_det.o.trs
 deterministic_simulations/rbc_det_stack_solve_algo_7_exo_lag.m.trs: deterministic_simulations/rbc_det.m.trs
 deterministic_simulations/rbc_det_stack_solve_algo_7_exo_lag.o.trs: deterministic_simulations/rbc_det.o.trs
-deterministic_simulations/lola_solve_one_boundary_mfs1.m.trs: deterministic_simulations/lola_solve_one_boundary.m.trs
-deterministic_simulations/lola_solve_one_boundary_mfs1.o.trs: deterministic_simulations/lola_solve_one_boundary.o.trs
-deterministic_simulations/lola_solve_one_boundary_mfs2.m.trs: deterministic_simulations/lola_solve_one_boundary.m.trs
-deterministic_simulations/lola_solve_one_boundary_mfs2.o.trs: deterministic_simulations/lola_solve_one_boundary.o.trs
-deterministic_simulations/lola_solve_one_boundary_mfs3.m.trs: deterministic_simulations/lola_solve_one_boundary.m.trs
-deterministic_simulations/lola_solve_one_boundary_mfs3.o.trs: deterministic_simulations/lola_solve_one_boundary.o.trs
 
+block_bytecode/lola_solve_one_boundary_mfs1.m.trs: block_bytecode/lola_solve_one_boundary.m.trs
+block_bytecode/lola_solve_one_boundary_mfs1.o.trs: block_bytecode/lola_solve_one_boundary.o.trs
+block_bytecode/lola_solve_one_boundary_mfs2.m.trs: block_bytecode/lola_solve_one_boundary.m.trs
+block_bytecode/lola_solve_one_boundary_mfs2.o.trs: block_bytecode/lola_solve_one_boundary.o.trs
+block_bytecode/lola_solve_one_boundary_mfs3.m.trs: block_bytecode/lola_solve_one_boundary.m.trs
+block_bytecode/lola_solve_one_boundary_mfs3.o.trs: block_bytecode/lola_solve_one_boundary.o.trs
+block_bytecode/lola_stochastic_block.m.trs: block_bytecode/lola_stochastic.m.trs
+block_bytecode/lola_stochastic_block.o.trs: block_bytecode/lola_stochastic.o.trs
 
 histval_initval_file/ramst_initval_file.m.trs: histval_initval_file/ramst_initval_file_data.m.tls histval_initval_file/ramst_data_generate.m.trs
 histval_initval_file/ramst_initval_file.o.trs: histval_initval_file/ramst_initval_file_data.o.tls histval_initval_file/ramst_data_generate.o.trs
@@ -1109,6 +1119,8 @@ EXTRA_DIST = \
 	AIM/data_ca1.m \
 	AIM/fsdat.m \
 	block_bytecode/run_ls2003.m \
+	block_bytecode/lola_data.mat \
+	block_bytecode/lola_common.inc \
 	bvar_a_la_sims/bvar_sample.m \
 	dates/fsdat_simul.m \
 	dates/data_uav.xlsx \
@@ -1199,7 +1211,6 @@ EXTRA_DIST = \
 	optimizers/fs2000.common.inc \
 	estimation/MH_recover/fs2000.common.inc \
 	prior_posterior_function/posterior_function_demo.m \
-	deterministic_simulations/lola_data.mat \
 	k_order_perturbation/fs2000k++.mod \
 	lmmcp/sw-common-header.inc \
 	lmmcp/sw-common-footer.inc \
@@ -1256,7 +1267,7 @@ check-matlab-ols: $(M_OLS_TRS_FILES)
 
 %.m.trs %.m.log: %.mod
 	@echo "`tput bold``tput setaf 8`MATLAB: $(CURDIR)/$*... `tput sgr0`"
-# The while loop is a workaround for this glibc bug: https://sourceware.org/bugzilla/show_bug.cgi?id=19329
+# The while loop is a workaround for this glibc bug: https://sourceware.org/bugzilla/show_bug.cgi?id=19329 (should be fixed in glibc 2.34, likely included in Debian “bookworm” 12)
 	@while ! [ -f $*.m.trs ] || head -1 $*.m.log | egrep -q '^Inconsistency detected by ld.so'; \
 	do \
 	DYNARE_VERSION="$(PACKAGE_VERSION)" TOP_TEST_DIR="$(CURDIR)" FILESTEM="$*" \
@@ -1275,7 +1286,7 @@ check-matlab-ols: $(M_OLS_TRS_FILES)
 
 %.m.trs %.m.log : %.m
 	@echo "`tput bold``tput setaf 8`MATLAB: $(CURDIR)/$*... `tput sgr0`"
-# The while loop is a workaround for this glibc bug: https://sourceware.org/bugzilla/show_bug.cgi?id=19329
+# The while loop is a workaround for this glibc bug: https://sourceware.org/bugzilla/show_bug.cgi?id=19329 (should be fixed in glibc 2.34, likely included in Debian “bookworm” 12)
 	@while ! [ -f $*.m.trs ] || head -1 $*.m.log | egrep -q '^Inconsistency detected by ld.so'; \
 	do \
 	DYNARE_VERSION="$(PACKAGE_VERSION)" TOP_TEST_DIR="$(CURDIR)" \
diff --git a/tests/deterministic_simulations/lola_solve_one_boundary.mod b/tests/block_bytecode/lola_common.inc
similarity index 88%
rename from tests/deterministic_simulations/lola_solve_one_boundary.mod
rename to tests/block_bytecode/lola_common.inc
index 6930b4d7240974198371bffd9a3a762c02e2c171..ae582cea5067d17918d37d4570a88ea15bdaf7a3 100644
--- a/tests/deterministic_simulations/lola_solve_one_boundary.mod
+++ b/tests/block_bytecode/lola_common.inc
@@ -1,5 +1,13 @@
-//  Based on Luca Marchiori/Olivier Pierrard(2012) LOLA 2.0: Luxembourg OverLapping generation model for policy Analysis
-// Involves a call to solve_one_boundary.m that is tested here
+//  Based on Luca Marchiori/Olivier Pierrard(2012) LOLA 2.0: Luxembourg
+//  OverLapping generation model for policy Analysis
+
+/* Snippet common to several test files.
+
+   Macro-variables that can be defined to tune the computations:
+   - block
+   - mfs
+   - deterministic
+*/
 
 load lola_data.mat
 
@@ -144,9 +152,12 @@ set_param_value('bsY_iss',bsY_iss);
 NBRYb=NBR_iss/(phii_iss*gdp_iss);
 
 
-% =======================================================
-model(block);              
-% ======================================================
+
+@#if block
+model(block, mfs=@{mfs});
+@#else
+model;
+@#endif
 
 %  Labor Market Variables in the home country
 %  ------------------------------------------ 
@@ -705,18 +716,16 @@ tauc2=tauc;
 
 end;
 
-%========================================================
-% compute initial steady state and check eigenvalues 
-%========================================================
+%==============================
+% compute initial steady state
+%==============================
 
 resid;
 steady(solve_algo=3);
-check;
 
-%========================================================
-endval;
-%========================================================
+@#if deterministic
 
+endval;
 @#for i in wg
 n@{i}=n@{i}_fss;
 n@{i}_f=n@{i}_f_fss;
@@ -918,54 +927,50 @@ NBR2=NBR;
 tauf2=tauf;
 tauw2=tauw;
 tauc2=tauc;
-
 end;
 
-%========================================================
-% compute final steady state and check eigenvalues 
-%========================================================
+%============================
+% compute final steady state
+%============================
 
 resid;
 steady(solve_algo=3);
-check;
 
 
-% ===================================================
 shocks;     
-% ===================================================
-
-var P00;																																																			
-   periods 1:99;																																																			
-   values (se_P00);																																																			
+  var P00;
+  periods 1:99;
+  values (se_P00);
  
 @#for i in nbg
-var beta@{i};																																																			
-   periods  1:99;																																																			
-   values (se_beta@{i});
-var beta@{i}_f;																																																			
-   periods  1:99;																																																			
-   values (se_beta@{i}_f);
-var PD@{i};																																																			
-   periods  1:99;																																																			
-   values (se_PD@{i});
+  var beta@{i};
+  periods  1:99;
+  values (se_beta@{i});
+
+  var beta@{i}_f;
+  periods  1:99;
+  values (se_beta@{i}_f);
+
+  var PD@{i};
+  periods  1:99;
+  values (se_PD@{i});
 @#endfor  
-																																																			
-var P00_foP00;																																																			
-   periods 1:99;																																																			
-   values (se_P00_foP00);	 																																																
 
-var eps_g;																																																			
-   periods 1:99;																																																			
-   values (se_eps_g);	
+  var P00_foP00;
+  periods 1:99;
+  values (se_P00_foP00);
 
-var eps_PensCorr_F;																																																			
-   periods 1:99;																																																			
-   values (se_eps_PensCorr_F);	
+  var eps_g;
+  periods 1:99;
+  values (se_eps_g);
 
-var eps_PensCorr_L;																																																			
-   periods 1:99;																																																			
-   values (se_eps_PensCorr_L);	
+  var eps_PensCorr_F;
+  periods 1:99;
+  values (se_eps_PensCorr_F);	
 
+  var eps_PensCorr_L;
+  periods 1:99;
+  values (se_eps_PensCorr_L);	
 end;
 
 % *******************************************
@@ -978,3 +983,7 @@ perfect_foresight_solver(maxit=100);
 if ~oo_.deterministic_simulation.status
    error('Perfect foresight simulation failed')
 end
+
+@#else // stochastic case, used by files under tests/block_bytecode/lola_*
+stoch_simul(order=1);
+@#endif
diff --git a/tests/deterministic_simulations/lola_data.mat b/tests/block_bytecode/lola_data.mat
similarity index 100%
rename from tests/deterministic_simulations/lola_data.mat
rename to tests/block_bytecode/lola_data.mat
diff --git a/tests/block_bytecode/lola_solve_one_boundary.mod b/tests/block_bytecode/lola_solve_one_boundary.mod
new file mode 100644
index 0000000000000000000000000000000000000000..f72b8afca9ac3aac152f1a08523d27fce687cdd8
--- /dev/null
+++ b/tests/block_bytecode/lola_solve_one_boundary.mod
@@ -0,0 +1,6 @@
+// Involves a call to solve_one_boundary.m that is tested here
+
+@#define deterministic = true
+@#define block = true
+@#define mfs = 0
+@#include "lola_common.inc"
diff --git a/tests/block_bytecode/lola_solve_one_boundary_mfs1.mod b/tests/block_bytecode/lola_solve_one_boundary_mfs1.mod
new file mode 100644
index 0000000000000000000000000000000000000000..553171257d26babb89eed0970c6387159fa33b1d
--- /dev/null
+++ b/tests/block_bytecode/lola_solve_one_boundary_mfs1.mod
@@ -0,0 +1,12 @@
+// Tests option mfs=1 with block
+
+@#define deterministic = true
+@#define block = true
+@#define mfs = 1
+@#include "lola_common.inc"
+
+mfs0=load('lola_solve_one_boundary_results');
+
+if max(max(oo_.endo_simul-mfs0.oo_.endo_simul)) > options_.dynatol.x
+   error('Inconsistency with mfs=0')
+end
diff --git a/tests/block_bytecode/lola_solve_one_boundary_mfs2.mod b/tests/block_bytecode/lola_solve_one_boundary_mfs2.mod
new file mode 100644
index 0000000000000000000000000000000000000000..f582ecc8440f3064de4795c1addcc861786f0694
--- /dev/null
+++ b/tests/block_bytecode/lola_solve_one_boundary_mfs2.mod
@@ -0,0 +1,12 @@
+// Tests option mfs=2 with block
+
+@#define deterministic = true
+@#define block = true
+@#define mfs = 2
+@#include "lola_common.inc"
+
+mfs0=load('lola_solve_one_boundary_results');
+
+if max(max(oo_.endo_simul-mfs0.oo_.endo_simul)) > options_.dynatol.x
+   error('Inconsistency with mfs=0')
+end
diff --git a/tests/block_bytecode/lola_solve_one_boundary_mfs3.mod b/tests/block_bytecode/lola_solve_one_boundary_mfs3.mod
new file mode 100644
index 0000000000000000000000000000000000000000..3eb2dc9f0178f6cdeaa03407891c8d3126007e9b
--- /dev/null
+++ b/tests/block_bytecode/lola_solve_one_boundary_mfs3.mod
@@ -0,0 +1,12 @@
+// Tests option mfs=3 with block
+
+@#define deterministic = true
+@#define block = true
+@#define mfs = 3
+@#include "lola_common.inc"
+
+mfs0=load('lola_solve_one_boundary_results');
+
+if max(max(oo_.endo_simul-mfs0.oo_.endo_simul)) > options_.dynatol.x
+   error('Inconsistency with mfs=0')
+end
diff --git a/tests/block_bytecode/lola_stochastic.mod b/tests/block_bytecode/lola_stochastic.mod
new file mode 100644
index 0000000000000000000000000000000000000000..e6c0071fc93c53c743392447132b52df66ac1ba8
--- /dev/null
+++ b/tests/block_bytecode/lola_stochastic.mod
@@ -0,0 +1,5 @@
+// Stochastic version of LOLA model, as benchmark for lola_stochastic_block.mod
+
+@#define deterministic = false
+@#define block = false
+@#include "lola_common.inc"
diff --git a/tests/block_bytecode/lola_stochastic_block.mod b/tests/block_bytecode/lola_stochastic_block.mod
new file mode 100644
index 0000000000000000000000000000000000000000..5f6fff6b858a465e5d0a19ebf611ba1c9be38d10
--- /dev/null
+++ b/tests/block_bytecode/lola_stochastic_block.mod
@@ -0,0 +1,26 @@
+/* Stochastic version of block decomposed LOLA model.
+   Check that policy functions are the same as in non-block version. */
+
+@#define deterministic = false
+@#define block = true
+@#define mfs = 0
+@#include "lola_common.inc"
+
+[~, state_reorder] = sort(oo_.dr.state_var);
+
+ref = load('lola_stochastic_results.mat');
+
+[~, ref_state_reorder] = sort(ref.oo_.dr.state_var);
+
+/* NB: With block, the rows of ghx and ghu are in declaration order (and not in
+       DR-order as in non-block mode) */
+
+if max(max(abs(oo_.dr.ghx(:, state_reorder) - ref.oo_.dr.ghx(ref.oo_.dr.inv_order_var, ref_state_reorder)))) > 2e-9
+    error('Error in ghx')
+end
+
+if max(max(abs(oo_.dr.ghu - ref.oo_.dr.ghu(ref.oo_.dr.inv_order_var, :)))) > 3e-8
+    error('Error in ghu')
+end
+
+
diff --git a/tests/deterministic_simulations/lola_solve_one_boundary_mfs1.mod b/tests/deterministic_simulations/lola_solve_one_boundary_mfs1.mod
deleted file mode 100644
index 2554539c56a44748aabdb6b94f8164f9e9a0c8c7..0000000000000000000000000000000000000000
--- a/tests/deterministic_simulations/lola_solve_one_boundary_mfs1.mod
+++ /dev/null
@@ -1,987 +0,0 @@
-// Tests option mfs=1 with block
-
-load lola_data.mat
-
-% ====================================================
-% declarations var -- varexo -- para   
-% ====================================================
-
-@#define nbr_work_generations=9
-@#define nbr_early_generations=2
-@#define nbr_generations=16
-
-parameters
-length_period age_early;
-
-length_period=5;
-age_early=55;
-
-@#define wt=[1]
-
-@#define wg=0:nbr_work_generations-1
-@#define ag=0:nbr_generations-1
-@#define fwg=0:nbr_work_generations-nbr_early_generations-1
-@#define nbwg=1:nbr_work_generations-1
-@#define nbg=1:nbr_generations-1
-@#define rg=nbr_work_generations:nbr_generations-1
-@#define erg=nbr_work_generations-nbr_early_generations:nbr_work_generations-1
-@#define endg=[nbr_generations-1]
-@#define endw=[nbr_work_generations-1]
-
-@#for i in wg
-var 
-n@{i} u@{i} Omega@{i} w@{i} dWHN@{i} dWFN@{i} 
-n@{i}_f u@{i}_f Omega@{i}_f w@{i}_f dWFN@{i}_f
-i@{i} lambda@{i} i@{i}_f lambda@{i}_f eta@{i};
-parameters
-Du@{i} Dn@{i} h@{i} h@{i}_f chi@{i} eta@{i}b; 
-varexo
-eps_eta@{i};
-@#endfor
-
-@#for i in ag
-var 
-c@{i} s@{i} P@{i} P@{i}_f;
-varexo
-beta@{i} beta@{i}_f PD@{i};
-@#endfor
-
-@#for i in erg
-var 
-WE@{i} De_@{i};
-parameters
-De_@{i}b; 
-varexo
-eps_De_@{i};
-@#endfor
-
-var   
-wb wb_f 
-Omega Omega_f Omega_hf 
-V M qq p 
-N N_f 
-Q RR H K Y gdp nx FH pi
-ct st wshare rr 
-gamma mc phii D DH DF X bs bsY P00_f
-
-rhou rhoe rhol tauw tauc tauf tauk g 
-TFP gh rrb
-theta tau1 om1 om2 om2s Ds phijs
-
-DepRatio DepRatio_n DepRatio_d ZARA Ptot Ptot_f sleep du de dl inA inB in 
-NBR NBRY NBR2 tauw2 tauf2 tauc2
-PensCorr_L PensCorr_F;
-
-parameters 
-rho phi delta alpha beta ann 
-fc nu aa 
-
-rhoub rhoeb rholb tauwb taucb taufb taukb gb  
-TFPb ghb rrbb
-thetab tau1b om1b om2b om2sb Dsb phijsb
-
-NBRYb bsY_iss;
-
-varexo
-P00 P00_foP00
-
-eps_rhol eps_tauw eps_tauf eps_tauc eps_tauk
-eps_rhoe eps_rhou eps_TFP eps_gh eps_theta eps_g
-eps_Ds eps_phijs eps_PensCorr_L eps_PensCorr_F;
-
-
-% ============================================================
-% initialization
-% ============================================================
-
-@#for i in wg
-set_param_value('Du@{i}',Du@{i});
-set_param_value('Dn@{i}',Dn@{i});
-set_param_value('h@{i}',h@{i});
-set_param_value('h@{i}_f',h@{i}_f);
-set_param_value('chi@{i}',chi@{i});
-set_param_value('eta@{i}b',eta@{i}b);
-@#endfor
-
-@#for i in erg
-set_param_value('De_@{i}b',De_@{i}b);
-@#endfor
-
-set_param_value('rho',rho);
-set_param_value('phi',phi);
-set_param_value('delta',delta);
-set_param_value('alpha',alpha);
-set_param_value('beta',beta);
-set_param_value('ann',ann);
-set_param_value('fc',fc);
-set_param_value('nu',nu);
-set_param_value('aa',aa);
-
-set_param_value('rhoub',rhoub);
-set_param_value('rhoeb',rhoeb);
-set_param_value('rholb',rholb);
-set_param_value('tauwb',tauwb);
-set_param_value('taucb',taucb);
-set_param_value('taufb',taufb);
-set_param_value('taukb',taukb);
-set_param_value('gb',gb);
-
-set_param_value('TFPb',TFPb);
-set_param_value('ghb',ghb);
-set_param_value('rrbb',rrbb);
-
-set_param_value('thetab',thetab);
-set_param_value('tau1b',tau1b);
-set_param_value('om1b',om1b);
-set_param_value('om2b',om2b);
-set_param_value('om2sb',om2sb);
-set_param_value('Dsb',Dsb);
-set_param_value('phijsb',phijsb);
-
-set_param_value('bsY_iss',bsY_iss);
-
-NBRYb=NBR_iss/(phii_iss*gdp_iss);
-
-
-% =======================================================
-model(block, mfs=1);
-% ======================================================
-
-%  Labor Market Variables in the home country
-%  ------------------------------------------ 
-
-@#for i in fwg
-0=lambda@{i};
-@#endfor
-
-@#for i in wg
-1=n@{i}+u@{i}+i@{i};
-@#endfor   
-                            
-i0=lambda0;
-@#for i in nbwg
-i@{i}=lambda@{i-1}(-1)+lambda@{i}*(1-lambda@{i-1}(-1));
-@#endfor 
-
-P0=beta0*P00+PD0;
-@#for i in nbg
-P@{i}=beta@{i}*P@{i-1}(-1)+PD@{i};
-@#endfor
-
-Omega0=P0;
-@#for i in nbwg
-Omega@{i}=(1-lambda@{i})*( 1-lambda@{i-1}(-1)-(1-chi@{i})*n@{i-1}(-1))*P@{i};
-@#endfor 
-
-n0=p;
-@#for i in nbwg
-n@{i}=(1-lambda@{i})*((1-p)*(1-chi@{i})*n@{i-1}(-1)+p*(1-lambda@{i-1}(-1)));
-@#endfor 
-
-N=
-@#for i in wg
-+n@{i}*P@{i}
-@#endfor
-;
-
-%  Labor Market Variables in the foreign country
-%  --------------------------------------------- 
-
-@#for i in wg
-1=n@{i}_f+u@{i}_f+i@{i}_f;
-@#endfor    
-                            
-i0_f=lambda0_f;
-@#for i in nbwg
-i@{i}_f=lambda@{i-1}_f(-1)+lambda@{i}_f*(1-lambda@{i-1}_f(-1));
-@#endfor
-
-% -----------  reproduction cross-border   --------------------
-
-P00_f=P00_foP00*P00;
-
-P0_f=beta0_f*P00_f;
-@#for i in nbg
-P@{i}_f=beta@{i}_f*P@{i-1}_f(-1);
-@#endfor 
-
-Omega0_f=P0_f;
-@#for i in nbwg
-Omega@{i}_f=(1-lambda@{i}_f)*(1-lambda@{i-1}_f(-1)-(1-chi@{i})*n@{i-1}_f(-1))*P@{i}_f;
-@#endfor 
-
-n0_f=p;
-@#for i in nbwg
-n@{i}_f=(1-lambda@{i}_f)*((1-p)*(1-chi@{i})*n@{i-1}_f(-1)+p*(1-lambda@{i-1}_f(-1)));
-@#endfor
-
-N_f=
-@#for i in wg
-+n@{i}_f*P@{i}_f
-@#endfor
-;
-
-%  Matching
-% ----------
-
-Omega=
-@#for i in wg
-+Omega@{i}
-@#endfor
-;
-
-Omega_f=
-@#for i in wg
-+Omega@{i}_f
-@#endfor
-;
-
-Omega_hf=Omega+Omega_f;
- 
-M=V*Omega_hf/(V^nu+Omega_hf^nu)^(1/nu);
-
-qq=M/V;
-p=M/Omega_hf;
-
-%  Flow Budget Constraints (no bequests)
-% --------------------------------------
-
-rhou*w0*u0+ (1-tauw)*w0*n0 = (1+tauc)*c0+s0;
-@#for i in nbwg
-(1+rr*(1-tauk))/(beta@{i})^ann*s@{i-1}(-1)/(1+gh)+rhou*w@{i}*u@{i}+rhoe*w@{i}*i@{i}+(1-tauw)*w@{i}*n@{i}=(1+tauc)*c@{i}+s@{i};
-@#endfor
-@#for i in rg
-(1+rr*(1-tauk))/(beta@{i})^ann*s@{i-1}(-1)/(1+gh)+rhol*wb=(1+tauc)*c@{i}+s@{i};
-@#endfor
-@#for i in endg
-s@{i}=0;
-@#endfor
-
-wb=
-@#for i in wg
-+w@{i}/@{nbr_work_generations}
-@#endfor
-;
-
-%  Euler Conditions
-% ------------------
-
-@#for i in nbg
-1/(1+tauc)/c@{i-1}=beta*(1+rr(+1)*(1-tauk(+1)))/(1+tauc(+1))/c@{i}(+1)*(beta@{i})^(1-ann)/(1+gh);
-@#endfor
-
-
-%  Optimal Participation Rates (Early Retirement)
-%  ----------------------------------------------
-
-@#for i in erg
-WE@{i} = 0;
-@#endfor 
-
-@#for i in erg
-@#if i in endw
-WE@{i} = ((De_@{i}b*i@{i}^(phi-1))+Du@{i}+(rhoe-rhou)*w@{i}/(1+tauc)/c@{i})*(1-i@{i})-((1-tauw-rhou)*w@{i}/(1+tauc)/c@{i}-(Dn@{i}-Du@{i}))*n@{i};
-@#else  
-WE@{i} = ((De_@{i}b*i@{i}^(phi-1))+Du@{i}+(rhoe-rhou)*w@{i}/(1+tauc)/c@{i})*(1-i@{i})-((1-tauw-rhou)*w@{i}/(1+tauc)/c@{i}-(Dn@{i}-Du@{i}))*n@{i}+ beta*beta@{i+1}*WE@{i+1};
-@#endif
-@#endfor  
-
-%  Household Surplus
-% -------------------
-
-@#for i in wg
-@#if i in endw
-dWHN@{i} = w@{i}*(1-tauw-rhou)/(1+tauc)-(Dn@{i}-Du@{i})*c@{i};
-@#else  
-dWHN@{i} = w@{i}*(1-tauw-rhou)/(1+tauc)-(Dn@{i}-Du@{i})*c@{i} + beta*beta@{i+1}*c@{i}/c@{i+1}(+1)*dWHN@{i+1}(+1)*(1-p(+1))*(1-chi@{i+1})*(1-lambda@{i+1}(+1));
-@#endif
-@#endfor    
-
-%  Foreign household
-%  ------------------------
-%  participation and wages
-% ........................
-
-@#for i in wg
-lambda@{i}=lambda@{i}_f;
-w@{i}_f=w@{i}; 
-@#endfor
-
-wb_f=
-@#for i in wg
-+w@{i}_f/@{nbr_work_generations}
-@#endfor
-;
-
-%  Firm's Behavior
-% -------------------
-
-H=
-@#for i in wg
-+h@{i}*n@{i}*P@{i}+h@{i}_f*n@{i}_f*P@{i}_f
-@#endfor
-;
-
-wshare=(1+tauf)*(
-@#for i in wg
-+w@{i}*n@{i}*P@{i}+w@{i}_f*n@{i}_f*P@{i}_f
-@#endfor
-)/(phii*gdp);
-
-Y=TFP*H^(1-alpha)*(K)^alpha;
-gdp= TFP*H^(1-alpha)*(K)^alpha-aa*V/phii-fc/phii;
-pi = phii*gdp - wshare*phii*gdp - (rr+delta)*K(-1)/(1+gh);
-(rr(+1)+delta)/(1+rr(+1)*(1-tauk(+1))) = mc*TFP*alpha*(H/(K))^(1-alpha);
-FH=TFP*(1-alpha)*((K)/H)^alpha;
-
-RR=1+rr*(1-tauk);
-rr=rrb+tau1*(exp(bsY_iss-bsY)-1);
-
-%  Firm's Surplus
-%  ---------------------
-
-@#for i in wg
-@#if i in endw
-dWFN@{i} = h@{i}*mc*FH-(1+tauf)*w@{i};
-dWFN@{i}_f = h@{i}_f*mc*FH-(1+tauf)*w@{i}_f;
-@#else  
-dWFN@{i} = h@{i}*mc*FH-(1+tauf)*w@{i} + beta@{i+1}/RR(+1)*dWFN@{i+1}(+1)*(1-chi@{i+1})*(1-lambda@{i+1}(+1))*(1+gh);
-dWFN@{i}_f = h@{i}_f*mc*FH-(1+tauf)*w@{i}_f + beta@{i+1}_f/RR(+1)*dWFN@{i+1}_f(+1)*(1-chi@{i+1})*(1-lambda@{i+1}_f(+1))*(1+gh);
-@#endif
-@#endfor
-
-%  Free Entry Condition
-%  ---------------------
-
-aa=qq/Omega_hf*(
-@#for i in wg
-+Omega@{i}*dWFN@{i}+Omega@{i}_f*dWFN@{i}_f
-@#endfor
-);
-
-%  Wage Determination (Rent Sharing)
-%  -----------------------------------
-
-@#for i in wg
-(1-eta@{i})*dWHN@{i}  =  eta@{i}*((1-tauw)/(1+tauf)/(1+tauc))*dWFN@{i};
-@#endfor
-
-%  Equilibrium Conditions
-%  ----------------------
-
-ct=
-@#for i in ag
-+c@{i}*P@{i}
-@#endfor
-;
-
-st=
-@#for i in ag
-+s@{i}*P@{i}
-@#endfor
-;
-
-%  Non-Arbitrage condition (physical capital-shares)
-% .........................
-
-Q(+1)+pi(+1)=(1+rr(+1))*Q/(1+gh);
-
-%  New Open Economy Macroeconomics (NOEM)
-%  ---------------------------------------
-
-phii=mc/theta;
-D= ct + K-(1-delta)*K(-1)/(1+gh)  + g*gdp*phii + fc+aa*V;
-DH=(1/om1*phii)^(1/(rho-1))*D;
-X=(1/om2s*phii/gamma)^(1/(rho-1))*Ds;
-DF=(1/om2*gamma*phijs)^(1/(rho-1))*D;
-nx=phii*X-phijs*gamma*DF;     
-bsY=bs/(phii*gdp);
-Y=DH+X;
-phii*gdp=ct + K-(1-delta)*K(-1)/(1+gh) + g*gdp*phii+nx;
-
-st=K+Q +bs ;
-
-%  Policies
-%  ----------
-
-rhou=rhoub*eps_rhou;
-rhoe=rhoeb*eps_rhoe;
-rhol=rholb*eps_rhol;
-g=gb*eps_g;
-
-@#for i in erg
-De_@{i}=De_@{i}b*eps_De_@{i};
-@#endfor
-
-TFP=TFPb*eps_TFP;
-gh=ghb*eps_gh;
-
-@#for i in wg
-eta@{i}=eta@{i}b*eps_eta@{i};
-@#endfor
-
-rrb=rrbb;
-
-theta=thetab*eps_theta;
-tau1=tau1b;
-om1=om1b;
-om2=om2b;
-om2s=om2sb;
-Ds=Dsb*eps_Ds;
-phijs=phijsb*eps_phijs;
-
-
-% ----------- RefDR scenario 
-
-DepRatio_n=
-@#for i in rg
-+P@{i}
-@#endfor
-;
-DepRatio_d=
-@#for i in wg
-+P@{i}
-@#endfor
-;
-
-DepRatio=DepRatio_n/DepRatio_d;
-
-ZARA=age_early+length_period*(
-@#for i in erg
-+1-i@{i}
-@#endfor  
-);
-
-% ----------- WGEM 
-
-Ptot=
-@#for i in ag
-+P@{i}
-@#endfor
-;
-Ptot_f=
-@#for i in ag
-+P@{i}_f
-@#endfor
-;
-
-sleep=(1+rr)*(
-@#for i in nbg
-+1/beta@{i}*(1-1/beta@{i}^(ann-1))*s@{i-1}(-1)*P@{i}
-@#endfor  
-)/(1+gh);
-
-du=rhou*(
-@#for i in wg
-+w@{i}*u@{i}*P@{i}
-@#endfor  
-);
-
-de=rhoe*(
-@#for i in erg
-+w@{i}*i@{i}*P@{i}+w@{i}_f*i@{i}_f*P@{i}_f
-@#endfor  
-);
-
-dl=
-@#for i in rg
-+rhol*wb*PensCorr_L*P@{i}+rhol*wb_f*(N_f/(N+N_f))*PensCorr_F*P@{i}_f
-@#endfor  
-;
-
-PensCorr_L=eps_PensCorr_L; 
-PensCorr_F=eps_PensCorr_F;   
-
-inA=(tauw+tauf)*(
-@#for i in wg
-+n@{i}*w@{i}*P@{i}+n@{i}_f*w@{i}_f*P@{i}_f
-@#endfor  
-);
-
-inB=tauk*rr*(
-@#for i in nbg
-+1/beta@{i}^ann*s@{i-1}(-1)*P@{i}
-@#endfor  
-)/(1+gh);
-
-in=tauc*ct+inA+inB+sleep;
-NBR=g*phii*gdp+(du+de+dl)-(in);
-NBRY=NBR/(phii*gdp);
-
-% ----------- WGEM Adjustment variable ---------------
-
-tauf2=tauf;
-tauw2=tauw;
-NBR2=NBR;
-tauc2=tauc;
-tauf=taufb*eps_tauf;
-%----- WGEM: adjustment through tauc
-tauc=taucb*eps_tauc;
-%----- WGEM: adjustment through tauk
-tauk=taukb*eps_tauk;
-%----- WGEM: adjustment through tauw
-tauw=tauwb*eps_tauw;
-
-end;
-
-%==================================
-initval;
-%==================================
-
-@#for i in wg
-n@{i}=n@{i}_iss;
-n@{i}_f=n@{i}_f_iss;
-u@{i}=u@{i}_iss;
-u@{i}_f=u@{i}_f_iss;
-Omega@{i}=Omega@{i}_iss;
-Omega@{i}_f=Omega@{i}_f_iss;
-w@{i}=w@{i}_iss;
-w@{i}_f=w@{i}_f_iss;
-dWHN@{i}=dWHN@{i}_iss;
-dWFN@{i}=dWFN@{i}_iss;
-dWFN@{i}_f=dWFN@{i}_f_iss;
-eps_eta@{i}=eps_eta@{i}_iss;
-eta@{i}=eta@{i}b*eps_eta@{i};
-@#endfor
-
-@#for i in erg
-i@{i}=i@{i}_iss;
-lambda@{i}=lambda@{i}_iss;
-i@{i}_f=i@{i}_f_iss;
-lambda@{i}_f=lambda@{i}_f_iss;
-WE@{i}=0;
-eps_De_@{i}=eps_De_@{i}_iss;
-De_@{i}=De_@{i}b*eps_De_@{i}_iss;
-@#endfor
-
-@#for i in fwg
-i@{i}=0;
-lambda@{i}=0;
-i@{i}_f=0;
-lambda@{i}_f=0;
-@#endfor
-
-@#for i in ag
-@#if i in endg
-s@{i}=0;
-@#else
-s@{i}=s@{i}_iss;
-@#endif
-@#endfor
-
-@#for i in ag
-beta@{i}=beta@{i}_iss;
-beta@{i}_f=beta@{i}_f_iss;
-PD@{i}=PD@{i}_iss;
-c@{i}=c@{i}_iss;
-P@{i}=P@{i}_iss;
-P@{i}_f=P@{i}_f_iss;
-@#endfor
-
-wb          =   wb_iss;
-wb_f        =   wb_f_iss;
-Omega       =   Omega_iss;
-Omega_f     =   Omega_f_iss;
-Omega_hf    =   Omega_hf_iss;
-V        	=	V_iss	;
-M        	=	M_iss	;
-qq       	=	qq_iss	;
-p        	=	p_iss	;
-N        	=	N_iss	;
-N_f      	=	N_f_iss	;
-Q        	=	Q_iss	;
-RR       	=	RR_iss	;
-H        	=	H_iss	;
-K        	=	K_iss	;
-Y        	=	Y_iss	;
-gdp      	=	gdp_iss	;
-nx       	=	nx_iss	;
-FH       	=	FH_iss	;
-pi       	=	pi_iss	;
-ct       	=	ct_iss	;
-st       	=	st_iss	;
-wshare   	=	wshare_iss	;
-rr       	=	rr_iss	;
-
-gamma    	=	gamma_iss	;
-mc       	=	mc_iss	;
-phii     	=	phii_iss	;
-D        	=	D_iss	;
-DH       	=	DH_iss	;
-DF       	=	DF_iss	;
-X        	=	X_iss	;
-bs       	=	bs_iss	;
-bsY      	=	bsY_iss	;
-P00_f       =   P00_f_iss;
-
-eps_rhol=eps_rhol_iss;
-eps_rhoe=eps_rhoe_iss;
-eps_rhou=eps_rhou_iss;
-rhou=rhoub*eps_rhou_iss;
-rhoe=rhoeb*eps_rhoe_iss;
-rhol=rholb*eps_rhol_iss;
-eps_tauc=eps_tauc_iss;
-eps_tauk=eps_tauk_iss;
-eps_tauw=eps_tauw_iss;
-eps_tauf=eps_tauf_iss;
-tauw=tauwb*eps_tauw;
-tauc=taucb*eps_tauc;
-tauf=taufb*eps_tauf;
-tauk=taukb*eps_tauk;
-eps_theta=eps_theta_iss;
-eps_gh=eps_gh_iss;
-eps_TFP=eps_TFP_iss;
-eps_g=eps_g_iss;
-g=gb*eps_g_iss;
-TFP=TFPb*eps_TFP_iss;
-gh=ghb*eps_gh_iss;
-theta=thetab*eps_theta_iss;
-
-rrb=rrbb;
-tau1=tau1b;
-om1=om1b;
-om2=om2b;
-om2s=om2sb;
-
-eps_Ds=1;
-eps_phijs=1;
-Ds=Dsb*eps_Ds;
-phijs=phijsb*eps_phijs;
-
-eps_PensCorr_F=eps_PensCorr_F_iss;
-eps_PensCorr_L=eps_PensCorr_L_iss;
-PensCorr_F=eps_PensCorr_F_iss;
-PensCorr_L=eps_PensCorr_L_iss;
-
-P00	        =	P00_iss	;
-P00_foP00	=	P00_foP00_iss	;
-
-DepRatio_n=
-@#for i in rg
-+P@{i}
-@#endfor
-;
-DepRatio_d=
-@#for i in wg
-+P@{i}
-@#endfor
-;
-
-DepRatio=DepRatio_n/DepRatio_d;
-
-ZARA=age_early+length_period*(
-@#for i in erg
-+1-i@{i}
-@#endfor  
-);
-
-Ptot=Ptot_iss;
-Ptot_f=Ptot_f_iss;
-sleep=sleep_iss;
-du=du_iss;
-de=de_iss;
-dl=dl_iss;
-
-inA=inA_iss;%(tauw+tauf)*(
-%@#for i in wg
-%+n@{i}*w@{i}*P@{i}+n@{i}_f*w@{i}_f*P@{i}_f
-%@#endfor  
-%);
-
-inB=inB_iss;%tauk*rr*(
-%@#for i in nbg
-%+1/beta@{i}^ann*s@{i-1}*P@{i}
-%@#endfor  
-%)/(1+gh);
-
-in=in_iss;%tauc*ct+inA+inB+sleep;
-NBR=NBR_iss;%g*phii*gdp+(du+de+dl)-(in);
-NBRY=NBR/(phii*gdp);
-NBR2=NBR;
-tauf2=tauf;
-tauw2=tauw;
-tauc2=tauc;
-
-end;
-
-%========================================================
-% compute initial steady state and check eigenvalues 
-%========================================================
-
-resid;
-steady(solve_algo=3);
-check;
-
-%========================================================
-endval;
-%========================================================
-
-@#for i in wg
-n@{i}=n@{i}_fss;
-n@{i}_f=n@{i}_f_fss;
-u@{i}=u@{i}_fss;
-u@{i}_f=u@{i}_f_fss;
-Omega@{i}=Omega@{i}_fss;
-Omega@{i}_f=Omega@{i}_f_fss;
-w@{i}=w@{i}_fss;
-w@{i}_f=w@{i}_f_fss;
-dWHN@{i}=dWHN@{i}_fss;
-dWFN@{i}=dWFN@{i}_fss;
-dWFN@{i}_f=dWFN@{i}_f_fss;
-eps_eta@{i}=eps_eta@{i}_fss;
-eta@{i}=eta@{i}b*eps_eta@{i};
-@#endfor
-
-@#for i in erg
-i@{i}=i@{i}_fss;
-lambda@{i}=lambda@{i}_fss;
-i@{i}_f=i@{i}_f_fss;
-lambda@{i}_f=lambda@{i}_f_fss;
-WE@{i}=0;
-eps_De_@{i}=eps_De_@{i}_fss;
-De_@{i}=De_@{i}b*eps_De_@{i}_fss;
-@#endfor
-
-@#for i in fwg
-i@{i}=0;
-lambda@{i}=0;
-i@{i}_f=0;
-lambda@{i}_f=0;
-@#endfor
-
-@#for i in ag
-@#if i in endg
-s@{i}=0;
-@#else
-s@{i}=s@{i}_fss;
-@#endif
-@#endfor
-
-@#for i in ag
-beta@{i}=beta@{i}_fss;
-beta@{i}_f=beta@{i}_f_fss;
-PD@{i}=PD@{i}_fss;
-c@{i}=c@{i}_fss;
-P@{i}=P@{i}_fss;
-P@{i}_f=P@{i}_f_fss;
-@#endfor
-
-wb          =   wb_fss;
-wb_f        =   wb_f_fss;
-Omega       =   Omega_fss;
-Omega_f     =   Omega_f_fss;
-Omega_hf    =   Omega_hf_fss;
-V        	=	V_fss	;
-M        	=	M_fss	;
-qq       	=	qq_fss	;
-p        	=	p_fss	;
-N        	=	N_fss	;
-N_f      	=	N_f_fss	;
-Q        	=	Q_fss	;
-RR       	=	RR_fss	;
-H        	=	H_fss	;
-K        	=	K_fss	;
-Y        	=	Y_fss	;
-gdp      	=	gdp_fss	;
-nx       	=	nx_fss	;
-FH       	=	FH_fss	;
-pi       	=	pi_fss	;
-ct       	=	ct_fss	;
-st       	=	st_fss	;
-wshare   	=	wshare_fss	;
-rr       	=	rr_fss	;
-
-gamma    	=	gamma_fss	;
-mc       	=	mc_fss	;
-phii     	=	phii_fss	;
-D        	=	D_fss	;
-DH       	=	DH_fss	;
-DF       	=	DF_fss	;
-X        	=	X_fss	;
-bs       	=	bs_fss	;
-bsY      	=	bsY_fss	;
-P00_f       =   P00_f_fss;
-
-eps_rhol=eps_rhol_fss;
-eps_rhoe=eps_rhoe_fss;
-eps_rhou=eps_rhou_fss;
-rhou=rhoub*eps_rhou_fss;
-rhoe=rhoeb*eps_rhoe_fss;
-rhol=rholb*eps_rhol_fss;
-eps_tauc=eps_tauc_fss;
-eps_tauk=eps_tauk_fss;
-eps_tauw=eps_tauw_fss;
-eps_tauf=eps_tauf_fss;
-tauw=tauwb*eps_tauw;
-tauc=taucb*eps_tauc;
-tauf=taufb*eps_tauf;
-tauk=taukb*eps_tauk;
-eps_theta=eps_theta_fss;
-eps_gh=eps_gh_fss;
-eps_TFP=eps_TFP_fss;
-eps_g=eps_g_fss;
-g=gb*eps_g_fss;
-TFP=TFPb*eps_TFP_fss;
-gh=ghb*eps_gh_fss;
-theta=thetab*eps_theta_fss;
-
-rrb=rrbb;
-tau1=tau1b;
-om1=om1b;
-om2=om2b;
-om2s=om2sb;
-
-eps_Ds=1;
-eps_phijs=1;
-Ds=Dsb*eps_Ds;
-phijs=phijsb*eps_phijs;
-
-eps_PensCorr_F=eps_PensCorr_F_fss;
-eps_PensCorr_L=eps_PensCorr_L_fss;
-PensCorr_F=eps_PensCorr_F_fss;
-PensCorr_L=eps_PensCorr_L_fss;
-
-P00	        =	P00_fss	;
-P00_foP00	=	P00_foP00_fss	;
-
-DepRatio_n=
-@#for i in rg
-+P@{i}
-@#endfor
-;
-DepRatio_d=
-@#for i in wg
-+P@{i}
-@#endfor
-;
-
-DepRatio=DepRatio_n/DepRatio_d;
-
-ZARA=age_early+length_period*(
-@#for i in erg
-+1-i@{i}
-@#endfor  
-);
-
-Ptot=
-@#for i in ag
-+P@{i}
-@#endfor
-;
-Ptot_f=
-@#for i in ag
-+P@{i}_f
-@#endfor
-;
-
-sleep=(1+rr)*(
-@#for i in nbg
-+1/beta@{i}*(1-1/beta@{i}^(ann-1))*s@{i-1}*P@{i}
-@#endfor  
-)/(1+gh);
-
-du=rhou*(
-@#for i in wg
-+w@{i}*u@{i}*P@{i}
-@#endfor  
-);
-
-de=rhoe*(
-@#for i in erg
-+w@{i}*i@{i}*P@{i}+w@{i}_f*i@{i}_f*P@{i}_f
-@#endfor  
-);
-
-dl=
-@#for i in rg
-+rhol*wb*PensCorr_L*P@{i}+rhol*wb_f*(N_f/(N+N_f))*PensCorr_F*P@{i}_f
-@#endfor  
-;
-
-inA=(tauw+tauf)*(
-@#for i in wg
-+n@{i}*w@{i}*P@{i}+n@{i}_f*w@{i}_f*P@{i}_f
-@#endfor  
-);
-
-inB=tauk*rr*(
-@#for i in nbg
-+1/beta@{i}^ann*s@{i-1}*P@{i}
-@#endfor  
-)/(1+gh);
-
-in=tauc*ct+inA+inB+sleep;
-NBR=g*phii*gdp+(du+de+dl)-(in);
-NBRY=NBR/(phii*gdp);
-NBR2=NBR;
-tauf2=tauf;
-tauw2=tauw;
-tauc2=tauc;
-
-end;
-
-%========================================================
-% compute final steady state and check eigenvalues 
-%========================================================
-
-resid;
-steady(solve_algo=3);
-check;
-
-
-% ===================================================
-shocks;     
-% ===================================================
-
-var P00;																																																			
-   periods 1:99;																																																			
-   values (se_P00);																																																			
- 
-@#for i in nbg
-var beta@{i};																																																			
-   periods  1:99;																																																			
-   values (se_beta@{i});
-var beta@{i}_f;																																																			
-   periods  1:99;																																																			
-   values (se_beta@{i}_f);
-var PD@{i};																																																			
-   periods  1:99;																																																			
-   values (se_PD@{i});
-@#endfor  
-																																																			
-var P00_foP00;																																																			
-   periods 1:99;																																																			
-   values (se_P00_foP00);	 																																																
-
-var eps_g;																																																			
-   periods 1:99;																																																			
-   values (se_eps_g);	
-
-var eps_PensCorr_F;																																																			
-   periods 1:99;																																																			
-   values (se_eps_PensCorr_F);	
-
-var eps_PensCorr_L;																																																			
-   periods 1:99;																																																			
-   values (se_eps_PensCorr_L);	
-
-end;
-
-% *******************************************
-% Numerical Simulation, Control Parameters
-% *******************************************
-
-model_info;
-
-perfect_foresight_setup(periods=125);
-perfect_foresight_solver(maxit=100);
-
-if ~oo_.deterministic_simulation.status
-   error('Perfect foresight simulation failed')
-end
-
-mfs0=load('lola_solve_one_boundary_results');
-
-if max(max(oo_.endo_simul-mfs0.oo_.endo_simul)) > options_.dynatol.x
-   error('Inconsistency with mfs=0')
-end
diff --git a/tests/deterministic_simulations/lola_solve_one_boundary_mfs2.mod b/tests/deterministic_simulations/lola_solve_one_boundary_mfs2.mod
deleted file mode 100644
index 7d569e20a1370a08c95bcde81c5f915a99712e8a..0000000000000000000000000000000000000000
--- a/tests/deterministic_simulations/lola_solve_one_boundary_mfs2.mod
+++ /dev/null
@@ -1,985 +0,0 @@
-// Tests option mfs=2 with block
-
-load lola_data.mat
-
-% ====================================================
-% declarations var -- varexo -- para   
-% ====================================================
-
-@#define nbr_work_generations=9
-@#define nbr_early_generations=2
-@#define nbr_generations=16
-
-parameters
-length_period age_early;
-
-length_period=5;
-age_early=55;
-
-@#define wt=[1]
-
-@#define wg=0:nbr_work_generations-1
-@#define ag=0:nbr_generations-1
-@#define fwg=0:nbr_work_generations-nbr_early_generations-1
-@#define nbwg=1:nbr_work_generations-1
-@#define nbg=1:nbr_generations-1
-@#define rg=nbr_work_generations:nbr_generations-1
-@#define erg=nbr_work_generations-nbr_early_generations:nbr_work_generations-1
-@#define endg=[nbr_generations-1]
-@#define endw=[nbr_work_generations-1]
-
-@#for i in wg
-var 
-n@{i} u@{i} Omega@{i} w@{i} dWHN@{i} dWFN@{i} 
-n@{i}_f u@{i}_f Omega@{i}_f w@{i}_f dWFN@{i}_f
-i@{i} lambda@{i} i@{i}_f lambda@{i}_f eta@{i};
-parameters
-Du@{i} Dn@{i} h@{i} h@{i}_f chi@{i} eta@{i}b; 
-varexo
-eps_eta@{i};
-@#endfor
-
-@#for i in ag
-var 
-c@{i} s@{i} P@{i} P@{i}_f;
-varexo
-beta@{i} beta@{i}_f PD@{i};
-@#endfor
-
-@#for i in erg
-var 
-WE@{i} De_@{i};
-parameters
-De_@{i}b; 
-varexo
-eps_De_@{i};
-@#endfor
-
-var   
-wb wb_f 
-Omega Omega_f Omega_hf 
-V M qq p 
-N N_f 
-Q RR H K Y gdp nx FH pi
-ct st wshare rr 
-gamma mc phii D DH DF X bs bsY P00_f
-
-rhou rhoe rhol tauw tauc tauf tauk g 
-TFP gh rrb
-theta tau1 om1 om2 om2s Ds phijs
-
-DepRatio DepRatio_n DepRatio_d ZARA Ptot Ptot_f sleep du de dl inA inB in 
-NBR NBRY NBR2 tauw2 tauf2 tauc2
-PensCorr_L PensCorr_F;
-
-parameters 
-rho phi delta alpha beta ann 
-fc nu aa 
-
-rhoub rhoeb rholb tauwb taucb taufb taukb gb  
-TFPb ghb rrbb
-thetab tau1b om1b om2b om2sb Dsb phijsb
-
-NBRYb bsY_iss;
-
-varexo
-P00 P00_foP00
-
-eps_rhol eps_tauw eps_tauf eps_tauc eps_tauk
-eps_rhoe eps_rhou eps_TFP eps_gh eps_theta eps_g
-eps_Ds eps_phijs eps_PensCorr_L eps_PensCorr_F;
-
-
-% ============================================================
-% initialization
-% ============================================================
-
-@#for i in wg
-set_param_value('Du@{i}',Du@{i});
-set_param_value('Dn@{i}',Dn@{i});
-set_param_value('h@{i}',h@{i});
-set_param_value('h@{i}_f',h@{i}_f);
-set_param_value('chi@{i}',chi@{i});
-set_param_value('eta@{i}b',eta@{i}b);
-@#endfor
-
-@#for i in erg
-set_param_value('De_@{i}b',De_@{i}b);
-@#endfor
-
-set_param_value('rho',rho);
-set_param_value('phi',phi);
-set_param_value('delta',delta);
-set_param_value('alpha',alpha);
-set_param_value('beta',beta);
-set_param_value('ann',ann);
-set_param_value('fc',fc);
-set_param_value('nu',nu);
-set_param_value('aa',aa);
-
-set_param_value('rhoub',rhoub);
-set_param_value('rhoeb',rhoeb);
-set_param_value('rholb',rholb);
-set_param_value('tauwb',tauwb);
-set_param_value('taucb',taucb);
-set_param_value('taufb',taufb);
-set_param_value('taukb',taukb);
-set_param_value('gb',gb);
-
-set_param_value('TFPb',TFPb);
-set_param_value('ghb',ghb);
-set_param_value('rrbb',rrbb);
-
-set_param_value('thetab',thetab);
-set_param_value('tau1b',tau1b);
-set_param_value('om1b',om1b);
-set_param_value('om2b',om2b);
-set_param_value('om2sb',om2sb);
-set_param_value('Dsb',Dsb);
-set_param_value('phijsb',phijsb);
-
-set_param_value('bsY_iss',bsY_iss);
-
-NBRYb=NBR_iss/(phii_iss*gdp_iss);
-
-
-% =======================================================
-model(block, mfs=2);
-% ======================================================
-
-%  Labor Market Variables in the home country
-%  ------------------------------------------ 
-
-@#for i in fwg
-0=lambda@{i};
-@#endfor
-
-@#for i in wg
-1=n@{i}+u@{i}+i@{i};
-@#endfor   
-                            
-i0=lambda0;
-@#for i in nbwg
-i@{i}=lambda@{i-1}(-1)+lambda@{i}*(1-lambda@{i-1}(-1));
-@#endfor 
-
-P0=beta0*P00+PD0;
-@#for i in nbg
-P@{i}=beta@{i}*P@{i-1}(-1)+PD@{i};
-@#endfor
-
-Omega0=P0;
-@#for i in nbwg
-Omega@{i}=(1-lambda@{i})*( 1-lambda@{i-1}(-1)-(1-chi@{i})*n@{i-1}(-1))*P@{i};
-@#endfor 
-
-n0=p;
-@#for i in nbwg
-n@{i}=(1-lambda@{i})*((1-p)*(1-chi@{i})*n@{i-1}(-1)+p*(1-lambda@{i-1}(-1)));
-@#endfor 
-
-N=
-@#for i in wg
-+n@{i}*P@{i}
-@#endfor
-;
-
-%  Labor Market Variables in the foreign country
-%  --------------------------------------------- 
-
-@#for i in wg
-1=n@{i}_f+u@{i}_f+i@{i}_f;
-@#endfor    
-                            
-i0_f=lambda0_f;
-@#for i in nbwg
-i@{i}_f=lambda@{i-1}_f(-1)+lambda@{i}_f*(1-lambda@{i-1}_f(-1));
-@#endfor
-
-% -----------  reproduction cross-border   --------------------
-
-P00_f=P00_foP00*P00;
-
-P0_f=beta0_f*P00_f;
-@#for i in nbg
-P@{i}_f=beta@{i}_f*P@{i-1}_f(-1);
-@#endfor 
-
-Omega0_f=P0_f;
-@#for i in nbwg
-Omega@{i}_f=(1-lambda@{i}_f)*(1-lambda@{i-1}_f(-1)-(1-chi@{i})*n@{i-1}_f(-1))*P@{i}_f;
-@#endfor 
-
-n0_f=p;
-@#for i in nbwg
-n@{i}_f=(1-lambda@{i}_f)*((1-p)*(1-chi@{i})*n@{i-1}_f(-1)+p*(1-lambda@{i-1}_f(-1)));
-@#endfor
-
-N_f=
-@#for i in wg
-+n@{i}_f*P@{i}_f
-@#endfor
-;
-
-%  Matching
-% ----------
-
-Omega=
-@#for i in wg
-+Omega@{i}
-@#endfor
-;
-
-Omega_f=
-@#for i in wg
-+Omega@{i}_f
-@#endfor
-;
-
-Omega_hf=Omega+Omega_f;
- 
-M=V*Omega_hf/(V^nu+Omega_hf^nu)^(1/nu);
-
-qq=M/V;
-p=M/Omega_hf;
-
-%  Flow Budget Constraints (no bequests)
-% --------------------------------------
-
-rhou*w0*u0+ (1-tauw)*w0*n0 = (1+tauc)*c0+s0;
-@#for i in nbwg
-(1+rr*(1-tauk))/(beta@{i})^ann*s@{i-1}(-1)/(1+gh)+rhou*w@{i}*u@{i}+rhoe*w@{i}*i@{i}+(1-tauw)*w@{i}*n@{i}=(1+tauc)*c@{i}+s@{i};
-@#endfor
-@#for i in rg
-(1+rr*(1-tauk))/(beta@{i})^ann*s@{i-1}(-1)/(1+gh)+rhol*wb=(1+tauc)*c@{i}+s@{i};
-@#endfor
-@#for i in endg
-s@{i}=0;
-@#endfor
-
-wb=
-@#for i in wg
-+w@{i}/@{nbr_work_generations}
-@#endfor
-;
-
-%  Euler Conditions
-% ------------------
-
-@#for i in nbg
-1/(1+tauc)/c@{i-1}=beta*(1+rr(+1)*(1-tauk(+1)))/(1+tauc(+1))/c@{i}(+1)*(beta@{i})^(1-ann)/(1+gh);
-@#endfor
-
-
-%  Optimal Participation Rates (Early Retirement)
-%  ----------------------------------------------
-
-@#for i in erg
-WE@{i} = 0;
-@#endfor 
-
-@#for i in erg
-@#if i in endw
-WE@{i} = ((De_@{i}b*i@{i}^(phi-1))+Du@{i}+(rhoe-rhou)*w@{i}/(1+tauc)/c@{i})*(1-i@{i})-((1-tauw-rhou)*w@{i}/(1+tauc)/c@{i}-(Dn@{i}-Du@{i}))*n@{i};
-@#else  
-WE@{i} = ((De_@{i}b*i@{i}^(phi-1))+Du@{i}+(rhoe-rhou)*w@{i}/(1+tauc)/c@{i})*(1-i@{i})-((1-tauw-rhou)*w@{i}/(1+tauc)/c@{i}-(Dn@{i}-Du@{i}))*n@{i}+ beta*beta@{i+1}*WE@{i+1};
-@#endif
-@#endfor  
-
-%  Household Surplus
-% -------------------
-
-@#for i in wg
-@#if i in endw
-dWHN@{i} = w@{i}*(1-tauw-rhou)/(1+tauc)-(Dn@{i}-Du@{i})*c@{i};
-@#else  
-dWHN@{i} = w@{i}*(1-tauw-rhou)/(1+tauc)-(Dn@{i}-Du@{i})*c@{i} + beta*beta@{i+1}*c@{i}/c@{i+1}(+1)*dWHN@{i+1}(+1)*(1-p(+1))*(1-chi@{i+1})*(1-lambda@{i+1}(+1));
-@#endif
-@#endfor    
-
-%  Foreign household
-%  ------------------------
-%  participation and wages
-% ........................
-
-@#for i in wg
-lambda@{i}=lambda@{i}_f;
-w@{i}_f=w@{i}; 
-@#endfor
-
-wb_f=
-@#for i in wg
-+w@{i}_f/@{nbr_work_generations}
-@#endfor
-;
-
-%  Firm's Behavior
-% -------------------
-
-H=
-@#for i in wg
-+h@{i}*n@{i}*P@{i}+h@{i}_f*n@{i}_f*P@{i}_f
-@#endfor
-;
-
-wshare=(1+tauf)*(
-@#for i in wg
-+w@{i}*n@{i}*P@{i}+w@{i}_f*n@{i}_f*P@{i}_f
-@#endfor
-)/(phii*gdp);
-
-Y=TFP*H^(1-alpha)*(K)^alpha;
-gdp= TFP*H^(1-alpha)*(K)^alpha-aa*V/phii-fc/phii;
-pi = phii*gdp - wshare*phii*gdp - (rr+delta)*K(-1)/(1+gh);
-(rr(+1)+delta)/(1+rr(+1)*(1-tauk(+1))) = mc*TFP*alpha*(H/(K))^(1-alpha);
-FH=TFP*(1-alpha)*((K)/H)^alpha;
-
-RR=1+rr*(1-tauk);
-rr=rrb+tau1*(exp(bsY_iss-bsY)-1);
-
-%  Firm's Surplus
-%  ---------------------
-
-@#for i in wg
-@#if i in endw
-dWFN@{i} = h@{i}*mc*FH-(1+tauf)*w@{i};
-dWFN@{i}_f = h@{i}_f*mc*FH-(1+tauf)*w@{i}_f;
-@#else  
-dWFN@{i} = h@{i}*mc*FH-(1+tauf)*w@{i} + beta@{i+1}/RR(+1)*dWFN@{i+1}(+1)*(1-chi@{i+1})*(1-lambda@{i+1}(+1))*(1+gh);
-dWFN@{i}_f = h@{i}_f*mc*FH-(1+tauf)*w@{i}_f + beta@{i+1}_f/RR(+1)*dWFN@{i+1}_f(+1)*(1-chi@{i+1})*(1-lambda@{i+1}_f(+1))*(1+gh);
-@#endif
-@#endfor
-
-%  Free Entry Condition
-%  ---------------------
-
-aa=qq/Omega_hf*(
-@#for i in wg
-+Omega@{i}*dWFN@{i}+Omega@{i}_f*dWFN@{i}_f
-@#endfor
-);
-
-%  Wage Determination (Rent Sharing)
-%  -----------------------------------
-
-@#for i in wg
-(1-eta@{i})*dWHN@{i}  =  eta@{i}*((1-tauw)/(1+tauf)/(1+tauc))*dWFN@{i};
-@#endfor
-
-%  Equilibrium Conditions
-%  ----------------------
-
-ct=
-@#for i in ag
-+c@{i}*P@{i}
-@#endfor
-;
-
-st=
-@#for i in ag
-+s@{i}*P@{i}
-@#endfor
-;
-
-%  Non-Arbitrage condition (physical capital-shares)
-% .........................
-
-Q(+1)+pi(+1)=(1+rr(+1))*Q/(1+gh);
-
-%  New Open Economy Macroeconomics (NOEM)
-%  ---------------------------------------
-
-phii=mc/theta;
-D= ct + K-(1-delta)*K(-1)/(1+gh)  + g*gdp*phii + fc+aa*V;
-DH=(1/om1*phii)^(1/(rho-1))*D;
-X=(1/om2s*phii/gamma)^(1/(rho-1))*Ds;
-DF=(1/om2*gamma*phijs)^(1/(rho-1))*D;
-nx=phii*X-phijs*gamma*DF;     
-bsY=bs/(phii*gdp);
-Y=DH+X;
-phii*gdp=ct + K-(1-delta)*K(-1)/(1+gh) + g*gdp*phii+nx;
-
-st=K+Q +bs ;
-
-%  Policies
-%  ----------
-
-rhou=rhoub*eps_rhou;
-rhoe=rhoeb*eps_rhoe;
-rhol=rholb*eps_rhol;
-g=gb*eps_g;
-
-@#for i in erg
-De_@{i}=De_@{i}b*eps_De_@{i};
-@#endfor
-
-TFP=TFPb*eps_TFP;
-gh=ghb*eps_gh;
-
-@#for i in wg
-eta@{i}=eta@{i}b*eps_eta@{i};
-@#endfor
-
-rrb=rrbb;
-
-theta=thetab*eps_theta;
-tau1=tau1b;
-om1=om1b;
-om2=om2b;
-om2s=om2sb;
-Ds=Dsb*eps_Ds;
-phijs=phijsb*eps_phijs;
-
-
-% ----------- RefDR scenario 
-
-DepRatio_n=
-@#for i in rg
-+P@{i}
-@#endfor
-;
-DepRatio_d=
-@#for i in wg
-+P@{i}
-@#endfor
-;
-
-DepRatio=DepRatio_n/DepRatio_d;
-
-ZARA=age_early+length_period*(
-@#for i in erg
-+1-i@{i}
-@#endfor  
-);
-
-% ----------- WGEM 
-
-Ptot=
-@#for i in ag
-+P@{i}
-@#endfor
-;
-Ptot_f=
-@#for i in ag
-+P@{i}_f
-@#endfor
-;
-
-sleep=(1+rr)*(
-@#for i in nbg
-+1/beta@{i}*(1-1/beta@{i}^(ann-1))*s@{i-1}(-1)*P@{i}
-@#endfor  
-)/(1+gh);
-
-du=rhou*(
-@#for i in wg
-+w@{i}*u@{i}*P@{i}
-@#endfor  
-);
-
-de=rhoe*(
-@#for i in erg
-+w@{i}*i@{i}*P@{i}+w@{i}_f*i@{i}_f*P@{i}_f
-@#endfor  
-);
-
-dl=
-@#for i in rg
-+rhol*wb*PensCorr_L*P@{i}+rhol*wb_f*(N_f/(N+N_f))*PensCorr_F*P@{i}_f
-@#endfor  
-;
-
-PensCorr_L=eps_PensCorr_L; 
-PensCorr_F=eps_PensCorr_F;   
-
-inA=(tauw+tauf)*(
-@#for i in wg
-+n@{i}*w@{i}*P@{i}+n@{i}_f*w@{i}_f*P@{i}_f
-@#endfor  
-);
-
-inB=tauk*rr*(
-@#for i in nbg
-+1/beta@{i}^ann*s@{i-1}(-1)*P@{i}
-@#endfor  
-)/(1+gh);
-
-in=tauc*ct+inA+inB+sleep;
-NBR=g*phii*gdp+(du+de+dl)-(in);
-NBRY=NBR/(phii*gdp);
-
-% ----------- WGEM Adjustment variable ---------------
-
-tauf2=tauf;
-tauw2=tauw;
-NBR2=NBR;
-tauc2=tauc;
-tauf=taufb*eps_tauf;
-%----- WGEM: adjustment through tauc
-tauc=taucb*eps_tauc;
-%----- WGEM: adjustment through tauk
-tauk=taukb*eps_tauk;
-%----- WGEM: adjustment through tauw
-tauw=tauwb*eps_tauw;
-
-end;
-
-%==================================
-initval;
-%==================================
-
-@#for i in wg
-n@{i}=n@{i}_iss;
-n@{i}_f=n@{i}_f_iss;
-u@{i}=u@{i}_iss;
-u@{i}_f=u@{i}_f_iss;
-Omega@{i}=Omega@{i}_iss;
-Omega@{i}_f=Omega@{i}_f_iss;
-w@{i}=w@{i}_iss;
-w@{i}_f=w@{i}_f_iss;
-dWHN@{i}=dWHN@{i}_iss;
-dWFN@{i}=dWFN@{i}_iss;
-dWFN@{i}_f=dWFN@{i}_f_iss;
-eps_eta@{i}=eps_eta@{i}_iss;
-eta@{i}=eta@{i}b*eps_eta@{i};
-@#endfor
-
-@#for i in erg
-i@{i}=i@{i}_iss;
-lambda@{i}=lambda@{i}_iss;
-i@{i}_f=i@{i}_f_iss;
-lambda@{i}_f=lambda@{i}_f_iss;
-WE@{i}=0;
-eps_De_@{i}=eps_De_@{i}_iss;
-De_@{i}=De_@{i}b*eps_De_@{i}_iss;
-@#endfor
-
-@#for i in fwg
-i@{i}=0;
-lambda@{i}=0;
-i@{i}_f=0;
-lambda@{i}_f=0;
-@#endfor
-
-@#for i in ag
-@#if i in endg
-s@{i}=0;
-@#else
-s@{i}=s@{i}_iss;
-@#endif
-@#endfor
-
-@#for i in ag
-beta@{i}=beta@{i}_iss;
-beta@{i}_f=beta@{i}_f_iss;
-PD@{i}=PD@{i}_iss;
-c@{i}=c@{i}_iss;
-P@{i}=P@{i}_iss;
-P@{i}_f=P@{i}_f_iss;
-@#endfor
-
-wb          =   wb_iss;
-wb_f        =   wb_f_iss;
-Omega       =   Omega_iss;
-Omega_f     =   Omega_f_iss;
-Omega_hf    =   Omega_hf_iss;
-V        	=	V_iss	;
-M        	=	M_iss	;
-qq       	=	qq_iss	;
-p        	=	p_iss	;
-N        	=	N_iss	;
-N_f      	=	N_f_iss	;
-Q        	=	Q_iss	;
-RR       	=	RR_iss	;
-H        	=	H_iss	;
-K        	=	K_iss	;
-Y        	=	Y_iss	;
-gdp      	=	gdp_iss	;
-nx       	=	nx_iss	;
-FH       	=	FH_iss	;
-pi       	=	pi_iss	;
-ct       	=	ct_iss	;
-st       	=	st_iss	;
-wshare   	=	wshare_iss	;
-rr       	=	rr_iss	;
-
-gamma    	=	gamma_iss	;
-mc       	=	mc_iss	;
-phii     	=	phii_iss	;
-D        	=	D_iss	;
-DH       	=	DH_iss	;
-DF       	=	DF_iss	;
-X        	=	X_iss	;
-bs       	=	bs_iss	;
-bsY      	=	bsY_iss	;
-P00_f       =   P00_f_iss;
-
-eps_rhol=eps_rhol_iss;
-eps_rhoe=eps_rhoe_iss;
-eps_rhou=eps_rhou_iss;
-rhou=rhoub*eps_rhou_iss;
-rhoe=rhoeb*eps_rhoe_iss;
-rhol=rholb*eps_rhol_iss;
-eps_tauc=eps_tauc_iss;
-eps_tauk=eps_tauk_iss;
-eps_tauw=eps_tauw_iss;
-eps_tauf=eps_tauf_iss;
-tauw=tauwb*eps_tauw;
-tauc=taucb*eps_tauc;
-tauf=taufb*eps_tauf;
-tauk=taukb*eps_tauk;
-eps_theta=eps_theta_iss;
-eps_gh=eps_gh_iss;
-eps_TFP=eps_TFP_iss;
-eps_g=eps_g_iss;
-g=gb*eps_g_iss;
-TFP=TFPb*eps_TFP_iss;
-gh=ghb*eps_gh_iss;
-theta=thetab*eps_theta_iss;
-
-rrb=rrbb;
-tau1=tau1b;
-om1=om1b;
-om2=om2b;
-om2s=om2sb;
-
-eps_Ds=1;
-eps_phijs=1;
-Ds=Dsb*eps_Ds;
-phijs=phijsb*eps_phijs;
-
-eps_PensCorr_F=eps_PensCorr_F_iss;
-eps_PensCorr_L=eps_PensCorr_L_iss;
-PensCorr_F=eps_PensCorr_F_iss;
-PensCorr_L=eps_PensCorr_L_iss;
-
-P00	        =	P00_iss	;
-P00_foP00	=	P00_foP00_iss	;
-
-DepRatio_n=
-@#for i in rg
-+P@{i}
-@#endfor
-;
-DepRatio_d=
-@#for i in wg
-+P@{i}
-@#endfor
-;
-
-DepRatio=DepRatio_n/DepRatio_d;
-
-ZARA=age_early+length_period*(
-@#for i in erg
-+1-i@{i}
-@#endfor  
-);
-
-Ptot=Ptot_iss;
-Ptot_f=Ptot_f_iss;
-sleep=sleep_iss;
-du=du_iss;
-de=de_iss;
-dl=dl_iss;
-
-inA=inA_iss;%(tauw+tauf)*(
-%@#for i in wg
-%+n@{i}*w@{i}*P@{i}+n@{i}_f*w@{i}_f*P@{i}_f
-%@#endfor  
-%);
-
-inB=inB_iss;%tauk*rr*(
-%@#for i in nbg
-%+1/beta@{i}^ann*s@{i-1}*P@{i}
-%@#endfor  
-%)/(1+gh);
-
-in=in_iss;%tauc*ct+inA+inB+sleep;
-NBR=NBR_iss;%g*phii*gdp+(du+de+dl)-(in);
-NBRY=NBR/(phii*gdp);
-NBR2=NBR;
-tauf2=tauf;
-tauw2=tauw;
-tauc2=tauc;
-
-end;
-
-%========================================================
-% compute initial steady state and check eigenvalues 
-%========================================================
-
-resid;
-steady(solve_algo=3);
-check;
-
-%========================================================
-endval;
-%========================================================
-
-@#for i in wg
-n@{i}=n@{i}_fss;
-n@{i}_f=n@{i}_f_fss;
-u@{i}=u@{i}_fss;
-u@{i}_f=u@{i}_f_fss;
-Omega@{i}=Omega@{i}_fss;
-Omega@{i}_f=Omega@{i}_f_fss;
-w@{i}=w@{i}_fss;
-w@{i}_f=w@{i}_f_fss;
-dWHN@{i}=dWHN@{i}_fss;
-dWFN@{i}=dWFN@{i}_fss;
-dWFN@{i}_f=dWFN@{i}_f_fss;
-eps_eta@{i}=eps_eta@{i}_fss;
-eta@{i}=eta@{i}b*eps_eta@{i};
-@#endfor
-
-@#for i in erg
-i@{i}=i@{i}_fss;
-lambda@{i}=lambda@{i}_fss;
-i@{i}_f=i@{i}_f_fss;
-lambda@{i}_f=lambda@{i}_f_fss;
-WE@{i}=0;
-eps_De_@{i}=eps_De_@{i}_fss;
-De_@{i}=De_@{i}b*eps_De_@{i}_fss;
-@#endfor
-
-@#for i in fwg
-i@{i}=0;
-lambda@{i}=0;
-i@{i}_f=0;
-lambda@{i}_f=0;
-@#endfor
-
-@#for i in ag
-@#if i in endg
-s@{i}=0;
-@#else
-s@{i}=s@{i}_fss;
-@#endif
-@#endfor
-
-@#for i in ag
-beta@{i}=beta@{i}_fss;
-beta@{i}_f=beta@{i}_f_fss;
-PD@{i}=PD@{i}_fss;
-c@{i}=c@{i}_fss;
-P@{i}=P@{i}_fss;
-P@{i}_f=P@{i}_f_fss;
-@#endfor
-
-wb          =   wb_fss;
-wb_f        =   wb_f_fss;
-Omega       =   Omega_fss;
-Omega_f     =   Omega_f_fss;
-Omega_hf    =   Omega_hf_fss;
-V        	=	V_fss	;
-M        	=	M_fss	;
-qq       	=	qq_fss	;
-p        	=	p_fss	;
-N        	=	N_fss	;
-N_f      	=	N_f_fss	;
-Q        	=	Q_fss	;
-RR       	=	RR_fss	;
-H        	=	H_fss	;
-K        	=	K_fss	;
-Y        	=	Y_fss	;
-gdp      	=	gdp_fss	;
-nx       	=	nx_fss	;
-FH       	=	FH_fss	;
-pi       	=	pi_fss	;
-ct       	=	ct_fss	;
-st       	=	st_fss	;
-wshare   	=	wshare_fss	;
-rr       	=	rr_fss	;
-
-gamma    	=	gamma_fss	;
-mc       	=	mc_fss	;
-phii     	=	phii_fss	;
-D        	=	D_fss	;
-DH       	=	DH_fss	;
-DF       	=	DF_fss	;
-X        	=	X_fss	;
-bs       	=	bs_fss	;
-bsY      	=	bsY_fss	;
-P00_f       =   P00_f_fss;
-
-eps_rhol=eps_rhol_fss;
-eps_rhoe=eps_rhoe_fss;
-eps_rhou=eps_rhou_fss;
-rhou=rhoub*eps_rhou_fss;
-rhoe=rhoeb*eps_rhoe_fss;
-rhol=rholb*eps_rhol_fss;
-eps_tauc=eps_tauc_fss;
-eps_tauk=eps_tauk_fss;
-eps_tauw=eps_tauw_fss;
-eps_tauf=eps_tauf_fss;
-tauw=tauwb*eps_tauw;
-tauc=taucb*eps_tauc;
-tauf=taufb*eps_tauf;
-tauk=taukb*eps_tauk;
-eps_theta=eps_theta_fss;
-eps_gh=eps_gh_fss;
-eps_TFP=eps_TFP_fss;
-eps_g=eps_g_fss;
-g=gb*eps_g_fss;
-TFP=TFPb*eps_TFP_fss;
-gh=ghb*eps_gh_fss;
-theta=thetab*eps_theta_fss;
-
-rrb=rrbb;
-tau1=tau1b;
-om1=om1b;
-om2=om2b;
-om2s=om2sb;
-
-eps_Ds=1;
-eps_phijs=1;
-Ds=Dsb*eps_Ds;
-phijs=phijsb*eps_phijs;
-
-eps_PensCorr_F=eps_PensCorr_F_fss;
-eps_PensCorr_L=eps_PensCorr_L_fss;
-PensCorr_F=eps_PensCorr_F_fss;
-PensCorr_L=eps_PensCorr_L_fss;
-
-P00	        =	P00_fss	;
-P00_foP00	=	P00_foP00_fss	;
-
-DepRatio_n=
-@#for i in rg
-+P@{i}
-@#endfor
-;
-DepRatio_d=
-@#for i in wg
-+P@{i}
-@#endfor
-;
-
-DepRatio=DepRatio_n/DepRatio_d;
-
-ZARA=age_early+length_period*(
-@#for i in erg
-+1-i@{i}
-@#endfor  
-);
-
-Ptot=
-@#for i in ag
-+P@{i}
-@#endfor
-;
-Ptot_f=
-@#for i in ag
-+P@{i}_f
-@#endfor
-;
-
-sleep=(1+rr)*(
-@#for i in nbg
-+1/beta@{i}*(1-1/beta@{i}^(ann-1))*s@{i-1}*P@{i}
-@#endfor  
-)/(1+gh);
-
-du=rhou*(
-@#for i in wg
-+w@{i}*u@{i}*P@{i}
-@#endfor  
-);
-
-de=rhoe*(
-@#for i in erg
-+w@{i}*i@{i}*P@{i}+w@{i}_f*i@{i}_f*P@{i}_f
-@#endfor  
-);
-
-dl=
-@#for i in rg
-+rhol*wb*PensCorr_L*P@{i}+rhol*wb_f*(N_f/(N+N_f))*PensCorr_F*P@{i}_f
-@#endfor  
-;
-
-inA=(tauw+tauf)*(
-@#for i in wg
-+n@{i}*w@{i}*P@{i}+n@{i}_f*w@{i}_f*P@{i}_f
-@#endfor  
-);
-
-inB=tauk*rr*(
-@#for i in nbg
-+1/beta@{i}^ann*s@{i-1}*P@{i}
-@#endfor  
-)/(1+gh);
-
-in=tauc*ct+inA+inB+sleep;
-NBR=g*phii*gdp+(du+de+dl)-(in);
-NBRY=NBR/(phii*gdp);
-NBR2=NBR;
-tauf2=tauf;
-tauw2=tauw;
-tauc2=tauc;
-
-end;
-
-%========================================================
-% compute final steady state and check eigenvalues 
-%========================================================
-
-resid;
-steady(solve_algo=3);
-check;
-
-
-% ===================================================
-shocks;     
-% ===================================================
-
-var P00;																																																			
-   periods 1:99;																																																			
-   values (se_P00);																																																			
- 
-@#for i in nbg
-var beta@{i};																																																			
-   periods  1:99;																																																			
-   values (se_beta@{i});
-var beta@{i}_f;																																																			
-   periods  1:99;																																																			
-   values (se_beta@{i}_f);
-var PD@{i};																																																			
-   periods  1:99;																																																			
-   values (se_PD@{i});
-@#endfor  
-																																																			
-var P00_foP00;																																																			
-   periods 1:99;																																																			
-   values (se_P00_foP00);	 																																																
-
-var eps_g;																																																			
-   periods 1:99;																																																			
-   values (se_eps_g);	
-
-var eps_PensCorr_F;																																																			
-   periods 1:99;																																																			
-   values (se_eps_PensCorr_F);	
-
-var eps_PensCorr_L;																																																			
-   periods 1:99;																																																			
-   values (se_eps_PensCorr_L);	
-
-end;
-
-% *******************************************
-% Numerical Simulation, Control Parameters
-% *******************************************
-
-perfect_foresight_setup(periods=125);
-perfect_foresight_solver(maxit=100);
-
-if ~oo_.deterministic_simulation.status
-   error('Perfect foresight simulation failed')
-end
-
-mfs0=load('lola_solve_one_boundary_results');
-
-if max(max(oo_.endo_simul-mfs0.oo_.endo_simul)) > options_.dynatol.x
-   error('Inconsistency with mfs=0')
-end
diff --git a/tests/deterministic_simulations/lola_solve_one_boundary_mfs3.mod b/tests/deterministic_simulations/lola_solve_one_boundary_mfs3.mod
deleted file mode 100644
index 2a64bd222a446c02abf07026d71939cf35e5feaa..0000000000000000000000000000000000000000
--- a/tests/deterministic_simulations/lola_solve_one_boundary_mfs3.mod
+++ /dev/null
@@ -1,985 +0,0 @@
-// Tests option mfs=3 with block
-
-load lola_data.mat
-
-% ====================================================
-% declarations var -- varexo -- para   
-% ====================================================
-
-@#define nbr_work_generations=9
-@#define nbr_early_generations=2
-@#define nbr_generations=16
-
-parameters
-length_period age_early;
-
-length_period=5;
-age_early=55;
-
-@#define wt=[1]
-
-@#define wg=0:nbr_work_generations-1
-@#define ag=0:nbr_generations-1
-@#define fwg=0:nbr_work_generations-nbr_early_generations-1
-@#define nbwg=1:nbr_work_generations-1
-@#define nbg=1:nbr_generations-1
-@#define rg=nbr_work_generations:nbr_generations-1
-@#define erg=nbr_work_generations-nbr_early_generations:nbr_work_generations-1
-@#define endg=[nbr_generations-1]
-@#define endw=[nbr_work_generations-1]
-
-@#for i in wg
-var 
-n@{i} u@{i} Omega@{i} w@{i} dWHN@{i} dWFN@{i} 
-n@{i}_f u@{i}_f Omega@{i}_f w@{i}_f dWFN@{i}_f
-i@{i} lambda@{i} i@{i}_f lambda@{i}_f eta@{i};
-parameters
-Du@{i} Dn@{i} h@{i} h@{i}_f chi@{i} eta@{i}b; 
-varexo
-eps_eta@{i};
-@#endfor
-
-@#for i in ag
-var 
-c@{i} s@{i} P@{i} P@{i}_f;
-varexo
-beta@{i} beta@{i}_f PD@{i};
-@#endfor
-
-@#for i in erg
-var 
-WE@{i} De_@{i};
-parameters
-De_@{i}b; 
-varexo
-eps_De_@{i};
-@#endfor
-
-var   
-wb wb_f 
-Omega Omega_f Omega_hf 
-V M qq p 
-N N_f 
-Q RR H K Y gdp nx FH pi
-ct st wshare rr 
-gamma mc phii D DH DF X bs bsY P00_f
-
-rhou rhoe rhol tauw tauc tauf tauk g 
-TFP gh rrb
-theta tau1 om1 om2 om2s Ds phijs
-
-DepRatio DepRatio_n DepRatio_d ZARA Ptot Ptot_f sleep du de dl inA inB in 
-NBR NBRY NBR2 tauw2 tauf2 tauc2
-PensCorr_L PensCorr_F;
-
-parameters 
-rho phi delta alpha beta ann 
-fc nu aa 
-
-rhoub rhoeb rholb tauwb taucb taufb taukb gb  
-TFPb ghb rrbb
-thetab tau1b om1b om2b om2sb Dsb phijsb
-
-NBRYb bsY_iss;
-
-varexo
-P00 P00_foP00
-
-eps_rhol eps_tauw eps_tauf eps_tauc eps_tauk
-eps_rhoe eps_rhou eps_TFP eps_gh eps_theta eps_g
-eps_Ds eps_phijs eps_PensCorr_L eps_PensCorr_F;
-
-
-% ============================================================
-% initialization
-% ============================================================
-
-@#for i in wg
-set_param_value('Du@{i}',Du@{i});
-set_param_value('Dn@{i}',Dn@{i});
-set_param_value('h@{i}',h@{i});
-set_param_value('h@{i}_f',h@{i}_f);
-set_param_value('chi@{i}',chi@{i});
-set_param_value('eta@{i}b',eta@{i}b);
-@#endfor
-
-@#for i in erg
-set_param_value('De_@{i}b',De_@{i}b);
-@#endfor
-
-set_param_value('rho',rho);
-set_param_value('phi',phi);
-set_param_value('delta',delta);
-set_param_value('alpha',alpha);
-set_param_value('beta',beta);
-set_param_value('ann',ann);
-set_param_value('fc',fc);
-set_param_value('nu',nu);
-set_param_value('aa',aa);
-
-set_param_value('rhoub',rhoub);
-set_param_value('rhoeb',rhoeb);
-set_param_value('rholb',rholb);
-set_param_value('tauwb',tauwb);
-set_param_value('taucb',taucb);
-set_param_value('taufb',taufb);
-set_param_value('taukb',taukb);
-set_param_value('gb',gb);
-
-set_param_value('TFPb',TFPb);
-set_param_value('ghb',ghb);
-set_param_value('rrbb',rrbb);
-
-set_param_value('thetab',thetab);
-set_param_value('tau1b',tau1b);
-set_param_value('om1b',om1b);
-set_param_value('om2b',om2b);
-set_param_value('om2sb',om2sb);
-set_param_value('Dsb',Dsb);
-set_param_value('phijsb',phijsb);
-
-set_param_value('bsY_iss',bsY_iss);
-
-NBRYb=NBR_iss/(phii_iss*gdp_iss);
-
-
-% =======================================================
-model(block, mfs=3);
-% ======================================================
-
-%  Labor Market Variables in the home country
-%  ------------------------------------------ 
-
-@#for i in fwg
-0=lambda@{i};
-@#endfor
-
-@#for i in wg
-1=n@{i}+u@{i}+i@{i};
-@#endfor   
-                            
-i0=lambda0;
-@#for i in nbwg
-i@{i}=lambda@{i-1}(-1)+lambda@{i}*(1-lambda@{i-1}(-1));
-@#endfor 
-
-P0=beta0*P00+PD0;
-@#for i in nbg
-P@{i}=beta@{i}*P@{i-1}(-1)+PD@{i};
-@#endfor
-
-Omega0=P0;
-@#for i in nbwg
-Omega@{i}=(1-lambda@{i})*( 1-lambda@{i-1}(-1)-(1-chi@{i})*n@{i-1}(-1))*P@{i};
-@#endfor 
-
-n0=p;
-@#for i in nbwg
-n@{i}=(1-lambda@{i})*((1-p)*(1-chi@{i})*n@{i-1}(-1)+p*(1-lambda@{i-1}(-1)));
-@#endfor 
-
-N=
-@#for i in wg
-+n@{i}*P@{i}
-@#endfor
-;
-
-%  Labor Market Variables in the foreign country
-%  --------------------------------------------- 
-
-@#for i in wg
-1=n@{i}_f+u@{i}_f+i@{i}_f;
-@#endfor    
-                            
-i0_f=lambda0_f;
-@#for i in nbwg
-i@{i}_f=lambda@{i-1}_f(-1)+lambda@{i}_f*(1-lambda@{i-1}_f(-1));
-@#endfor
-
-% -----------  reproduction cross-border   --------------------
-
-P00_f=P00_foP00*P00;
-
-P0_f=beta0_f*P00_f;
-@#for i in nbg
-P@{i}_f=beta@{i}_f*P@{i-1}_f(-1);
-@#endfor 
-
-Omega0_f=P0_f;
-@#for i in nbwg
-Omega@{i}_f=(1-lambda@{i}_f)*(1-lambda@{i-1}_f(-1)-(1-chi@{i})*n@{i-1}_f(-1))*P@{i}_f;
-@#endfor 
-
-n0_f=p;
-@#for i in nbwg
-n@{i}_f=(1-lambda@{i}_f)*((1-p)*(1-chi@{i})*n@{i-1}_f(-1)+p*(1-lambda@{i-1}_f(-1)));
-@#endfor
-
-N_f=
-@#for i in wg
-+n@{i}_f*P@{i}_f
-@#endfor
-;
-
-%  Matching
-% ----------
-
-Omega=
-@#for i in wg
-+Omega@{i}
-@#endfor
-;
-
-Omega_f=
-@#for i in wg
-+Omega@{i}_f
-@#endfor
-;
-
-Omega_hf=Omega+Omega_f;
- 
-M=V*Omega_hf/(V^nu+Omega_hf^nu)^(1/nu);
-
-qq=M/V;
-p=M/Omega_hf;
-
-%  Flow Budget Constraints (no bequests)
-% --------------------------------------
-
-rhou*w0*u0+ (1-tauw)*w0*n0 = (1+tauc)*c0+s0;
-@#for i in nbwg
-(1+rr*(1-tauk))/(beta@{i})^ann*s@{i-1}(-1)/(1+gh)+rhou*w@{i}*u@{i}+rhoe*w@{i}*i@{i}+(1-tauw)*w@{i}*n@{i}=(1+tauc)*c@{i}+s@{i};
-@#endfor
-@#for i in rg
-(1+rr*(1-tauk))/(beta@{i})^ann*s@{i-1}(-1)/(1+gh)+rhol*wb=(1+tauc)*c@{i}+s@{i};
-@#endfor
-@#for i in endg
-s@{i}=0;
-@#endfor
-
-wb=
-@#for i in wg
-+w@{i}/@{nbr_work_generations}
-@#endfor
-;
-
-%  Euler Conditions
-% ------------------
-
-@#for i in nbg
-1/(1+tauc)/c@{i-1}=beta*(1+rr(+1)*(1-tauk(+1)))/(1+tauc(+1))/c@{i}(+1)*(beta@{i})^(1-ann)/(1+gh);
-@#endfor
-
-
-%  Optimal Participation Rates (Early Retirement)
-%  ----------------------------------------------
-
-@#for i in erg
-WE@{i} = 0;
-@#endfor 
-
-@#for i in erg
-@#if i in endw
-WE@{i} = ((De_@{i}b*i@{i}^(phi-1))+Du@{i}+(rhoe-rhou)*w@{i}/(1+tauc)/c@{i})*(1-i@{i})-((1-tauw-rhou)*w@{i}/(1+tauc)/c@{i}-(Dn@{i}-Du@{i}))*n@{i};
-@#else  
-WE@{i} = ((De_@{i}b*i@{i}^(phi-1))+Du@{i}+(rhoe-rhou)*w@{i}/(1+tauc)/c@{i})*(1-i@{i})-((1-tauw-rhou)*w@{i}/(1+tauc)/c@{i}-(Dn@{i}-Du@{i}))*n@{i}+ beta*beta@{i+1}*WE@{i+1};
-@#endif
-@#endfor  
-
-%  Household Surplus
-% -------------------
-
-@#for i in wg
-@#if i in endw
-dWHN@{i} = w@{i}*(1-tauw-rhou)/(1+tauc)-(Dn@{i}-Du@{i})*c@{i};
-@#else  
-dWHN@{i} = w@{i}*(1-tauw-rhou)/(1+tauc)-(Dn@{i}-Du@{i})*c@{i} + beta*beta@{i+1}*c@{i}/c@{i+1}(+1)*dWHN@{i+1}(+1)*(1-p(+1))*(1-chi@{i+1})*(1-lambda@{i+1}(+1));
-@#endif
-@#endfor    
-
-%  Foreign household
-%  ------------------------
-%  participation and wages
-% ........................
-
-@#for i in wg
-lambda@{i}=lambda@{i}_f;
-w@{i}_f=w@{i}; 
-@#endfor
-
-wb_f=
-@#for i in wg
-+w@{i}_f/@{nbr_work_generations}
-@#endfor
-;
-
-%  Firm's Behavior
-% -------------------
-
-H=
-@#for i in wg
-+h@{i}*n@{i}*P@{i}+h@{i}_f*n@{i}_f*P@{i}_f
-@#endfor
-;
-
-wshare=(1+tauf)*(
-@#for i in wg
-+w@{i}*n@{i}*P@{i}+w@{i}_f*n@{i}_f*P@{i}_f
-@#endfor
-)/(phii*gdp);
-
-Y=TFP*H^(1-alpha)*(K)^alpha;
-gdp= TFP*H^(1-alpha)*(K)^alpha-aa*V/phii-fc/phii;
-pi = phii*gdp - wshare*phii*gdp - (rr+delta)*K(-1)/(1+gh);
-(rr(+1)+delta)/(1+rr(+1)*(1-tauk(+1))) = mc*TFP*alpha*(H/(K))^(1-alpha);
-FH=TFP*(1-alpha)*((K)/H)^alpha;
-
-RR=1+rr*(1-tauk);
-rr=rrb+tau1*(exp(bsY_iss-bsY)-1);
-
-%  Firm's Surplus
-%  ---------------------
-
-@#for i in wg
-@#if i in endw
-dWFN@{i} = h@{i}*mc*FH-(1+tauf)*w@{i};
-dWFN@{i}_f = h@{i}_f*mc*FH-(1+tauf)*w@{i}_f;
-@#else  
-dWFN@{i} = h@{i}*mc*FH-(1+tauf)*w@{i} + beta@{i+1}/RR(+1)*dWFN@{i+1}(+1)*(1-chi@{i+1})*(1-lambda@{i+1}(+1))*(1+gh);
-dWFN@{i}_f = h@{i}_f*mc*FH-(1+tauf)*w@{i}_f + beta@{i+1}_f/RR(+1)*dWFN@{i+1}_f(+1)*(1-chi@{i+1})*(1-lambda@{i+1}_f(+1))*(1+gh);
-@#endif
-@#endfor
-
-%  Free Entry Condition
-%  ---------------------
-
-aa=qq/Omega_hf*(
-@#for i in wg
-+Omega@{i}*dWFN@{i}+Omega@{i}_f*dWFN@{i}_f
-@#endfor
-);
-
-%  Wage Determination (Rent Sharing)
-%  -----------------------------------
-
-@#for i in wg
-(1-eta@{i})*dWHN@{i}  =  eta@{i}*((1-tauw)/(1+tauf)/(1+tauc))*dWFN@{i};
-@#endfor
-
-%  Equilibrium Conditions
-%  ----------------------
-
-ct=
-@#for i in ag
-+c@{i}*P@{i}
-@#endfor
-;
-
-st=
-@#for i in ag
-+s@{i}*P@{i}
-@#endfor
-;
-
-%  Non-Arbitrage condition (physical capital-shares)
-% .........................
-
-Q(+1)+pi(+1)=(1+rr(+1))*Q/(1+gh);
-
-%  New Open Economy Macroeconomics (NOEM)
-%  ---------------------------------------
-
-phii=mc/theta;
-D= ct + K-(1-delta)*K(-1)/(1+gh)  + g*gdp*phii + fc+aa*V;
-DH=(1/om1*phii)^(1/(rho-1))*D;
-X=(1/om2s*phii/gamma)^(1/(rho-1))*Ds;
-DF=(1/om2*gamma*phijs)^(1/(rho-1))*D;
-nx=phii*X-phijs*gamma*DF;     
-bsY=bs/(phii*gdp);
-Y=DH+X;
-phii*gdp=ct + K-(1-delta)*K(-1)/(1+gh) + g*gdp*phii+nx;
-
-st=K+Q +bs ;
-
-%  Policies
-%  ----------
-
-rhou=rhoub*eps_rhou;
-rhoe=rhoeb*eps_rhoe;
-rhol=rholb*eps_rhol;
-g=gb*eps_g;
-
-@#for i in erg
-De_@{i}=De_@{i}b*eps_De_@{i};
-@#endfor
-
-TFP=TFPb*eps_TFP;
-gh=ghb*eps_gh;
-
-@#for i in wg
-eta@{i}=eta@{i}b*eps_eta@{i};
-@#endfor
-
-rrb=rrbb;
-
-theta=thetab*eps_theta;
-tau1=tau1b;
-om1=om1b;
-om2=om2b;
-om2s=om2sb;
-Ds=Dsb*eps_Ds;
-phijs=phijsb*eps_phijs;
-
-
-% ----------- RefDR scenario 
-
-DepRatio_n=
-@#for i in rg
-+P@{i}
-@#endfor
-;
-DepRatio_d=
-@#for i in wg
-+P@{i}
-@#endfor
-;
-
-DepRatio=DepRatio_n/DepRatio_d;
-
-ZARA=age_early+length_period*(
-@#for i in erg
-+1-i@{i}
-@#endfor  
-);
-
-% ----------- WGEM 
-
-Ptot=
-@#for i in ag
-+P@{i}
-@#endfor
-;
-Ptot_f=
-@#for i in ag
-+P@{i}_f
-@#endfor
-;
-
-sleep=(1+rr)*(
-@#for i in nbg
-+1/beta@{i}*(1-1/beta@{i}^(ann-1))*s@{i-1}(-1)*P@{i}
-@#endfor  
-)/(1+gh);
-
-du=rhou*(
-@#for i in wg
-+w@{i}*u@{i}*P@{i}
-@#endfor  
-);
-
-de=rhoe*(
-@#for i in erg
-+w@{i}*i@{i}*P@{i}+w@{i}_f*i@{i}_f*P@{i}_f
-@#endfor  
-);
-
-dl=
-@#for i in rg
-+rhol*wb*PensCorr_L*P@{i}+rhol*wb_f*(N_f/(N+N_f))*PensCorr_F*P@{i}_f
-@#endfor  
-;
-
-PensCorr_L=eps_PensCorr_L; 
-PensCorr_F=eps_PensCorr_F;   
-
-inA=(tauw+tauf)*(
-@#for i in wg
-+n@{i}*w@{i}*P@{i}+n@{i}_f*w@{i}_f*P@{i}_f
-@#endfor  
-);
-
-inB=tauk*rr*(
-@#for i in nbg
-+1/beta@{i}^ann*s@{i-1}(-1)*P@{i}
-@#endfor  
-)/(1+gh);
-
-in=tauc*ct+inA+inB+sleep;
-NBR=g*phii*gdp+(du+de+dl)-(in);
-NBRY=NBR/(phii*gdp);
-
-% ----------- WGEM Adjustment variable ---------------
-
-tauf2=tauf;
-tauw2=tauw;
-NBR2=NBR;
-tauc2=tauc;
-tauf=taufb*eps_tauf;
-%----- WGEM: adjustment through tauc
-tauc=taucb*eps_tauc;
-%----- WGEM: adjustment through tauk
-tauk=taukb*eps_tauk;
-%----- WGEM: adjustment through tauw
-tauw=tauwb*eps_tauw;
-
-end;
-
-%==================================
-initval;
-%==================================
-
-@#for i in wg
-n@{i}=n@{i}_iss;
-n@{i}_f=n@{i}_f_iss;
-u@{i}=u@{i}_iss;
-u@{i}_f=u@{i}_f_iss;
-Omega@{i}=Omega@{i}_iss;
-Omega@{i}_f=Omega@{i}_f_iss;
-w@{i}=w@{i}_iss;
-w@{i}_f=w@{i}_f_iss;
-dWHN@{i}=dWHN@{i}_iss;
-dWFN@{i}=dWFN@{i}_iss;
-dWFN@{i}_f=dWFN@{i}_f_iss;
-eps_eta@{i}=eps_eta@{i}_iss;
-eta@{i}=eta@{i}b*eps_eta@{i};
-@#endfor
-
-@#for i in erg
-i@{i}=i@{i}_iss;
-lambda@{i}=lambda@{i}_iss;
-i@{i}_f=i@{i}_f_iss;
-lambda@{i}_f=lambda@{i}_f_iss;
-WE@{i}=0;
-eps_De_@{i}=eps_De_@{i}_iss;
-De_@{i}=De_@{i}b*eps_De_@{i}_iss;
-@#endfor
-
-@#for i in fwg
-i@{i}=0;
-lambda@{i}=0;
-i@{i}_f=0;
-lambda@{i}_f=0;
-@#endfor
-
-@#for i in ag
-@#if i in endg
-s@{i}=0;
-@#else
-s@{i}=s@{i}_iss;
-@#endif
-@#endfor
-
-@#for i in ag
-beta@{i}=beta@{i}_iss;
-beta@{i}_f=beta@{i}_f_iss;
-PD@{i}=PD@{i}_iss;
-c@{i}=c@{i}_iss;
-P@{i}=P@{i}_iss;
-P@{i}_f=P@{i}_f_iss;
-@#endfor
-
-wb          =   wb_iss;
-wb_f        =   wb_f_iss;
-Omega       =   Omega_iss;
-Omega_f     =   Omega_f_iss;
-Omega_hf    =   Omega_hf_iss;
-V        	=	V_iss	;
-M        	=	M_iss	;
-qq       	=	qq_iss	;
-p        	=	p_iss	;
-N        	=	N_iss	;
-N_f      	=	N_f_iss	;
-Q        	=	Q_iss	;
-RR       	=	RR_iss	;
-H        	=	H_iss	;
-K        	=	K_iss	;
-Y        	=	Y_iss	;
-gdp      	=	gdp_iss	;
-nx       	=	nx_iss	;
-FH       	=	FH_iss	;
-pi       	=	pi_iss	;
-ct       	=	ct_iss	;
-st       	=	st_iss	;
-wshare   	=	wshare_iss	;
-rr       	=	rr_iss	;
-
-gamma    	=	gamma_iss	;
-mc       	=	mc_iss	;
-phii     	=	phii_iss	;
-D        	=	D_iss	;
-DH       	=	DH_iss	;
-DF       	=	DF_iss	;
-X        	=	X_iss	;
-bs       	=	bs_iss	;
-bsY      	=	bsY_iss	;
-P00_f       =   P00_f_iss;
-
-eps_rhol=eps_rhol_iss;
-eps_rhoe=eps_rhoe_iss;
-eps_rhou=eps_rhou_iss;
-rhou=rhoub*eps_rhou_iss;
-rhoe=rhoeb*eps_rhoe_iss;
-rhol=rholb*eps_rhol_iss;
-eps_tauc=eps_tauc_iss;
-eps_tauk=eps_tauk_iss;
-eps_tauw=eps_tauw_iss;
-eps_tauf=eps_tauf_iss;
-tauw=tauwb*eps_tauw;
-tauc=taucb*eps_tauc;
-tauf=taufb*eps_tauf;
-tauk=taukb*eps_tauk;
-eps_theta=eps_theta_iss;
-eps_gh=eps_gh_iss;
-eps_TFP=eps_TFP_iss;
-eps_g=eps_g_iss;
-g=gb*eps_g_iss;
-TFP=TFPb*eps_TFP_iss;
-gh=ghb*eps_gh_iss;
-theta=thetab*eps_theta_iss;
-
-rrb=rrbb;
-tau1=tau1b;
-om1=om1b;
-om2=om2b;
-om2s=om2sb;
-
-eps_Ds=1;
-eps_phijs=1;
-Ds=Dsb*eps_Ds;
-phijs=phijsb*eps_phijs;
-
-eps_PensCorr_F=eps_PensCorr_F_iss;
-eps_PensCorr_L=eps_PensCorr_L_iss;
-PensCorr_F=eps_PensCorr_F_iss;
-PensCorr_L=eps_PensCorr_L_iss;
-
-P00	        =	P00_iss	;
-P00_foP00	=	P00_foP00_iss	;
-
-DepRatio_n=
-@#for i in rg
-+P@{i}
-@#endfor
-;
-DepRatio_d=
-@#for i in wg
-+P@{i}
-@#endfor
-;
-
-DepRatio=DepRatio_n/DepRatio_d;
-
-ZARA=age_early+length_period*(
-@#for i in erg
-+1-i@{i}
-@#endfor  
-);
-
-Ptot=Ptot_iss;
-Ptot_f=Ptot_f_iss;
-sleep=sleep_iss;
-du=du_iss;
-de=de_iss;
-dl=dl_iss;
-
-inA=inA_iss;%(tauw+tauf)*(
-%@#for i in wg
-%+n@{i}*w@{i}*P@{i}+n@{i}_f*w@{i}_f*P@{i}_f
-%@#endfor  
-%);
-
-inB=inB_iss;%tauk*rr*(
-%@#for i in nbg
-%+1/beta@{i}^ann*s@{i-1}*P@{i}
-%@#endfor  
-%)/(1+gh);
-
-in=in_iss;%tauc*ct+inA+inB+sleep;
-NBR=NBR_iss;%g*phii*gdp+(du+de+dl)-(in);
-NBRY=NBR/(phii*gdp);
-NBR2=NBR;
-tauf2=tauf;
-tauw2=tauw;
-tauc2=tauc;
-
-end;
-
-%========================================================
-% compute initial steady state and check eigenvalues 
-%========================================================
-
-resid;
-steady(solve_algo=3);
-check;
-
-%========================================================
-endval;
-%========================================================
-
-@#for i in wg
-n@{i}=n@{i}_fss;
-n@{i}_f=n@{i}_f_fss;
-u@{i}=u@{i}_fss;
-u@{i}_f=u@{i}_f_fss;
-Omega@{i}=Omega@{i}_fss;
-Omega@{i}_f=Omega@{i}_f_fss;
-w@{i}=w@{i}_fss;
-w@{i}_f=w@{i}_f_fss;
-dWHN@{i}=dWHN@{i}_fss;
-dWFN@{i}=dWFN@{i}_fss;
-dWFN@{i}_f=dWFN@{i}_f_fss;
-eps_eta@{i}=eps_eta@{i}_fss;
-eta@{i}=eta@{i}b*eps_eta@{i};
-@#endfor
-
-@#for i in erg
-i@{i}=i@{i}_fss;
-lambda@{i}=lambda@{i}_fss;
-i@{i}_f=i@{i}_f_fss;
-lambda@{i}_f=lambda@{i}_f_fss;
-WE@{i}=0;
-eps_De_@{i}=eps_De_@{i}_fss;
-De_@{i}=De_@{i}b*eps_De_@{i}_fss;
-@#endfor
-
-@#for i in fwg
-i@{i}=0;
-lambda@{i}=0;
-i@{i}_f=0;
-lambda@{i}_f=0;
-@#endfor
-
-@#for i in ag
-@#if i in endg
-s@{i}=0;
-@#else
-s@{i}=s@{i}_fss;
-@#endif
-@#endfor
-
-@#for i in ag
-beta@{i}=beta@{i}_fss;
-beta@{i}_f=beta@{i}_f_fss;
-PD@{i}=PD@{i}_fss;
-c@{i}=c@{i}_fss;
-P@{i}=P@{i}_fss;
-P@{i}_f=P@{i}_f_fss;
-@#endfor
-
-wb          =   wb_fss;
-wb_f        =   wb_f_fss;
-Omega       =   Omega_fss;
-Omega_f     =   Omega_f_fss;
-Omega_hf    =   Omega_hf_fss;
-V        	=	V_fss	;
-M        	=	M_fss	;
-qq       	=	qq_fss	;
-p        	=	p_fss	;
-N        	=	N_fss	;
-N_f      	=	N_f_fss	;
-Q        	=	Q_fss	;
-RR       	=	RR_fss	;
-H        	=	H_fss	;
-K        	=	K_fss	;
-Y        	=	Y_fss	;
-gdp      	=	gdp_fss	;
-nx       	=	nx_fss	;
-FH       	=	FH_fss	;
-pi       	=	pi_fss	;
-ct       	=	ct_fss	;
-st       	=	st_fss	;
-wshare   	=	wshare_fss	;
-rr       	=	rr_fss	;
-
-gamma    	=	gamma_fss	;
-mc       	=	mc_fss	;
-phii     	=	phii_fss	;
-D        	=	D_fss	;
-DH       	=	DH_fss	;
-DF       	=	DF_fss	;
-X        	=	X_fss	;
-bs       	=	bs_fss	;
-bsY      	=	bsY_fss	;
-P00_f       =   P00_f_fss;
-
-eps_rhol=eps_rhol_fss;
-eps_rhoe=eps_rhoe_fss;
-eps_rhou=eps_rhou_fss;
-rhou=rhoub*eps_rhou_fss;
-rhoe=rhoeb*eps_rhoe_fss;
-rhol=rholb*eps_rhol_fss;
-eps_tauc=eps_tauc_fss;
-eps_tauk=eps_tauk_fss;
-eps_tauw=eps_tauw_fss;
-eps_tauf=eps_tauf_fss;
-tauw=tauwb*eps_tauw;
-tauc=taucb*eps_tauc;
-tauf=taufb*eps_tauf;
-tauk=taukb*eps_tauk;
-eps_theta=eps_theta_fss;
-eps_gh=eps_gh_fss;
-eps_TFP=eps_TFP_fss;
-eps_g=eps_g_fss;
-g=gb*eps_g_fss;
-TFP=TFPb*eps_TFP_fss;
-gh=ghb*eps_gh_fss;
-theta=thetab*eps_theta_fss;
-
-rrb=rrbb;
-tau1=tau1b;
-om1=om1b;
-om2=om2b;
-om2s=om2sb;
-
-eps_Ds=1;
-eps_phijs=1;
-Ds=Dsb*eps_Ds;
-phijs=phijsb*eps_phijs;
-
-eps_PensCorr_F=eps_PensCorr_F_fss;
-eps_PensCorr_L=eps_PensCorr_L_fss;
-PensCorr_F=eps_PensCorr_F_fss;
-PensCorr_L=eps_PensCorr_L_fss;
-
-P00	        =	P00_fss	;
-P00_foP00	=	P00_foP00_fss	;
-
-DepRatio_n=
-@#for i in rg
-+P@{i}
-@#endfor
-;
-DepRatio_d=
-@#for i in wg
-+P@{i}
-@#endfor
-;
-
-DepRatio=DepRatio_n/DepRatio_d;
-
-ZARA=age_early+length_period*(
-@#for i in erg
-+1-i@{i}
-@#endfor  
-);
-
-Ptot=
-@#for i in ag
-+P@{i}
-@#endfor
-;
-Ptot_f=
-@#for i in ag
-+P@{i}_f
-@#endfor
-;
-
-sleep=(1+rr)*(
-@#for i in nbg
-+1/beta@{i}*(1-1/beta@{i}^(ann-1))*s@{i-1}*P@{i}
-@#endfor  
-)/(1+gh);
-
-du=rhou*(
-@#for i in wg
-+w@{i}*u@{i}*P@{i}
-@#endfor  
-);
-
-de=rhoe*(
-@#for i in erg
-+w@{i}*i@{i}*P@{i}+w@{i}_f*i@{i}_f*P@{i}_f
-@#endfor  
-);
-
-dl=
-@#for i in rg
-+rhol*wb*PensCorr_L*P@{i}+rhol*wb_f*(N_f/(N+N_f))*PensCorr_F*P@{i}_f
-@#endfor  
-;
-
-inA=(tauw+tauf)*(
-@#for i in wg
-+n@{i}*w@{i}*P@{i}+n@{i}_f*w@{i}_f*P@{i}_f
-@#endfor  
-);
-
-inB=tauk*rr*(
-@#for i in nbg
-+1/beta@{i}^ann*s@{i-1}*P@{i}
-@#endfor  
-)/(1+gh);
-
-in=tauc*ct+inA+inB+sleep;
-NBR=g*phii*gdp+(du+de+dl)-(in);
-NBRY=NBR/(phii*gdp);
-NBR2=NBR;
-tauf2=tauf;
-tauw2=tauw;
-tauc2=tauc;
-
-end;
-
-%========================================================
-% compute final steady state and check eigenvalues 
-%========================================================
-
-resid;
-steady(solve_algo=3);
-check;
-
-
-% ===================================================
-shocks;     
-% ===================================================
-
-var P00;																																																			
-   periods 1:99;																																																			
-   values (se_P00);																																																			
- 
-@#for i in nbg
-var beta@{i};																																																			
-   periods  1:99;																																																			
-   values (se_beta@{i});
-var beta@{i}_f;																																																			
-   periods  1:99;																																																			
-   values (se_beta@{i}_f);
-var PD@{i};																																																			
-   periods  1:99;																																																			
-   values (se_PD@{i});
-@#endfor  
-																																																			
-var P00_foP00;																																																			
-   periods 1:99;																																																			
-   values (se_P00_foP00);	 																																																
-
-var eps_g;																																																			
-   periods 1:99;																																																			
-   values (se_eps_g);	
-
-var eps_PensCorr_F;																																																			
-   periods 1:99;																																																			
-   values (se_eps_PensCorr_F);	
-
-var eps_PensCorr_L;																																																			
-   periods 1:99;																																																			
-   values (se_eps_PensCorr_L);	
-
-end;
-
-% *******************************************
-% Numerical Simulation, Control Parameters
-% *******************************************
-
-perfect_foresight_setup(periods=125);
-perfect_foresight_solver(maxit=100);
-
-if ~oo_.deterministic_simulation.status
-   error('Perfect foresight simulation failed')
-end
-
-mfs0=load('lola_solve_one_boundary_results');
-
-if max(max(oo_.endo_simul-mfs0.oo_.endo_simul)) > options_.dynatol.x
-   error('Inconsistency with mfs=0')
-end
diff --git a/tests/deterministic_simulations/ramst_block_mfs1.mod b/tests/deterministic_simulations/ramst_block_mfs1.mod
index 7b9017cfa79257db6efecb457eb67ae939492d84..4044bd51604cf69916e5033aecd8bc222da0e9f7 100644
--- a/tests/deterministic_simulations/ramst_block_mfs1.mod
+++ b/tests/deterministic_simulations/ramst_block_mfs1.mod
@@ -30,8 +30,6 @@ end;
 
 steady;
 
-check;
-
 shocks;
   var x;
   periods 1;
diff --git a/tests/ecb/SURGibbs/fulton_fish.mod b/tests/ecb/SURGibbs/fulton_fish.mod
index 23e8c8f41150e7d37afe1a5af5fefac830c2b820..cfa6c9b639e58215cef2a8eb65b357f424baa136 100644
--- a/tests/ecb/SURGibbs/fulton_fish.mod
+++ b/tests/ecb/SURGibbs/fulton_fish.mod
@@ -38,6 +38,12 @@ good = [6.791587808530124
   -0.530200000000000
   -0.397400000000000];
 
-if sum(abs(M_.params-good)) > 2e-14
+if isoctave
+  tolerance = 1e-2;
+else
+  tolerance = 2e-14;
+end
+
+if sum(abs(M_.params-good)) > tolerance
     error(['sum of M_.params - good was: ' num2str(sum(abs(M_.params-good)))]);
 end
diff --git a/tests/estimation/fs2000_init_from_previous.mod b/tests/estimation/fs2000_init_from_previous.mod
new file mode 100644
index 0000000000000000000000000000000000000000..f53948cddf83dbf868c5b6c6dd58f93bfc7d2852
--- /dev/null
+++ b/tests/estimation/fs2000_init_from_previous.mod
@@ -0,0 +1,92 @@
+// Test the mh_initialize_from_previous_mcmc option of the estimation command
+
+var m P c e W R k d n l gy_obs gp_obs y dA;
+varexo e_a e_m e_b;
+
+parameters alp bet gam mst rho psi del;
+
+alp = 0.33;
+bet = 0.99;
+gam = 0.003;
+mst = 1.011;
+rho = 0.7;
+psi = 0.787;
+del = 0.02;
+
+model;
+dA = exp(gam+e_a);
+log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
+-P/(c(+1)*P(+1)*m)+bet*exp(e_b)*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
+W = l/n;
+-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
+R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
+1/(c*P)-bet*exp(e_b)*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
+c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
+P*c = m;
+m-1+d = l;
+e = exp(e_a);
+y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
+gy_obs = dA*y/y(-1);
+gp_obs = (P/P(-1))*m(-1)/dA;
+end;
+
+steady_state_model;
+  dA = exp(gam);
+  gst = 1/dA;
+  m = mst;
+  khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1));
+  xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1);
+  nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp );
+  n  = xist/(nust+xist);
+  P  = xist + nust;
+  k  = khst*n;
+
+  l  = psi*mst*n/( (1-psi)*(1-n) );
+  c  = mst/P;
+  d  = l - mst + 1;
+  y  = k^alp*n^(1-alp)*gst^alp;
+  R  = mst/bet;
+  W  = l/n;
+  ist  = y-c;
+  q  = 1 - d;
+
+  e = 1;
+  
+  gp_obs = m/dA;
+  gy_obs = dA;
+end;
+
+
+shocks;
+var e_a; stderr 0.014;
+var e_m; stderr 0.005;
+end;
+
+steady;
+
+check;
+
+estimated_params;
+alp, beta_pdf, 0.356, 0.02;
+bet, beta_pdf, 0.993, 0.002;
+gam, normal_pdf, 0.0085, 0.003;
+mst, normal_pdf, 1.0002, 0.007;
+rho, beta_pdf, 0.129, 0.223;
+psi, beta_pdf, 0.65, 0.05;
+del, beta_pdf, 0.01, 0.005;
+stderr e_a, inv_gamma_pdf, 0.035449, inf;
+stderr e_b, inv_gamma_pdf, 0.01, inf;
+stderr e_m, inv_gamma_pdf, 0.008862, inf;
+end;
+
+varobs gp_obs gy_obs;
+
+% Test mh_initialize_from_previous_mcmc option
+estimation(order=1, datafile=fsdat_simul, nobs=192, loglinear, mh_replic=50, mh_nblocks=1, posterior_sampling_method='slice',
+           mh_initialize_from_previous_mcmc, mh_initialize_from_previous_mcmc_directory='./fs2000');
+
+% Provide record file and prior file
+estimation(order=1, datafile=fsdat_simul, nobs=192, loglinear, mh_replic=50, mh_nblocks=1, posterior_sampling_method='slice',
+           mh_initialize_from_previous_mcmc, mh_initialize_from_previous_mcmc_directory='',
+           mh_initialize_from_previous_mcmc_record = './fs2000/metropolis/fs2000_mh_history_0.mat',
+           mh_initialize_from_previous_mcmc_prior = './fs2000/prior/definition.mat');
diff --git a/tests/estimation/no_init_estimation_check_first_obs/fs2000_init_check.mod b/tests/estimation/no_init_estimation_check_first_obs/fs2000_init_check.mod
new file mode 100644
index 0000000000000000000000000000000000000000..2324043a30363c82ec5198e11b4092cc558d9049
--- /dev/null
+++ b/tests/estimation/no_init_estimation_check_first_obs/fs2000_init_check.mod
@@ -0,0 +1,74 @@
+// Test for the no_init_estimation_check_first_obs option
+
+var m P c e W R k d n l gy_obs gp_obs y dA;
+varexo e_a e_m;
+
+parameters alp bet gam mst rho psi del;
+
+alp = 0.33;
+bet = 0.99;
+gam = 0.003;
+mst = 1.011;
+rho = 0.7;
+psi = 0.787;
+del = 0.02;
+
+model;
+dA = exp(gam+e_a);
+log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
+-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
+W = l/n;
+-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
+R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
+1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
+c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
+P*c = m;
+m-1+d = l;
+e = exp(e_a);
+y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
+gy_obs = dA*y/y(-1);
+gp_obs = (P/P(-1))*m(-1)/dA;
+end;
+
+initval;
+k = 6;
+m = mst;
+P = 2.25;
+c = 0.45;
+e = 1;
+W = 4;
+R = 1.02;
+d = 0.85;
+n = 0.19;
+l = 0.86;
+y = 0.6;
+gy_obs = exp(gam);
+gp_obs = exp(-gam);
+dA = exp(gam);
+end;
+
+shocks;
+var e_a; stderr 0.014;
+var e_m; stderr 0.005;
+end;
+
+steady;
+
+check;
+
+estimated_params;
+alp, beta_pdf, 0.356, 0.02;
+bet, beta_pdf, 0.993, 0.002;
+gam, normal_pdf, 0.0085, 0.003;
+mst, normal_pdf, 1.0002, 0.007;
+rho, beta_pdf, 0.129, 0.223;
+psi, beta_pdf, 0.65, 0.05;
+del, beta_pdf, 0.01, 0.005;
+stderr e_a, inv_gamma_pdf, 0.035449, inf;
+stderr e_m, inv_gamma_pdf, 0.008862, inf;
+end;
+
+varobs gp_obs gy_obs k;
+
+options_.solve_tolf = 1e-12;
+estimation(order=1,datafile=fsdat_mat,nobs=192,loglinear,mh_replic=0,use_univariate_filters_if_singularity_is_detected=0, smoother, consider_all_endogenous, no_init_estimation_check_first_obs);
diff --git a/tests/estimation/no_init_estimation_check_first_obs/fs2000_init_check_XFAIL.mod b/tests/estimation/no_init_estimation_check_first_obs/fs2000_init_check_XFAIL.mod
new file mode 100644
index 0000000000000000000000000000000000000000..7557812393f9eaebd899a46a9f608d48aa5ff102
--- /dev/null
+++ b/tests/estimation/no_init_estimation_check_first_obs/fs2000_init_check_XFAIL.mod
@@ -0,0 +1,74 @@
+// Test for the (absence of) the no_init_estimation_check_first_obs option
+
+var m P c e W R k d n l gy_obs gp_obs y dA;
+varexo e_a e_m;
+
+parameters alp bet gam mst rho psi del;
+
+alp = 0.33;
+bet = 0.99;
+gam = 0.003;
+mst = 1.011;
+rho = 0.7;
+psi = 0.787;
+del = 0.02;
+
+model;
+dA = exp(gam+e_a);
+log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
+-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
+W = l/n;
+-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
+R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
+1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
+c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
+P*c = m;
+m-1+d = l;
+e = exp(e_a);
+y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
+gy_obs = dA*y/y(-1);
+gp_obs = (P/P(-1))*m(-1)/dA;
+end;
+
+initval;
+k = 6;
+m = mst;
+P = 2.25;
+c = 0.45;
+e = 1;
+W = 4;
+R = 1.02;
+d = 0.85;
+n = 0.19;
+l = 0.86;
+y = 0.6;
+gy_obs = exp(gam);
+gp_obs = exp(-gam);
+dA = exp(gam);
+end;
+
+shocks;
+var e_a; stderr 0.014;
+var e_m; stderr 0.005;
+end;
+
+steady;
+
+check;
+
+estimated_params;
+alp, beta_pdf, 0.356, 0.02;
+bet, beta_pdf, 0.993, 0.002;
+gam, normal_pdf, 0.0085, 0.003;
+mst, normal_pdf, 1.0002, 0.007;
+rho, beta_pdf, 0.129, 0.223;
+psi, beta_pdf, 0.65, 0.05;
+del, beta_pdf, 0.01, 0.005;
+stderr e_a, inv_gamma_pdf, 0.035449, inf;
+stderr e_m, inv_gamma_pdf, 0.008862, inf;
+end;
+
+varobs gp_obs gy_obs k;
+
+options_.solve_tolf = 1e-12;
+estimation(order=1,datafile=fsdat_mat,nobs=192,loglinear,mh_replic=0,use_univariate_filters_if_singularity_is_detected=0, consider_all_endogenous, smoother);
diff --git a/tests/estimation/no_init_estimation_check_first_obs/fs2000_init_check_XFAIL2.mod b/tests/estimation/no_init_estimation_check_first_obs/fs2000_init_check_XFAIL2.mod
new file mode 100644
index 0000000000000000000000000000000000000000..b4df214cc37082ba59568c3c0f4bd6e39a15a410
--- /dev/null
+++ b/tests/estimation/no_init_estimation_check_first_obs/fs2000_init_check_XFAIL2.mod
@@ -0,0 +1,74 @@
+// Test for the no_init_estimation_check_first_obs option
+
+var m P c e W R k d n l gy_obs gp_obs y dA;
+varexo e_a e_m;
+
+parameters alp bet gam mst rho psi del;
+
+alp = 0.33;
+bet = 0.99;
+gam = 0.003;
+mst = 1.011;
+rho = 0.7;
+psi = 0.787;
+del = 0.02;
+
+model;
+dA = exp(gam+e_a);
+log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
+-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
+W = l/n;
+-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
+R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
+1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
+c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
+P*c = m;
+m-1+d = l;
+e = exp(e_a);
+y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
+gy_obs = dA*y/y(-1);
+gp_obs = (P/P(-1))*m(-1)/dA;
+end;
+
+initval;
+k = 6;
+m = mst;
+P = 2.25;
+c = 0.45;
+e = 1;
+W = 4;
+R = 1.02;
+d = 0.85;
+n = 0.19;
+l = 0.86;
+y = 0.6;
+gy_obs = exp(gam);
+gp_obs = exp(-gam);
+dA = exp(gam);
+end;
+
+shocks;
+var e_a; stderr 0.014;
+var e_m; stderr 0.005;
+end;
+
+steady;
+
+check;
+
+estimated_params;
+alp, beta_pdf, 0.356, 0.02;
+bet, beta_pdf, 0.993, 0.002;
+gam, normal_pdf, 0.0085, 0.003;
+mst, normal_pdf, 1.0002, 0.007;
+rho, beta_pdf, 0.129, 0.223;
+psi, beta_pdf, 0.65, 0.05;
+del, beta_pdf, 0.01, 0.005;
+stderr e_a, inv_gamma_pdf, 0.035449, inf;
+stderr e_m, inv_gamma_pdf, 0.008862, inf;
+end;
+
+varobs gp_obs gy_obs k;
+
+options_.solve_tolf = 1e-12;
+estimation(order=1,datafile=fsdat_mat_XFAIL,nobs=192,loglinear,mh_replic=0,use_univariate_filters_if_singularity_is_detected=0, smoother, consider_all_endogenous, no_init_estimation_check_first_obs);
diff --git a/tests/estimation/no_init_estimation_check_first_obs/fsdat_mat.m b/tests/estimation/no_init_estimation_check_first_obs/fsdat_mat.m
new file mode 100644
index 0000000000000000000000000000000000000000..423aa4f6862772389cb17e0deec7241d52584a63
--- /dev/null
+++ b/tests/estimation/no_init_estimation_check_first_obs/fsdat_mat.m
@@ -0,0 +1,3 @@
+fsdat_simul;
+k=nan(size(gp_obs));
+k(1) = 20;
diff --git a/tests/estimation/no_init_estimation_check_first_obs/fsdat_mat_XFAIL.m b/tests/estimation/no_init_estimation_check_first_obs/fsdat_mat_XFAIL.m
new file mode 100644
index 0000000000000000000000000000000000000000..7f3e6ffd9685d33084e68ca8eb98a9f78ac135da
--- /dev/null
+++ b/tests/estimation/no_init_estimation_check_first_obs/fsdat_mat_XFAIL.m
@@ -0,0 +1,3 @@
+fsdat_simul;
+k=nan(size(gp_obs));
+k(1:2) = 20;
diff --git a/tests/estimation/no_init_estimation_check_first_obs/fsdat_simul.m b/tests/estimation/no_init_estimation_check_first_obs/fsdat_simul.m
new file mode 100644
index 0000000000000000000000000000000000000000..d4f4a8066f17ba49faad004256693ebc1b9b01e9
--- /dev/null
+++ b/tests/estimation/no_init_estimation_check_first_obs/fsdat_simul.m
@@ -0,0 +1,828 @@
+gy_obs          =[
+      1.0030045
+     0.99990934
+      1.0172778
+     0.99464043
+      1.0253423
+      1.0150215
+     0.97772557
+     0.97832186
+      1.0159561
+      1.0085937
+      1.0102649
+      1.0007604
+      1.0112596
+      1.0163279
+      1.0173204
+      1.0103896
+      1.0006493
+     0.99447124
+      1.0196405
+      1.0089304
+     0.99650737
+      1.0139707
+     0.97865842
+      1.0192225
+     0.99139628
+      1.0141362
+      1.0196612
+     0.97483476
+     0.99686151
+     0.99594464
+      1.0000642
+      1.0172243
+      1.0025773
+     0.97199728
+      1.0217815
+      1.0219949
+     0.99490252
+      1.0190728
+      1.0111337
+      1.0003792
+     0.98969164
+       1.010438
+      1.0216309
+      1.0016671
+      1.0357588
+     0.98803787
+      1.0093457
+      1.0177035
+     0.98548204
+      1.0274294
+      1.0141377
+      1.0091174
+     0.96427632
+      1.0083272
+      1.0007882
+     0.99038262
+      1.0031336
+     0.99500213
+     0.98203716
+      0.9889452
+       1.011632
+     0.99451949
+     0.97291047
+     0.98750871
+     0.99992418
+     0.97657318
+     0.99930448
+      1.0008515
+      1.0044064
+     0.98133792
+      1.0091702
+      1.0087023
+      1.0119876
+      1.0143019
+      1.0311061
+     0.99340471
+      1.0057428
+     0.99197259
+      1.0071019
+     0.99448853
+      1.0061819
+      1.0070088
+      0.9950913
+      1.0302318
+      0.9817693
+      1.0072885
+     0.97355282
+     0.98782586
+      1.0136674
+     0.99863956
+      1.0205668
+     0.99611384
+      1.0073805
+     0.99691529
+      1.0089194
+      1.0030467
+      1.0112006
+      1.0260523
+     0.97803331
+     0.99423374
+      1.0043727
+      1.0140173
+      1.0111473
+     0.99524348
+     0.99775943
+      0.9958619
+      0.9982344
+      1.0210212
+      1.0022288
+      1.0014801
+       1.011456
+      1.0124871
+     0.99843599
+     0.99324886
+     0.99912838
+       1.003327
+      1.0072071
+      1.0115223
+       1.009266
+      1.0070554
+      1.0129916
+      1.0053413
+      1.0051638
+     0.99212952
+      1.0214422
+     0.98716707
+     0.99905788
+     0.98877357
+     0.98568476
+     0.99767393
+      1.0061791
+     0.98423439
+     0.99492949
+     0.98786999
+     0.99754239
+      1.0168619
+     0.99472384
+      1.0041658
+     0.98123181
+      1.0112882
+     0.99245422
+      1.0010255
+      1.0017799
+      1.0089968
+      1.0072824
+     0.99768475
+      1.0044726
+      1.0118678
+      1.0056385
+      1.0276965
+      1.0025122
+      1.0065161
+      1.0234338
+     0.99760167
+     0.98922272
+      1.0101918
+       1.011615
+      1.0085286
+      1.0074455
+     0.98866757
+     0.99959012
+      1.0129881
+     0.99127881
+     0.97971901
+      1.0185314
+       1.020054
+      1.0132605
+     0.98063643
+     0.99490253
+      1.0101531
+      1.0004526
+      1.0059109
+     0.98974491
+      1.0062391
+      1.0216488
+     0.99398446
+     0.97786609
+      1.0019274
+     0.99587153
+      1.0095881
+      1.0111887
+     0.99457649
+     0.97896734
+       1.000172
+      1.0142951
+      1.0034224
+      1.0037242
+      1.0016059
+       1.016556
+     0.99687023
+      1.0117844
+      1.0059212
+     0.98083159
+     0.98638851
+      1.0128713
+      1.0096232
+      1.0115891
+      1.0011213
+      1.0147105
+      1.0066344
+      1.0164429
+     0.99825038
+     0.99403411
+
+];
+
+gp_obs          =[
+      1.0079715
+      1.0074573
+      1.0153107
+      1.0152677
+      1.0011653
+     0.99950061
+      1.0328311
+      1.0192317
+       1.009827
+     0.99588916
+       1.007474
+      1.0113061
+     0.98696624
+     0.99978663
+     0.98240542
+     0.98861723
+     0.99008763
+      1.0185076
+      1.0052452
+     0.99447194
+      1.0092685
+        1.01208
+      1.0105237
+     0.98513875
+      1.0165628
+     0.99485934
+      1.0050255
+      1.0140756
+      1.0093128
+      1.0155868
+      1.0107023
+     0.99212762
+      1.0095465
+      1.0028435
+      1.0069437
+      1.0070473
+      1.0145902
+      1.0186922
+      1.0059917
+      1.0113072
+      1.0107386
+     0.99769196
+     0.99793444
+      1.0050791
+     0.98307821
+      1.0107594
+     0.99689982
+     0.98667064
+      0.9991662
+     0.98274722
+     0.98422032
+     0.99393016
+      1.0118567
+     0.99912781
+      1.0023744
+      1.0086662
+      1.0164773
+      1.0169327
+      1.0372478
+      1.0314242
+      1.0004256
+      1.0110541
+      1.0076575
+      1.0119851
+      1.0055188
+      1.0213959
+      1.0234416
+      1.0264917
+      1.0292725
+      1.0385184
+      1.0200999
+      1.0107697
+       1.008583
+      1.0200332
+      1.0030413
+      1.0108659
+      1.0185145
+      1.0168619
+      1.0180462
+      1.0239657
+      1.0205509
+      1.0189973
+      1.0246446
+      1.0135089
+      1.0352973
+      1.0099289
+      1.0266474
+      1.0279829
+      1.0101653
+       1.041216
+      1.0103861
+      1.0114727
+      1.0054605
+      1.0190722
+      1.0114837
+      1.0179213
+       1.006082
+      1.0049696
+      1.0143629
+      0.9971036
+      1.0005602
+      1.0078403
+      1.0240222
+      1.0195063
+      1.0355136
+      1.0218743
+      1.0171331
+      1.0049817
+      1.0140974
+      1.0168431
+      1.0049966
+      1.0045568
+      1.0156414
+      1.0273055
+      1.0197653
+      1.0030624
+      1.0154993
+     0.99782084
+     0.99711648
+       1.014408
+      1.0057417
+     0.99936837
+      1.0096934
+      1.0095138
+      1.0057734
+      1.0114497
+      1.0059784
+      1.0328889
+      1.0098032
+      1.0041114
+      1.0101247
+      1.0181588
+      1.0115712
+      1.0227509
+      1.0065104
+      1.0110902
+      1.0298169
+      1.0089532
+      1.0368733
+      1.0123033
+      1.0060763
+      1.0150937
+      1.0239325
+     0.99555536
+     0.99861271
+      1.0076201
+     0.99941535
+      1.0119522
+      1.0129183
+     0.99288924
+      1.0260784
+      1.0144982
+      1.0121985
+      1.0234916
+        1.02215
+      1.0190118
+      1.0172679
+      1.0118398
+      1.0002123
+      1.0092124
+      1.0071943
+     0.99508468
+      1.0019303
+      1.0030733
+      0.9964198
+      1.0027298
+     0.99797614
+       1.006942
+     0.99793928
+      1.0083214
+      1.0283732
+      1.0111102
+       1.016936
+      1.0229061
+     0.98846454
+      1.0015387
+      1.0201769
+      1.0079822
+      1.0064007
+      1.0095543
+      1.0092207
+      1.0135485
+      1.0198974
+      1.0140252
+      1.0128686
+      1.0092903
+      1.0141974
+      1.0023492
+     0.99731455
+      1.0026598
+     0.99303643
+      1.0036469
+      1.0160975
+      1.0368378
+      1.0139625
+        1.01493
+      1.0113531
+      1.0114548
+     0.99833441
+     0.99648401
+     0.97645361
+      1.0154053
+        1.01703
+
+];
+
+Y_obs           =[
+              1
+     0.99690484
+      1.0111781
+      1.0028141
+      1.0251518
+      1.0371688
+      1.0118899
+     0.98720726
+      1.0001589
+      1.0057481
+      1.0130085
+      1.0107643
+      1.0190194
+      1.0323428
+      1.0466587
+      1.0540438
+      1.0516886
+      1.0431553
+      1.0597913
+      1.0657172
+      1.0592201
+      1.0701863
+      1.0458402
+      1.0620582
+      1.0504499
+      1.0615817
+      1.0782384
+      1.0500687
+      1.0439257
+      1.0368658
+      1.0339255
+      1.0481453
+      1.0477181
+      1.0167109
+      1.0354878
+      1.0544782
+      1.0463762
+      1.0624445
+      1.0705737
+      1.0679484
+      1.0546356
+      1.0620691
+      1.0806955
+      1.0793581
+      1.1121124
+      1.0971458
+      1.1034869
+      1.1181859
+      1.1006634
+      1.1250883
+      1.1362214
+      1.1423343
+      1.1036061
+      1.1089288
+      1.1067125
+      1.0940906
+      1.0942197
+      1.0862174
+        1.06525
+      1.0511907
+      1.0598182
+      1.0513331
+      1.0212391
+      1.0057433
+       1.002663
+     0.97623167
+     0.97253165
+     0.97037865
+     0.97178055
+     0.95011397
+     0.95627969
+     0.96197747
+     0.97096053
+     0.98225794
+      1.0103595
+      1.0007597
+       1.003498
+     0.99246608
+     0.99656347
+     0.98804749
+     0.99122491
+     0.99522926
+     0.98731605
+      1.0145434
+     0.99330816
+     0.99759216
+     0.96814048
+     0.95296183
+     0.96362471
+     0.95925977
+     0.97682205
+     0.96993138
+      0.9743074
+     0.96821818
+     0.97413308
+      0.9741753
+     0.98237142
+      1.0054193
+     0.98044807
+      0.9716773
+      0.9730455
+     0.98405828
+     0.99220103
+     0.98444001
+     0.97919493
+     0.97205233
+     0.96728223
+     0.98529893
+     0.98452324
+     0.98299888
+     0.99145042
+       1.000933
+     0.99636447
+     0.98660883
+     0.98273271
+     0.98305518
+     0.98725774
+     0.99577549
+       1.002037
+      1.0060879
+       1.016075
+      1.0184118
+      1.0205711
+      1.0096961
+      1.0281337
+      1.0122963
+      1.0083497
+     0.99411874
+       0.976799
+     0.97146842
+     0.97464304
+     0.95587292
+     0.94779791
+     0.93266339
+     0.92720128
+     0.94105864
+     0.93277798
+     0.93393927
+     0.91216657
+     0.92045028
+         0.9099
+     0.90792098
+     0.90669634
+     0.91268867
+     0.91696661
+     0.91164685
+     0.91311495
+     0.92197825
+     0.92461222
+     0.94930422
+      0.9488119
+     0.95232353
+     0.97275278
+     0.96734995
+     0.95356817
+     0.96075548
+     0.96936594
+     0.97489002
+     0.97933106
+     0.96499412
+     0.96157973
+     0.97156334
+     0.95983765
+     0.93655215
+     0.95207909
+     0.96912862
+     0.97938462
+     0.95701655
+     0.94891457
+     0.95606317
+     0.95351125
+     0.95641767
+     0.94315807
+     0.94639265
+     0.96503697
+     0.95601693
+     0.93087851
+     0.92980141
+     0.92266844
+     0.92925206
+     0.93743628
+     0.92900826
+      0.9049711
+     0.90213859
+     0.91342916
+     0.91384707
+     0.91456681
+     0.91316822
+     0.92671976
+     0.92058549
+     0.92936541
+     0.93228212
+     0.91010921
+     0.89349322
+     0.90336005
+     0.90997873
+     0.91856328
+     0.91668007
+     0.92838606
+       0.932016
+     0.94545438
+     0.94070026
+     0.93172987
+
+];
+
+P_obs           =[
+              1
+     0.99948573
+      1.0068249
+      1.0141211
+      1.0073149
+     0.99884398
+      1.0237035
+      1.0349636
+       1.036819
+      1.0247366
+      1.0242391
+      1.0275737
+      1.0065684
+     0.99838346
+     0.97281734
+     0.95346302
+      0.9355791
+      0.9461152
+     0.94338882
+     0.92988921
+      0.9311862
+     0.93529467
+     0.93784681
+     0.91501401
+     0.92360522
+     0.91049302
+     0.90754698
+     0.91365103
+     0.91499228
+     0.92260749
+     0.92533824
+     0.90949431
+     0.91106924
+     0.90594116
+     0.90491334
+      0.9039891
+     0.91060772
+     0.92132842
+     0.91934854
+     0.92268418
+     0.92545127
+     0.91517169
+     0.90513459
+     0.90224212
+     0.87734878
+     0.88013667
+     0.86906494
+     0.84776403
+     0.83895869
+     0.81373437
+     0.78998314
+     0.77594176
+     0.77982695
+     0.77098321
+     0.76538611
+     0.76608075
+     0.77458654
+     0.78354767
+     0.81282389
+     0.83627649
+     0.82873051
+     0.83181309
+     0.83149903
+     0.83551261
+     0.83305985
+     0.84648418
+     0.86195421
+     0.88047436
+     0.90177533
+     0.93232215
+     0.94445051
+      0.9472487
+     0.94786015
+     0.95992178
+     0.95499149
+     0.95788581
+      0.9684288
+     0.97731917
+     0.98739379
+      1.0033879
+      1.0159673
+      1.0269931
+      1.0436661
+      1.0492034
+      1.0765292
+      1.0784865
+      1.0971624
+      1.1171737
+      1.1193675
+      1.1526119
+      1.1550265
+      1.1585277
+      1.1560166
+      1.1671172
+      1.1706294
+      1.1805791
+      1.1786896
+      1.1756876
+      1.1820789
+       1.171211
+      1.1637997
+      1.1636684
+       1.179719
+      1.1912538
+      1.2187959
+      1.2326986
+      1.2418602
+      1.2388704
+      1.2449963
+      1.2538678
+      1.2508929
+      1.2474781
+       1.255148
+       1.274482
+      1.2862757
+      1.2813665
+      1.2888943
+      1.2787436
+      1.2678886
+       1.274325
+      1.2720952
+       1.263492
+      1.2652139
+      1.2667561
+       1.264558
+      1.2680362
+      1.2660431
+      1.2909605
+      1.2927921
+       1.288932
+      1.2910852
+      1.3012725
+      1.3048721
+      1.3196515
+      1.3181903
+       1.321309
+      1.3431543
+       1.344136
+      1.3730377
+      1.3773695
+      1.3754742
+      1.3825964
+      1.3985574
+      1.3861412
+      1.3767823
+      1.3764309
+      1.3678747
+      1.3718554
+      1.3768022
+      1.3617199
+      1.3798267
+      1.3863533
+      1.3905803
+      1.4061004
+      1.4202788
+      1.4313191
+      1.4406155
+      1.4444837
+      1.4367244
+      1.4379653
+      1.4371881
+      1.4243012
+        1.41826
+      1.4133617
+        1.40181
+      1.3965683
+      1.3865729
+      1.3855433
+      1.3755111
+      1.3758609
+      1.3962625
+      1.3994012
+      1.4083656
+      1.4233002
+      1.4037932
+      1.3973604
+      1.4095657
+      1.4095764
+      1.4080055
+      1.4095882
+      1.4108374
+      1.4164143
+      1.4283402
+      1.4343939
+      1.4392909
+      1.4406097
+      1.4468355
+      1.4412132
+      1.4305562
+      1.4252445
+      1.4103094
+      1.4059847
+      1.4141106
+      1.4429769
+      1.4489679
+      1.4559263
+      1.4593079
+      1.4627911
+       1.453154
+      1.4416665
+      1.4101485
+      1.4175823
+      1.4266407
+
+];
+