diff --git a/doc/manual/source/the-model-file.rst b/doc/manual/source/the-model-file.rst
index 07101fdb6..fb077a848 100644
--- a/doc/manual/source/the-model-file.rst
+++ b/doc/manual/source/the-model-file.rst
@@ -1769,7 +1769,7 @@ in this case ``initval`` is used to specify the terminal conditions.
     in the last ``initval`` or ``endval`` block (or the steady state
     file if you provided one, see :ref:`st-st`).
 
-.. command:: initval_file (filename = FILENAME);
+.. command:: initval_file (OPTIONS...);
 
     |br| In a deterministic setup, this command is used to specify a
     path for all endogenous and exogenous variables. The length of
@@ -1786,33 +1786,534 @@ in this case ``initval`` is used to specify the terminal conditions.
           by the path for endogenous variables for the simulation
           periods (excluding initial and terminal conditions)
 
-    The command accepts three file formats:
+    In perfect foresight and stochastic contexts, ``steady`` uses the
+    first observation loaded by ``initval_file`` as guess value to
+    solve for the steady state of the model. This first observation is
+    determined by the ``first_obs`` option when it is used.
+
+    Don’t mix ``initval_file`` with ``initval`` statements. However,
+    after ``initval_file``, you can modify the historical initial
+    values with ``histval`` or ``histval_file`` statement.
+
+    There can be several ``initval_file`` statements in a model
+    file. Each statement resets ``oo_.initval_series``.
+
+    *Options*
+
+    .. option:: datafile = FILENAME
+                filename = FILENAME (deprecated)
+		
+        The name of the file containing the data. It must be included in quotes if the filename
+        contains a path or an extension. The command accepts the following file formats:
 
         * M-file (extension ``.m``): for each endogenous and exogenous
           variable, the file must contain a row or column vector of
-          the same name. Their length must be ``periods +
-          M_.maximum_lag + M_.maximum_lead``
+          the same name.
         * MAT-file (extension ``.mat``): same as for M-files.
         * Excel file (extension ``.xls`` or ``.xlsx``): for each
-          endogenous and exogenous, the file must contain a column of
-          the same name. NB: Octave only supports the ``.xlsx`` file
-          extension and must have the `io`_ package installed (easily
-          done via octave by typing ‘``pkg install -forge io``’).
+          endogenous and exogenous variable, the file must contain a
+          column of the same name. NB: Octave only supports the
+          ``.xlsx`` file extension and must have the `io`_ package
+          installed (easily done via octave by typing ‘``pkg
+          install -forge io``’). The first column may contain the date
+	  of each observation.
+        * CSV files (extension ``.csv``): for each endogenous and
+          exogenous variable, the file must contain a column of the
+          same name. The first column may contain the date of each
+          observation.
 
-    .. warning:: The extension must be omitted in the command
-                 argument. Dynare will automatically figure out the
-                 extension and select the appropriate file type. If
-                 there are several files with the same name but different
-                 extensions, then the order of precedence is as follows:
-                 first ``.m``, then ``.mat``, ``.xls`` and finally ``.xlsx``.
+    .. option:: first_obs = {INTEGER | DATE}
+	       
+        The observation number or the date (see
+	:ref:`dates-members`) of the first observation to be used in the file
 
+    .. option:: first_simulation_period = {INTEGER | DATE}
 
-.. command:: histval_file (filename = FILENAME);
+	The observation number in the file or the date (see
+	:ref:`dates <dates-members>`) at which the simulation (or the forecast) is
+	starting. This option avoids to have to compute the maximum
+	number of lags in the model.  The observation corresponding to
+	the first period of simulation doesn’t need to exist in the
+	file as the only dates necessary for initialization are before
+	that date.
+	
+    .. option:: last_obs = {INTEGER | DATE}
+	       
+        The observaton number or the date (see
+	:ref:`dates-members`) of the last observation to be used in
+	the file.
+
+    .. option:: nobs = INTEGER
+
+	The number of observations to be used in the file (starting
+	with first of ``first_obs`` observation).	 
+
+    .. option:: series = DSERIES NAME
+
+        The name of a DSERIES containing the data (see :ref:`dseries-members`)		
+
+    *Example 1*
+
+        ::
+
+            var c x;
+            varexo e;
+	    parameters a b c d;
+
+	    a = 1.5;
+	    b = -0,6;
+	    c = 0.5;
+	    d = 0.5;
+	    
+            model;
+            x = a*x(-1) + b*x(-2) + e;
+            log(c) = c*x + d*log(c(+1));
+            end;
+
+	    initval_file(datafile=mydata.csv);
+
+	    perfect_foresight_setup(periods=200);
+	    perfect_foresight_solver;
+
+	The initial and terminal values are taken from file
+	``mydata.csv`` (nothing guarantees that these vales are the
+	steady state of the model). The guess value for the
+	trajectories are also taken from the file. The file must
+	contain at least 203 observations of variables ``c``, ``x``
+	and ``e``. If there are more than 203 observations available
+	in the file, the first 203 are used by
+	``perfect_foresight_setup(periods=200)``.
+	Note that the values for the auxiliary variable corresponding
+	to ``x(-2)`` are automatically computed by ``initval_file``.
+	
+    *Example 2*
+
+        ::
+
+            var c x;
+            varexo e;
+	    parameters a b c d;
+
+	    a = 1.5;
+	    b = -0,6;
+	    c = 0.5;
+	    d = 0.5;
+	    
+            model;
+            x = a*x(-1) + b*x(-2) + e;
+            log(c) = c*x + d*log(c(+1));
+            end;
+
+	    initval_file(datafile=mydata.csv,
+	                 first_obs=10);
+
+	    perfect_foresight_setup(periods=200);
+	    perfect_foresight_solver;
+
+	The initial and terminal values are taken from file
+	``mydata.csv`` starting with the 10th observation in the
+	file. There must be at least 212 observations in the file.
+	
+    *Example 3*
+
+        ::
+
+            var c x;
+            varexo e;
+	    parameters a b c d;
+
+	    a = 1.5;
+	    b = -0,6;
+	    c = 0.5;
+	    d = 0.5;
+	    
+            model;
+            x = a*x(-1) + b*x(-2) + e;
+            log(c) = c*x + d*log(c(+1));
+            end;
+
+	    ds = dseries(mydata.csv);
+	    lds = log(ds);
+	    
+	    initval_file(series=lds,
+	                 first_obs=2010Q1);
+
+	    perfect_foresight_setup(periods=200);
+	    perfect_foresight_solver;
+
+	The initial and terminal values are taken from dseries
+	``lds``. All observations are loaded starting with the 1st quarter of
+	2010 until the end of the file. There must be data available
+	at least until 2050Q3.
+	
+    *Example 4*
+
+        ::
+
+            var c x;
+            varexo e;
+	    parameters a b c d;
+
+	    a = 1.5;
+	    b = -0,6;
+	    c = 0.5;
+	    d = 0.5;
+	    
+            model;
+            x = a*x(-1) + b*x(-2) + e;
+            log(c) = c*x + d*log(c(+1));
+            end;
+
+	    initval_file(datafile=mydata.csv,
+	                 first_simulation_period=2010Q1);
+
+	    perfect_foresight_setup(periods=200);
+	    perfect_foresight_solver;
+
+	The initial and terminal values are taken from file
+	``mydata.csv``. The observations in the file must have
+	dates. All observations are loaded from the 3rd quarter of
+	2009 until the end of the file. There must be data available
+	in the file at least until 2050Q1.
+	
+    *Example 5*
+
+        ::
+
+            var c x;
+            varexo e;
+	    parameters a b c d;
+
+	    a = 1.5;
+	    b = -0,6;
+	    c = 0.5;
+	    d = 0.5;
+	    
+            model;
+            x = a*x(-1) + b*x(-2) + e;
+            log(c) = c*x + d*log(c(+1));
+            end;
+
+	    initval_file(datafile=mydata.csv,
+	                 last_obs = 212);
+
+	    perfect_foresight_setup(periods=200);
+	    perfect_foresight_solver;
+
+	The initial and terminal values are taken from file
+	``mydata.csv``. The first 212 observations are loaded and the
+	first 203 observations will be used by
+	``perfect_foresight_setup(periods=200)``.
+	
+    *Example 6*
+
+        ::
+
+            var c x;
+            varexo e;
+	    parameters a b c d;
+
+	    a = 1.5;
+	    b = -0,6;
+	    c = 0.5;
+	    d = 0.5;
+	    
+            model;
+            x = a*x(-1) + b*x(-2) + e;
+            log(c) = c*x + d*log(c(+1));
+            end;
+
+	    initval_file(datafile=mydata.csv,
+	                 first_obs = 10,
+			 nobs = 203);
+
+	    perfect_foresight_setup(periods=200);
+	    perfect_foresight_solver;
+
+	The initial and terminal values are taken from file
+	``mydata.csv``. Observations 10 to 212 are loaded.
+	   
+    *Example 7*
+
+        ::
+
+            var c x;
+            varexo e;
+	    parameters a b c d;
+
+	    a = 1.5;
+	    b = -0,6;
+	    c = 0.5;
+	    d = 0.5;
+	    
+            model;
+            x = a*x(-1) + b*x(-2) + e;
+            log(c) = c*x + d*log(c(+1));
+            end;
+
+	    initval_file(datafile=mydata.csv,
+	                 first_obs = 10);
+
+	    steady;
+
+        The values of the 10th observation of ``mydata.csv`` are used
+	as guess value to compute the steady state. The exogenous
+	variables are set to values found in the file or zero if these
+	variables aren't present.
+
+.. command:: histval_file (OPTIONS...);
 
     |br| This command is equivalent to ``histval``, except that it
     reads its input from a file, and is typically used in conjunction
     with ``smoother2histval``.
 
+        *Options*
+
+    .. option:: datafile = FILENAME
+                filename = FILENAME (deprecated)
+
+        The name of the file containing the data. The command accepts
+	the following file formats:
+
+        * M-file (extension ``.m``): for each endogenous and exogenous
+          variable, the file must contain a row or column vector of
+          the same name.
+        * MAT-file (extension ``.mat``): same as for M-files.
+        * Excel file (extension ``.xls`` or ``.xlsx``): for each
+          endogenous and exogenous variable, the file must contain a
+          column of the same name. NB: Octave only supports the
+          ``.xlsx`` file extension and must have the `io`_ package
+          installed (easily done via octave by typing ‘``pkg
+          install -forge io``’).  The first column may contain the
+          date of each observation.
+        * CSV files (extension ``.csv``): for each endogenous and
+          exogenous variable, the file must contain a column of the
+          same name. The first column may contain the date of each
+          observation.
+
+    .. option:: first_obs = {INTEGER | DATE}
+	       
+        The observation number or the date (see :ref:`dates-members`) of
+	the first observation to be used in the file
+
+    .. option:: first_simulation_period = {INTEGER | DATE}
+
+	The observation number in the file or the date (see
+	:ref:`dates-members`) at which the simulation (or the forecast) is
+	starting. This option avoids to have to compute the maximum
+	number of lags in the model.  The observation corresponding to
+	the first period of simulation doesn’t need to exist in the
+	file as the only dates necessary for initialization are before
+	that date.
+	
+    .. option:: last_obs = {INTEGER | DATE}
+	       
+        The observation number or the date (see :ref:`dates-members`) of the
+	last observation to be used in the file.
+
+    .. option:: nobs = INTEGER
+
+	The number of observations to be used in the file (starting
+	with first of ``first_obs`` observation).	 
+
+    .. option:: series = DSERIES NAME
+
+        The name of a DSERIES containing the data (see :ref:`dseries-members`)		
+
+    *Example 1*
+
+        ::
+
+            var c x;
+            varexo e;
+	    parameters a b c d;
+
+	    a = 1.5;
+	    b = -0,6;
+	    c = 0.5;
+	    d = 0.5;
+	    
+            model;
+            x = a*x(-1) + b*x(-2) + e;
+            log(c) = c*x + d*log(c(+1));
+            end;
+
+	    steady_state_model;
+	    x = 0;
+	    c = exp(c*x/(1 - d));
+	    end;
+	    
+	    histval_file(datafile=mydata.csv);
+
+	    stoch_simul(order=1,periods=100);
+	    
+	The initial values for the stochastic simulation are taken
+	from the two first rows of file ``mydata.csv``.
+	   
+    *Example 2*
+
+        ::
+
+            var c x;
+            varexo e;
+	    parameters a b c d;
+
+	    a = 1.5;
+	    b = -0,6;
+	    c = 0.5;
+	    d = 0.5;
+	    
+            model;
+            x = a*x(-1) + b*x(-2) + e;
+            log(c) = c*x + d*log(c(+1));
+            end;
+
+	    histval_file(datafile=mydata.csv,
+	                 first_obs=10);
+
+	    stoch_simul(order=1,periods=100);
+
+	The initial values for the stochastic simulation are taken
+	from rows 10 and 11 of file ``mydata.csv``.
+
+	
+    *Example 3*
+
+        ::
+
+            var c x;
+            varexo e;
+	    parameters a b c d;
+
+	    a = 1.5;
+	    b = -0,6;
+	    c = 0.5;
+	    d = 0.5;
+	    
+            model;
+            x = a*x(-1) + b*x(-2) + e;
+            log(c) = c*x + d*log(c(+1));
+            end;
+
+	    histval_file(datafile=mydata.csv,
+	                 first_obs=2010Q1);
+
+	    stoch_simul(order=1,periods=100);
+
+	The initial values for the stochastic simulation are taken
+	from observations 2010Q1 and 2010Q2 of file ``mydata.csv``.
+	
+    *Example 4*
+
+        ::
+
+            var c x;
+            varexo e;
+	    parameters a b c d;
+
+	    a = 1.5;
+	    b = -0,6;
+	    c = 0.5;
+	    d = 0.5;
+	    
+            model;
+            x = a*x(-1) + b*x(-2) + e;
+            log(c) = c*x + d*log(c(+1));
+            end;
+
+	    histval_file(datafile=mydata.csv,
+	                 first_simulation_period=2010Q1)
+
+	    stoch_simul(order=1,periods=100);
+
+	The initial values for the stochastic simulation are taken
+	from observations 2009Q3 and 2009Q4 of file ``mydata.csv``.
+	
+    *Example 5*
+
+        ::
+
+            var c x;
+            varexo e;
+	    parameters a b c d;
+
+	    a = 1.5;
+	    b = -0,6;
+	    c = 0.5;
+	    d = 0.5;
+	    
+            model;
+            x = a*x(-1) + b*x(-2) + e;
+            log(c) = c*x + d*log(c(+1));
+            end;
+
+	    histval_file(datafile=mydata.csv,
+	                 last_obs = 4);
+
+	    stoch_simul(order=1,periods=100);
+
+	The initial values for the stochastic simulation are taken
+	from the two first rows of file ``mydata.csv``.
+	
+    *Example 6*
+
+        ::
+
+            var c x;
+            varexo e;
+	    parameters a b c d;
+
+	    a = 1.5;
+	    b = -0,6;
+	    c = 0.5;
+	    d = 0.5;
+	    
+            model;
+            x = a*x(-1) + b*x(-2) + e;
+            log(c) = c*x + d*log(c(+1));
+            end;
+
+	    initval_file(datafile=mydata.csv,
+	                 first_obs = 10,
+			 nobs = 4);
+
+	    stoch_simul(order=1,periods=100);
+
+	The initial values for the stochastic simulation are taken
+	from rows 10 and 11 of file ``mydata.csv``.
+	   
+    *Example 7*
+
+        ::
+
+            var c x;
+            varexo e;
+	    parameters a b c d;
+
+	    a = 1.5;
+	    b = -0,6;
+	    c = 0.5;
+	    d = 0.5;
+	    
+            model;
+            x = a*x(-1) + b*x(-2) + e;
+            log(c) = c*x + d*log(c(+1));
+            end;
+
+	    initval_file(datafile=mydata.csv,
+	                 first_obs=10);
+
+            histval_file(datafile=myotherdata.csv);
+	    
+	    perfect_foresight_setup(periods=200);
+	    perfect_foresight_solver;
+
+	Historical initial values for the simulation are taken from
+	the two first rows of file ``myotherdata.csv``.
+	
+        Terminal values and guess values for the simulation are taken
+	from file ``mydata.csv`` starting with the 12th observation in
+	the file. There must be at least 212 observations in the file.
 
 .. _shocks-exo:
 
diff --git a/matlab/global_initialization.m b/matlab/global_initialization.m
index 1b35a37e5..8f5f909d5 100644
--- a/matlab/global_initialization.m
+++ b/matlab/global_initialization.m
@@ -92,6 +92,7 @@ oo_.dr = [];
 oo_.exo_steady_state = [];
 oo_.exo_det_steady_state = [];
 oo_.exo_det_simul = [];
+oo_.initval_series = dseries();
 
 oo_.gui.ran_estimation = false;
 oo_.gui.ran_stoch_simul = false;
diff --git a/matlab/histvalf.m b/matlab/histvalf.m
index f1f9fb098..f15e72511 100644
--- a/matlab/histvalf.m
+++ b/matlab/histvalf.m
@@ -1,5 +1,5 @@
-function histvalf(fname)
-%function histvalf(fname)
+function [endo_histval, exo_histval, exo_det_histval] = histvalf(M, options)
+%function [endo_histval, exo_histval, exo_det_histval] = histvalf(M, options)
 % Sets initial values for simulation using values contained in `fname`, a
 % file possibly created by a call to `smoother2histval`
 %
@@ -13,7 +13,7 @@ function histvalf(fname)
 %    none
 
 
-% Copyright (C) 2014-2019 Dynare Team
+% Copyright (C) 2014-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -30,72 +30,41 @@ function histvalf(fname)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-global M_ oo_ ex0_
-
-if ~exist(fname, 'file')
-    error(['Can''t find datafile: ' fname ]);
+if ~isfield(options, 'nobs') || isempty(options.nobs)
+    options.nobs = M.orig_maximum_lag;
 end
 
-M_.endo_histval = repmat(oo_.steady_state, 1, M_.maximum_endo_lag);
-
-% Also fill in oo_.exo_simul: necessary if we are in deterministic context,
-% since aux vars for lagged exo are not created in this case
-if isempty(oo_.exo_simul)
-    if isempty(ex0_)
-        oo_.exo_simul = repmat(oo_.exo_steady_state',M_.maximum_lag,1);
+if ~isfield(options, 'first_obs') || isempty(options.first_obs)
+    if isfield(options, 'first_simulation_period')
+        options.first_obs = options.first_simulation_period ...
+            - options.nobs;
     else
-        oo_.exo_simul = repmat(ex0_',M_.maximum_lag,1);
+        options.first_obs = 1;
+    end
+elseif isfield(options, 'first_simulation_period')
+    nobs = options.first_simulation_period - opions_.first_obs;
+    if options.nobs ~= nobs
+        error(sprintf(['HISTVALF: first_obs = %d and', ...
+                       ' first_simulation_period = %d', ...
+                       ' don''t provide for the number of' ...
+                       ' lags in the model.'], ...
+                      options.first_obs, ...
+                      options.first_simulation_period))
     end
 end
 
-S = load(fname);
+series = histvalf_initvalf('HISTVAL', M, options);
+% capture the difference between stochastic and
+% perfect foresight setup
+k = M.orig_maximum_lag - M.maximum_lag + 1;
+endo_histval  = series{M.endo_names{:}}.data(k:end, :)';
 
-outvars = fieldnames(S);
-
-for i = 1:length(outvars)
-    ov_ = outvars{i};
-    if ov_(end) == '_'
-        ov = ov_(1:end-1);
-        j = strmatch(ov, M_.endo_names, 'exact');
-        if isempty(j)
-            warning(['smoother2histval: output variable ' ov ' does not exist.'])
-        end
-    else
-        % Lagged endogenous or exogenous, search through aux vars
-        undidx = find(ov_ == '_', 1, 'last'); % Index of last underscore in name
-        ov = ov_(1:(undidx-1));
-        lead_lag = ov_((undidx+1):end);
-        lead_lag = regexprep(lead_lag,'l','-');
-        lead_lag = str2num(lead_lag);
-        j = [];
-        for i = 1:length(M_.aux_vars)
-            if M_.aux_vars(i).type ~= 1 && M_.aux_vars(i).type ~= 3
-                continue
-            end
-            if M_.aux_vars(i).type == 1
-                % Endogenous
-                orig_var = M_.endo_names{M_.aux_vars(i).orig_index};
-            else
-                % Exogenous
-                orig_var = M_.exo_names{M_.aux_vars(i).orig_index};
-            end
-            if strcmp(orig_var, ov) && M_.aux_vars(i).orig_lead_lag == lead_lag
-                j = M_.aux_vars(i).endo_index;
-            end
-        end
-        if isempty(j)
-            % There is no aux var corresponding to (orig_var, lead_lag).
-            % If this is an exogenous variable, then it means we should put
-            % the value in oo_.exo_simul (we are probably in deterministic
-            % context).
-            k = strmatch(ov, M_.exo_names);
-            if isempty(k)
-                warning(['smoother2histval: output variable ' ov '(' lead_lag ') does not exist.'])
-            else
-                oo_.exo_simul((M_.maximum_lag-M_.maximum_endo_lag+1):M_.maximum_lag, k) = S.(ov_);
-            end
-            continue
-        end
-    end
-    M_.endo_histval(j, :) = S.(ov_);
+exo_histval  = [];
+if M.exo_nbr
+    exo_histval  = series{M.exo_names{:}}.data(k:end, :)';
 end
+exo_det_histval  = [];
+if M.exo_det_nbr
+    exo_det_histval  = series{M.exo_names{:}}.data(k:end, :)';
+end
+
diff --git a/matlab/histvalf_initvalf.m b/matlab/histvalf_initvalf.m
new file mode 100644
index 000000000..e65779c49
--- /dev/null
+++ b/matlab/histvalf_initvalf.m
@@ -0,0 +1,223 @@
+function series = histvalf_initvalf(caller, M, options) 
+% function initvalf(M)
+%
+% handles options for histvalf_initvalf() and initvalf()
+%
+% INPUTS
+%    caller:           string, name of calling function
+%    M:                model structure
+%    options:          options specific to initivalf
+%
+% OUTPUTS
+%    series:           dseries containing selected data from a file or a dseries
+%
+
+% Copyright (C) 2003-2020 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+
+
+% dseries
+if isfield(options, 'series')
+    series = options.series;
+    dseries_ispresent = true;
+else
+    dseries_ispresent = false;
+end
+
+% file
+datafile = '';
+if isfield(options, 'filename')
+    warning([caller, '_FILE: option FILENAME is deprecated, please use', ...
+                 ' option DATAFILE'])
+    if dseries_ispresent
+        error([caller, '_FILE: you can''t use option FILENAME and option SERIES', ...
+               ' at the same time'])
+    end
+    if isfield(options, 'datafile')
+        error([caller, '_FILE: you can''t use option DATAFILE and option FILENAME', ...
+               ' at the same time'])
+    end
+    datafile = options.filename;
+end
+
+if isfield(options, 'datafile')
+    if dseries_ispresent
+        error([caller, '_FILE: you can''t use option DATAFILE and option SERIES', ...
+               ' at the same time'])
+    end
+    datafile = options.datafile;
+end
+
+if datafile
+    [directory,basename,extension] = fileparts(datafile);
+    % Auto-detect extension if not provided
+    if isempty(extension)
+        if exist([basename '.m'],'file')
+            extension = '.m';
+        elseif exist([basename '.mat'],'file')
+            extension = '.mat';
+        elseif exist([basename '.xls'],'file')
+            extension = '.xls';
+        elseif exist([basename '.xlsx'],'file')
+            extension = '.xlsx';
+        else
+            error([caller, '_FILE: Can''t find datafile: ' basename '.{m,mat,xls,xlsx}']);
+        end
+    end
+
+    fullname = [basename extension];
+    series = dseries(fullname);
+end
+
+% checking that all variable are present
+error_flag = false;
+for i = 1:M.orig_endo_nbr
+    if ~series.exist(M.endo_names{i})
+        disp(sprintf('%s_FILE: endogenous variable %s is missing', ...
+                     caller, M.endo_names{i}))
+        error_flag = true;
+    end
+end
+
+for i = 1:M.exo_nbr
+    if ~series.exist(M.exo_names{i})
+        disp(sprintf('%s_FILE: exogenous variable %s is missing', ...
+                     caller, M.exo_names{i}))
+        error_flag = true;
+    end
+end
+
+for i = 1:M.exo_det_nbr
+    if ~series.exist(M.exo_det_names{i})
+        disp(sprintf('%s_FILE: exo_det variable %s is missing', ...
+                     caller, M.exo_det_names{i}))
+        error_flag = true;
+    end
+end
+
+if error_flag
+    error([caller, '_FILE: some variables are missing'])
+end
+
+if exist(sprintf('+%s/dynamic_set_auxiliary_series', M.fname), 'file')
+    series = feval(sprintf('%s.dynamic_set_auxiliary_series', M.fname), series, M.params);
+end
+
+% selecting observations
+if isfield(options, 'nobs')
+    nobs = options.nobs;
+else
+    nobs = 0;
+end
+
+periods = series.dates;
+nobs0 = series.nobs;
+
+first_obs_ispresent = false;
+last_obs_ispresent = false;
+if isfield(options, 'first_obs')
+    i = options.first_obs;
+    if i < 1
+        error([caller, '_FILE: the first requested period is before available', ...
+               ' data.'])
+    elseif i > nobs0
+        error([caller, '_FILE: the first requested period is after available', ...
+               ' data.'])
+    end
+    first_obs = periods(i);
+    if nobs > 0
+        last_obs = first_obs + nobs - 1;
+        last_obs_ispresent = true;
+    end
+    first_obs_ispresent = true;
+elseif isfield(options, 'firstobs')
+    first_obs = options.firstobs; 
+    if nobs > 0
+        last_obs = first_obs + nobs - 1;
+        last_obs_ispresent = true;
+    end
+    first_obs_ispresent = true;
+end
+
+if last_obs_ispresent
+    if isfield(options, 'last_obs')
+        i = options.last_obs;
+        if i < 1
+            error([caller, '_FILE: the last requested period is before available', ...
+                   ' data.'])
+        elseif i > nobs0
+            error([caller, '_FILE: the last requested period is after available', ...
+                   ' data.'])
+        end
+        if last_obs ~= periods(i) 
+            error([caller, '_FILE: FIST_OBS, LAST_OBS and NOBS contain', ...
+                   ' inconsistent information. Use only two of these', ...
+                   ' options.'])
+        end    
+    elseif isfield(options, 'lastobs')
+        if last_obs ~= options.lastobs 
+            error([caller, '_FILE: FIST_OBS, LAST_OBS and NOBS contain', ...
+                   ' inconsistent information. Use only two of these', ...
+                   ' options.'])
+        end    
+    end
+elseif isfield(options, 'last_obs')
+    i = options.last_obs;
+    if i < 1
+        error([caller, '_FILE: the last requested period is before available', ...
+               ' data.'])
+    elseif i > nobs0
+        error([caller, '_FILE: the last requested period is after available', ...
+               ' data.'])
+    end
+    last_obs = periods(i);
+    if nobs > 0
+        first_obs = last_obs - nobs + 1;
+        first_obs_ispresent = true;
+    end
+    last_obs_ispresent = true;
+elseif isfield(options, 'lastobs')
+    last_obs = options.lastobs; 
+    if nobs > 0
+        first_obs = last_obs - nobs + 1;
+        first_obs_ispresent = true;
+    end
+    last_obs_ispresent = true;
+end
+
+if ~first_obs_ispresent
+    first_obs = periods(1);
+end
+
+if ~last_obs_ispresent
+    if nobs > 0
+        last_obs = first_obs + nobs - 1;
+    else
+        last_obs = periods(end);
+    end
+end
+
+if first_obs < series.init
+    error([caller, '_FILE: the first requested period is before available', ...
+                    ' data.'])
+elseif last_obs > series.last
+    error([caller, '_FILE: the last requested period is after available', ...
+                    ' data.'])
+else
+    series = series(first_obs:last_obs);
+end
+
diff --git a/matlab/initvalf.m b/matlab/initvalf.m
index f301e045f..17a5ff54a 100644
--- a/matlab/initvalf.m
+++ b/matlab/initvalf.m
@@ -1,19 +1,18 @@
-function initvalf(fname_)
-% function initvalf(fname_)
+function series = initvalf(M, options) 
+% function initvalf(M)
 %
-% Reads an initial path from the 'fname_' file for exogenous and endogenous variables
+% handles options for histvalf() and initvalf()
 %
 % INPUTS
-%    fname_:         name of the function or file containing the data
+%    caller:           string, name of calling function
+%    M:                model structure
+%    options:          options specific to initivalf
 %
 % OUTPUTS
-%    none
+%    series:           dseries containing selected data from a file or a dseries
 %
-% SPECIAL REQUIREMENTS
-%    All variables local to this function have an underscore appended to
-%    their name, to minimize clashes with model variables loaded by this function.
 
-% Copyright (C) 2003-2018 Dynare Team
+% Copyright (C) 2003-2020 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -30,81 +29,6 @@ function initvalf(fname_)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-global M_ oo_ options_
+series = histvalf_initvalf('INITVALF', M, options);
 
-series_ = 1;
 
-[directory,basename,extension] = fileparts(fname_);
-
-% Auto-detect extension if not provided
-if isempty(extension)
-    if exist([basename '.m'],'file')
-        extension = '.m';
-    elseif exist([basename '.mat'],'file')
-        extension = '.mat';
-    elseif exist([basename '.xls'],'file')
-        extension = '.xls';
-    elseif exist([basename '.xlsx'],'file')
-        extension = '.xlsx';
-    else
-        error(['Can''t find datafile: ' basename '.{m,mat,xls,xlsx}']);
-    end
-end
-
-fullname = [basename extension];
-
-if ~exist(fullname)
-    error(['Can''t find datafile: ' fullname ]);
-end
-
-switch (extension)
-  case '.m'
-    eval(basename);
-  case '.mat'
-    load(basename);
-  case { '.xls', '.xlsx' }
-    [data_,names_v_]=xlsread(fullname); % Octave needs the extension explicitly
-    series_=0;
-  otherwise
-    error(['Unsupported extension for datafile: ' extension])
-end
-
-options_.initval_file = true;
-oo_.endo_simul = [];
-oo_.exo_simul = [];
-
-for i_=1:length(M_.endo_names)
-    if series_ == 1
-        x_ = eval(M_.endo_names{i_});
-        if size(x_,2)>size(x_,1) %oo_.endo_simul must be collection of row vectors
-            oo_.endo_simul = [oo_.endo_simul; x_];
-        else %transpose if column vector
-            oo_.endo_simul = [oo_.endo_simul; x_'];
-        end
-    else
-        k_ = strmatch(M_.endo_names{i_}, names_v_, 'exact');
-        if isempty(k_)
-            error(['INITVAL_FILE: ' M_.endo_names{i_} ' not found'])
-        end
-        x_ = data_(:,k_);
-        oo_.endo_simul = [oo_.endo_simul; x_'];
-    end
-end
-
-for i_=1:length(M_.exo_names)
-    if series_ == 1
-        x_ = eval(M_.exo_names{i_});
-        if size(x_,2)>size(x_,1) %oo_.endo_simul must be collection of row vectors
-            oo_.exo_simul = [oo_.exo_simul x_'];
-        else %if column vector
-            oo_.exo_simul = [oo_.exo_simul x_];
-        end
-    else
-        k_ = strmatch(M_.exo_names{i_}, names_v_, 'exact');
-        if isempty(k_)
-            error(['INITVAL_FILE: ' M_.exo_names{i_} ' not found'])
-        end
-        x_ = data_(:,k_);
-        oo_.exo_simul = [oo_.exo_simul x_];
-    end
-end
diff --git a/matlab/perfect-foresight-models/make_ex_.m b/matlab/perfect-foresight-models/make_ex_.m
index 4f7259048..909ba238d 100644
--- a/matlab/perfect-foresight-models/make_ex_.m
+++ b/matlab/perfect-foresight-models/make_ex_.m
@@ -41,17 +41,33 @@ if M_.exo_det_nbr > 1 && isempty(oo_.exo_det_steady_state)
 end
 
 % Initialize oo_.exo_simul
-if isempty(M_.exo_histval)
-    if isempty(ex0_)
-        oo_.exo_simul = repmat(oo_.exo_steady_state',M_.maximum_lag+options_.periods+M_.maximum_lead,1);
+if isempty(oo_.initval_series)
+    if isempty(M_.exo_histval)
+        if isempty(ex0_)
+            oo_.exo_simul = repmat(oo_.exo_steady_state',M_.maximum_lag+options_.periods+M_.maximum_lead,1);
+        else
+            oo_.exo_simul = [ repmat(ex0_',M_.maximum_lag,1) ; repmat(oo_.exo_steady_state',options_.periods+M_.maximum_lead,1) ];
+        end
     else
-        oo_.exo_simul = [ repmat(ex0_',M_.maximum_lag,1) ; repmat(oo_.exo_steady_state',options_.periods+M_.maximum_lead,1) ];
+        if isempty(ex0_)
+            oo_.exo_simul = [M_.exo_histval'; repmat(oo_.exo_steady_state',options_.periods+M_.maximum_lead,1)];
+        else
+            error('histval and endval cannot be used simultaneously')
+        end
     end
-else
-    if isempty(ex0_)
-        oo_.exo_simul = [M_.exo_histval'; repmat(oo_.exo_steady_state',options_.periods+M_.maximum_lead,1)];
-    else
-        error('histval and endval cannot be used simultaneously')
+elseif M_.exo_nbr > 0
+    x = oo_.initval_series{M_.exo_names{:}}.data;
+    oo_.exo_simul = x(1:M_.maximum_lag + options_.periods + M_.maximum_lead,:);
+    if ~isempty(M_.exo_histval)
+        oo_.exo_simul(1:M_.maximum_lag, :) ...
+            = M_.exo_histval(:, 1:M_.maximum_lag)';
+    end
+elseif M_.exo_det_nbr > 0
+    x_det = oo_.initval_series{M_.exo_det_names{:}}.data;
+    oo_.exo_det_simul = x_det(1:M_.maximum_lag + options_.periods + M_.maximum_lead,:);
+    if ~isempty(M_.exo_det_histval)
+        oo_.exo_det_simul(1:M_.maximum_lag, :) ...
+            = M_.exo_det_histval(:, 1:M_.maximum_lag)';
     end
 end
 
diff --git a/matlab/perfect-foresight-models/make_y_.m b/matlab/perfect-foresight-models/make_y_.m
index abd48a84f..d36dd746c 100644
--- a/matlab/perfect-foresight-models/make_y_.m
+++ b/matlab/perfect-foresight-models/make_y_.m
@@ -45,18 +45,31 @@ if isempty(oo_.steady_state)
     oo_.steady_state = zeros(M_.endo_nbr,1);
 end
 
-if isempty(M_.endo_histval)
-    if isempty(ys0_)
-        oo_.endo_simul = repmat(oo_.steady_state, 1, M_.maximum_lag+options_.periods+M_.maximum_lead);
+if isempty(oo_.initval_series)
+    if isempty(M_.endo_histval)
+        if isempty(ys0_)
+            oo_.endo_simul = repmat(oo_.steady_state, 1, M_.maximum_lag+options_.periods+M_.maximum_lead);
+        else
+            oo_.endo_simul = [repmat(ys0_, 1, M_.maximum_lag) repmat(oo_.steady_state, 1,options_.periods+M_.maximum_lead)];
+        end
     else
-        oo_.endo_simul = [repmat(ys0_, 1, M_.maximum_lag) repmat(oo_.steady_state, 1,options_.periods+M_.maximum_lead)];
+        if ~isempty(ys0_)
+            error('histval and endval cannot be used simultaneously')
+        end
+        % the first NaNs take care of the case where there are lags > 1 on
+        % exogenous variables
+        oo_.endo_simul = [M_.endo_histval ...
+                          repmat(oo_.steady_state, 1, options_.periods+M_.maximum_lead)];
     end
 else
-    if ~isempty(ys0_)
-        error('histval and endval cannot be used simultaneously')
+    y = oo_.initval_series{M_.endo_names{:}}.data;
+    oo_.endo_simul = y(1:M_.maximum_lag + options_.periods + ...
+                       M_.maximum_lead, :)';
+    if ~isempty(M_.endo_histval)
+        if ~isempty(ys0_)
+            error('histval and endval cannot be used simultaneously')
+        end
+        oo_.endo_simul(:,1:M_.maximum_lag) ...
+            = M_.endo_histval(:, 1:M_.maximum_lag);
     end
-    % the first NaNs take care of the case where there are lags > 1 on
-    % exogenous variables
-    oo_.endo_simul = [M_.endo_histval ...
-                      repmat(oo_.steady_state, 1, options_.periods+M_.maximum_lead)];
-end
+end
\ No newline at end of file
diff --git a/matlab/perfect-foresight-models/perfect_foresight_setup.m b/matlab/perfect-foresight-models/perfect_foresight_setup.m
index 45f6f2789..98d8e12fc 100644
--- a/matlab/perfect-foresight-models/perfect_foresight_setup.m
+++ b/matlab/perfect-foresight-models/perfect_foresight_setup.m
@@ -64,7 +64,6 @@ if ~isempty(M_.det_shocks) && options_.periods<max([M_.det_shocks.periods])
     error('PERFECT_FORESIGHT_SETUP: Please check the declaration of the shocks or increase the value of the periods option.')
 end
 
-if ~options_.initval_file
-    oo_ = make_ex_(M_,options_,oo_);
-    oo_ = make_y_(M_,options_,oo_);
-end
+oo_ = make_ex_(M_,options_,oo_);
+oo_ = make_y_(M_,options_,oo_);
+
diff --git a/matlab/perfect-foresight-models/sim1.m b/matlab/perfect-foresight-models/sim1.m
index 3bd956c91..d37afbf7e 100644
--- a/matlab/perfect-foresight-models/sim1.m
+++ b/matlab/perfect-foresight-models/sim1.m
@@ -119,7 +119,9 @@ if options.endogenous_terminal_period
 end
 
 if stop
-    if any(any(isnan(endogenousvariables))) || any(any(isinf(endogenousvariables)))
+    % initial or terminal observations may contain
+    % harmless NaN or Inf. We test only values computed above
+    if any(any(isnan(y))) || any(any(isinf(y)))
         info.status = false;% NaN or Inf occurred
         info.error = err;
         info.iterations = iter;
diff --git a/matlab/smoother2histval.m b/matlab/smoother2histval.m
index 6e1eb2c43..3ef7e438e 100644
--- a/matlab/smoother2histval.m
+++ b/matlab/smoother2histval.m
@@ -163,13 +163,15 @@ end
 % Initialize outputs
 if ~isfield(opts, 'outfile')
     % Output to M_.endo_histval
-    M_.endo_histval = repmat(oo_.steady_state, 1, M_.maximum_endo_lag);
+    M_.endo_histval = repmat(oo_.steady_state, 1, M_.maximum_lag);
 else
     % Output to a file
-    o = struct();
+    o = dseries();
 end
 
 % Handle all endogenous variables to be copied
+data = zeros(M_.orig_maximum_endo_lag, length(invars));
+k = M_.orig_maximum_endo_lag - M_.maximum_endo_lag + 1: M_.orig_maximum_lag
 for i = 1:length(invars)
     if isempty(strmatch(invars{i}, M_.endo_names))
         % Skip exogenous
@@ -177,61 +179,68 @@ for i = 1:length(invars)
     end
     s = smoothedvars.(invars{i});
     j = strmatch(invars{i}, M_.endo_names, 'exact');
-    v = s((period-M_.maximum_endo_lag+1):period);% + steady_state(j);
+    v = s((period-M_.orig_maximum_endo_lag+1):period);% + steady_state(j);
     if ~isfield(opts, 'outfile')
         j = strmatch(outvars{i}, M_.endo_names, 'exact');
         if isempty(j)
             error(['smoother2histval: output variable ' outvars{i} ' does not exist.'])
         else
-            M_.endo_histval(j, :) = v;
+            M_.endo_histval(j, :) = v(k);
         end
     else
-        % When saving to a file, x(-1) is in the variable called "x_"
-        o.([ outvars{i} '_' ]) = v;
+        data(:, i) = v';
     end
 end
+if isfield(opts, 'outfile')
+    o = dseries(data, '1Y', invars);
+end
 
-% Handle auxiliary variables for lags (both on endogenous and exogenous)
-for i = 1:length(M_.aux_vars)
-    if ~ ismember(M_.endo_names{M_.aux_vars(i).endo_index},invars)
-        if M_.aux_vars(i).type ~= 1 && M_.aux_vars(i).type ~= 3
-            continue
-        end
-        if M_.aux_vars(i).type == 1
-            % Endogenous
-            orig_var = M_.endo_names{M_.aux_vars(i).orig_index};
-        else
-            % Exogenous
-            orig_var = M_.exo_names{M_.aux_vars(i).orig_index};
-        end
-        [m, k] = ismember(orig_var, outvars);
-        if m
-            if ~isempty(strmatch(invars{k}, M_.endo_names))
-                s = smoothedvars.(invars{k});
-            else
-                s = smoothedshocks.(invars{k});
-            end
-            l = M_.aux_vars(i).orig_lead_lag;
-            if period-M_.maximum_endo_lag+1+l < 1
-                error('The period that you indicated is too small to construct initial conditions')
-            end
-            j = M_.aux_vars(i).endo_index;
-            v = s((period-M_.maximum_endo_lag+1+l):(period+l)); %+steady_state(j);
-            if ~isfield(opts, 'outfile')
-                M_.endo_histval(j, :) = v;
-            else
-                % When saving to a file, x(-2) is in the variable called "x_l2"
-                lead_lag = num2str(l);
-                lead_lag = regexprep(lead_lag, '-', 'l');
-                o.([ orig_var '_' lead_lag ]) = v;
-            end
-        end
-    end
-end
+% $$$ % Handle auxiliary variables for lags (both on endogenous and exogenous)
+% $$$ for i = 1:length(M_.aux_vars)
+% $$$     if ~ ismember(M_.endo_names{M_.aux_vars(i).endo_index},invars)
+% $$$         if M_.aux_vars(i).type ~= 1 && M_.aux_vars(i).type ~= 3
+% $$$             continue
+% $$$         end
+% $$$         if M_.aux_vars(i).type == 1
+% $$$             % Endogenous
+% $$$             orig_var = M_.endo_names{M_.aux_vars(i).orig_index};
+% $$$         else
+% $$$             % Exogenous
+% $$$             orig_var = M_.exo_names{M_.aux_vars(i).orig_index};
+% $$$         end
+% $$$         [m, k] = ismember(orig_var, outvars);
+% $$$         if m
+% $$$             if ~isempty(strmatch(invars{k}, M_.endo_names))
+% $$$                 s = smoothedvars.(invars{k});
+% $$$             else
+% $$$                 s = smoothedshocks.(invars{k});
+% $$$             end
+% $$$             l = M_.aux_vars(i).orig_lead_lag;
+% $$$             if period-M_.maximum_endo_lag+1+l < 1
+% $$$                 error('The period that you indicated is too small to construct initial conditions')
+% $$$             end
+% $$$             j = M_.aux_vars(i).endo_index;
+% $$$             v = s((period-M_.maximum_endo_lag+1+l):(period+l)); %+steady_state(j);
+% $$$             if ~isfield(opts, 'outfile')
+% $$$                 M_.endo_histval(j, :) = v;
+% $$$             else
+% $$$                 % When saving to a file, x(-2) is in the variable called "x_l2"
+% $$$                 lead_lag = num2str(l);
+% $$$                 lead_lag = regexprep(lead_lag, '-', 'l');
+% $$$                 o.([ orig_var '_' lead_lag ]) = v;
+% $$$             end
+% $$$         end
+% $$$     end
+% $$$ end
 
 % Finalize output
 if isfield(opts, 'outfile')
-    save(opts.outfile, '-struct', 'o')
+    [dir, fname, ext] = fileparts(opts.outfile);
+    if ~strcmp(ext,'.mat') && ~isempty(ext)
+        error(['smoother2hisvtval: if outfile has an extension, it must ' ...
+               'be .mat'])
+    end
+    o.save([dir fname]);
 end
 
 end
diff --git a/preprocessor b/preprocessor
index a2bea00fe..d05ffde63 160000
--- a/preprocessor
+++ b/preprocessor
@@ -1 +1 @@
-Subproject commit a2bea00fee97dd4cdc2db32690e0d492fda2edfb
+Subproject commit d05ffde63ece3fa94882be784e0fdf9aa4ec8982
diff --git a/tests/.gitignore b/tests/.gitignore
index 35ff3c844..501fbc71d 100644
--- a/tests/.gitignore
+++ b/tests/.gitignore
@@ -25,6 +25,9 @@ wsOct
 /run_test_matlab_output.txt
 
 /block_bytecode/ls2003_tmp.mod
+/histval_initval_file/data.csv
+/histval_initval_file/data.xls
+/histval_initval_file/data.xlsx
 /partial_information/PItest3aHc0PCLsimModPiYrVarobsAll_PCL*
 /partial_information/PItest3aHc0PCLsimModPiYrVarobsCNR_PCL*
 
@@ -67,10 +70,11 @@ wsOct
 !/gsa/ls2003scr_mode.mat
 !/gsa/ls2003scr_results.mat
 !/gsa/morris/nk_est_data.m
+!/histval_initval_file/histval_initval_file_unit_tests.m
+!/histval_initval_file/ramst_initval_file_data.m
 !/identification/as2007/as2007_steadystate.m
 !/identification/as2007/G_QT.mat
 !/identification/kim/kim2_steadystate.m
-!/initval_file/ramst_initval_file_data.m
 !/internals/tests.m
 !/k_order_perturbation/run_fs2000kplusplus.m
 !/kalman/likelihood/compare_kalman_routines.m
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 39caeb54d..6cda69122 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -104,9 +104,12 @@ MODFILES = \
 	discretionary_policy/dennis_1.mod \
 	discretionary_policy/dennis_1_estim.mod \
 	discretionary_policy/Gali_discretion.mod \
-	initval_file/ramst_initval_file.mod \
-	initval_file/ramst_datafile.mod \
-	ramst_normcdf_and_friends.mod \
+	histval_initval_file/ramst_initval_file.mod \
+	histval_initval_file/ramst_data.mod \
+	histval_initval_file/ramst_datafile.mod \
+	histval_initval_file/sim_exo_lead_lag.mod \
+	histval_initval_file/sim_exo_lead_lag_initvalf.mod \
+        ramst_normcdf_and_friends.mod \
 	ramst_vec.mod \
 	example1_varexo_det.mod \
 	predetermined_variables.mod \
@@ -405,7 +408,8 @@ XFAIL_MODFILES = ramst_xfail.mod \
 	estimation/tune_mh_jscale/fs2000_1_xfail.mod \
 	estimation/tune_mh_jscale/fs2000_2_xfail.mod
 
-MFILES = initval_file/ramst_initval_file_data.m
+MFILES = histval_initval_file/ramst_initval_file_data.m \
+	 histval_initval_file_unit_tests.m
 
 # Dependencies
 example1_use_dll.m.trs: example1.m.trs
@@ -503,10 +507,15 @@ deterministic_simulations/rbc_det_stack_solve_algo_7_exo_lead.o.trs: determinist
 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
 
-initval_file/ramst_initval_file.m.trs: initval_file/ramst_initval_file_data.m.tls
-initval_file/ramst_initval_file.o.trs: initval_file/ramst_initval_file_data.o.tls
-initval_file/ramst_datafile.m.trs: initval_file/ramst_initval_file_data.m.tls
-initval_file/ramst_datafile.o.trs: initval_file/ramst_initval_file_data.o.tls
+histval_initval_file/ramst_initval_file.m.trs: histval_initval_file/ramst_initval_file_data.m.tls histval_initval_file/ramst_data.m.trs
+histval_initval_file/ramst_initval_file.o.trs: histval_initval_file/ramst_initval_file_data.o.tls histval_initval_file/ramst_data.o.trs
+histval_initval_file/ramst_datafile.m.trs: histval_initval_file/ramst_initval_file_data.m.tls
+histval_initval_file/ramst_datafile.o.trs: histval_initval_file/ramst_initval_file_data.o.tls
+histval_initval_file/sim_exo_lead_lag_initvalf.m.trs: histval_initval_file/sim_exo_lead_lag.m.trs
+histval_initval_file/sim_exo_lead_lag_initvalf.o.trs: histval_initval_file/sim_exo_lead_lag.o.trs
+histval_initval_file_unit_tests.m.trs: histval_initval_file/ramst_data.m.trs
+histval_initval_file_unit_tests.o.trs: histval_initval_file/ramst_data.o.trs
+
 
 identification/rbc_ident/rbc_ident_varexo_only.m.trs: identification/rbc_ident/rbc_ident_std_as_structural_par.m.trs
 identification/rbc_ident/rbc_ident_varexo_only.o.trs: identification/rbc_ident/rbc_ident_std_as_structural_par.o.trs
@@ -782,12 +791,20 @@ o/particle: $(patsubst %.mod, %.o.trs, $(PARTICLEFILES))
 
 # Matlab TRS Files
 M_TRS_FILES = $(patsubst %.mod, %.m.trs, $(MODFILES))
-M_TRS_FILES += run_block_byte_tests_matlab.m.trs run_reporting_test_matlab.m.trs run_all_unitary_tests.m.trs
+M_TRS_FILES += 	run_block_byte_tests_matlab.m.trs \
+		run_reporting_test_matlab.m.trs \
+		run_all_unitary_tests.m.trs \
+		histval_initval_file_unit_tests.m.trs
+
 M_XFAIL_TRS_FILES = $(patsubst %.mod, %.m.trs, $(XFAIL_MODFILES))
 
 # Octave TRS Files
 O_TRS_FILES = $(patsubst %.mod, %.o.trs, $(MODFILES))
-O_TRS_FILES += run_block_byte_tests_octave.o.trs run_reporting_test_octave.o.trs run_all_unitary_tests.o.trs
+O_TRS_FILES += 	run_block_byte_tests_octave.o.trs \
+		run_reporting_test_octave.o.trs \
+		run_all_unitary_tests.o.trs \
+	 	histval_initval_file_unit_tests.o.trs
+
 O_XFAIL_TRS_FILES = $(patsubst %.mod, %.o.trs, $(XFAIL_MODFILES))
 
 # Matlab TLS Files
@@ -920,8 +937,12 @@ EXTRA_DIST = \
 	k_order_perturbation/fs2000k++.mod \
 	lmmcp/sw-common-header.inc \
 	lmmcp/sw-common-footer.inc \
-	estimation/tune_mh_jscale/fs2000.inc
-
+	estimation/tune_mh_jscale/fs2000.inc \
+	histval_initval_file_unit_tests.m \
+	histval_initval_file/my_assert.m \
+	histval_initval_file/ramst_data.xls \
+	histval_initval_file/ramst_data.xlsx \
+	histval_initval_file/ramst_initval_file_data.m
 
 if ENABLE_MATLAB
 check-local: check-matlab
@@ -1073,7 +1094,7 @@ clean-local:
 
 	rm -f estimation/test_matrix.mat
 
-	rm -f initval_file/ramst_initval_file_data_col_vec_mat.mat initval_file/ramst_initval_file_data_row_vec_mat.mat initval_file/ramst_initval_file_excel.xls
+	rm -f histval_initval_file/ramst_initval_file_data_col_vec_mat.mat histval_initval_file/ramst_initval_file_data_row_vec_mat.mat histval_initval_file/ramst_initval_file_excel.xls
 
 	rm -f loglinear/results_exp_histval.mat loglinear/results_exp.mat
 
diff --git a/tests/histval_initval_file/my_assert.m b/tests/histval_initval_file/my_assert.m
new file mode 100644
index 000000000..2ac132e12
--- /dev/null
+++ b/tests/histval_initval_file/my_assert.m
@@ -0,0 +1,4 @@
+function failed_tests = my_assert(failed_tests, success, test_name)
+if ~success
+    failed_tests = cat(1, test_failed, test_name);
+end
\ No newline at end of file
diff --git a/tests/histval_initval_file/ramst_data.mod b/tests/histval_initval_file/ramst_data.mod
new file mode 100644
index 000000000..c65b9fc81
--- /dev/null
+++ b/tests/histval_initval_file/ramst_data.mod
@@ -0,0 +1,68 @@
+/* Verify that the “datafile” option of “perfect_foresight_setup” behaves as
+   “initval_file” (see #1663) */
+
+var c k;
+varexo x;
+
+parameters alph gam delt bet aa;
+alph=0.5;
+gam=0.5;
+delt=0.02;
+bet=0.05;
+aa=0.5;
+
+
+model;
+c + k - aa*x*k(-1)^alph - (1-delt)*k(-1);
+c^(-gam) - (1+bet)^(-1)*(aa*alph*x(+1)*k^(alph-1) + 1 - delt)*c(+1)^(-gam);
+end;
+
+initval;
+x = 1;
+k = ((delt+bet)/(1.0*aa*alph))^(1/(alph-1));
+c = aa*k^alph-delt*k;
+end;
+
+steady;
+
+shocks;
+  var x;
+  periods 2;
+  values 1.1;
+end;
+
+perfect_foresight_setup(periods=200);
+perfect_foresight_solver;
+
+fh = fopen('ramst_data.m', 'w');
+fprintf(fh, 'INIT__ = ''1Y'';\n');
+fprintf(fh, 'NAMES__ = {''c'', ''k'', ''x''};\n');
+fprintf(fh, 'TEX__ = {''c'', ''k'', ''x''};\n');
+fprintf(fh, 'c = [');
+fprintf(fh, '%f ', oo_.endo_simul(1,:));
+fprintf(fh, '];\n');
+fprintf(fh, 'k = [');
+fprintf(fh, '%f ', oo_.endo_simul(2,:));
+fprintf(fh, '];\n');
+fprintf(fh, 'x = [');
+fprintf(fh, '%f ', oo_.exo_simul);
+fprintf(fh, '];\n');
+fclose(fh);
+
+INIT__ = '1Y';
+NAMES__ = {'c', 'k', 'x'};
+TEX__  = {'c', 'k', 'x'};
+eval('c = oo_.endo_simul(1,:);');
+eval('k = oo_.endo_simul(2,:);');
+eval('x = oo_.exo_simul'';');
+save('ramst_data.mat', 'INIT__', 'NAMES__', ...
+     'TEX__', 'c', 'k', 'x');
+
+fh = fopen('ramst_data.csv', 'w');
+fprintf(fh, 'c,k,x\n');
+for i = 1:size(oo_.endo_simul, 2);
+  fprintf(fh, '%f, ', oo_.endo_simul(:, i));
+  fprintf(fh, '%f\n', oo_.exo_simul(i));
+end;
+fclose(fh);
+
diff --git a/tests/histval_initval_file/ramst_data.xls b/tests/histval_initval_file/ramst_data.xls
new file mode 100644
index 000000000..050e2166c
Binary files /dev/null and b/tests/histval_initval_file/ramst_data.xls differ
diff --git a/tests/histval_initval_file/ramst_data.xlsx b/tests/histval_initval_file/ramst_data.xlsx
new file mode 100644
index 000000000..849953eac
Binary files /dev/null and b/tests/histval_initval_file/ramst_data.xlsx differ
diff --git a/tests/initval_file/ramst_datafile.mod b/tests/histval_initval_file/ramst_datafile.mod
similarity index 100%
rename from tests/initval_file/ramst_datafile.mod
rename to tests/histval_initval_file/ramst_datafile.mod
diff --git a/tests/histval_initval_file/ramst_histval_file.mod b/tests/histval_initval_file/ramst_histval_file.mod
new file mode 100644
index 000000000..8971ba04e
--- /dev/null
+++ b/tests/histval_initval_file/ramst_histval_file.mod
@@ -0,0 +1,96 @@
+/* Test for the histval_file() command. This file needs ramst_histval_data.m. */
+
+var c k;
+varexo x;
+
+parameters alph gam delt bet aa;
+alph=0.5;
+gam=0.5;
+delt=0.02;
+bet=0.05;
+aa=0.5;
+
+
+model;
+c + k - aa*x*k(-1)^alph - (1-delt)*k(-1);
+c^(-gam) - (1+bet)^(-1)*(aa*alph*x(+1)*k^(alph-1) + 1 - delt)*c(+1)^(-gam);
+end;
+
+initval;
+x = 1;
+k = ((delt+bet)/(1.0*aa*alph))^(1/(alph-1));
+c = aa*k^alph-delt*k;
+end;
+
+steady;
+
+histval_file(datafile = 'ramst_histval_data.m');
+perfect_foresight_setup(periods = 200);
+perfect_foresight_solver;
+if ~oo_.deterministic_simulation.status
+   error('Perfect foresight simulation failed')
+end
+oo_.exo_simul = [];
+oo_.endo_simul = [];
+
+initval_file(filename = ramst_initval_file_data_row_vec_mat);
+perfect_foresight_setup(periods=200);
+if oo_.exo_simul(2) ~= 1.2
+  error('initval_file problem with exogenous variable');
+end
+if oo_.endo_simul(2, 2) ~= 13
+  error('initval_file option problem with endogenous variable');
+end
+perfect_foresight_solver;
+if ~oo_.deterministic_simulation.status
+   error('Perfect foresight simulation failed')
+end
+
+initval_file(datafile = 'ramst_data.m');
+perfect_foresight_setup(periods = 200);
+perfect_foresight_solver;
+if ~oo_.deterministic_simulation.status
+   error('Perfect foresight simulation failed')
+end
+oo_.exo_simul = [];
+oo_.endo_simul = [];
+
+initval_file(datafile = 'ramst_data.mat');
+perfect_foresight_setup(periods = 200);
+perfect_foresight_solver;
+if ~oo_.deterministic_simulation.status
+   error('Perfect foresight simulation failed')
+end
+oo_.exo_simul = [];
+oo_.endo_simul = [];
+
+initval_file(datafile = 'ramst_data.csv');
+perfect_foresight_setup(periods = 200);
+perfect_foresight_solver;
+if ~oo_.deterministic_simulation.status
+   error('Perfect foresight simulation failed')
+end
+oo_.exo_simul = [];
+oo_.endo_simul = [];
+
+initval_file(datafile = 'ramst_data.xlsx');
+perfect_foresight_setup(periods = 200);
+perfect_foresight_solver;
+if ~oo_.deterministic_simulation.status
+   error('Perfect foresight simulation failed')
+end
+oo_.exo_simul = [];
+oo_.endo_simul = [];
+
+if ispc;
+  initval_file(datafile = 'ramst_data.cls');
+  perfect_foresight_setup(periods = 200);
+  perfect_foresight_solver;
+  oo_.exo_simul = [];
+  oo_.endo_simul = [];
+
+  if ~oo_.deterministic_simulation.status
+    error('Perfect foresight simulation failed')
+  end;
+end;
+
diff --git a/tests/histval_initval_file/ramst_initval_file.mod b/tests/histval_initval_file/ramst_initval_file.mod
new file mode 100644
index 000000000..8e72fe750
--- /dev/null
+++ b/tests/histval_initval_file/ramst_initval_file.mod
@@ -0,0 +1,113 @@
+/* Test for the initval_file() command. This file needs ramst_initval_file_data.m. It should give results similar to those of ramst.mod */
+
+var c k;
+varexo x;
+
+parameters alph gam delt bet aa;
+alph=0.5;
+gam=0.5;
+delt=0.02;
+bet=0.05;
+aa=0.5;
+
+
+model;
+c + k - aa*x*k(-1)^alph - (1-delt)*k(-1);
+c^(-gam) - (1+bet)^(-1)*(aa*alph*x(+1)*k^(alph-1) + 1 - delt)*c(+1)^(-gam);
+end;
+
+initval;
+x = 1;
+k = ((delt+bet)/(1.0*aa*alph))^(1/(alph-1));
+c = aa*k^alph-delt*k;
+end;
+
+steady;
+
+initval_file(filename = ramst_initval_file_data_row_vec_mat);
+perfect_foresight_setup(periods=200);
+if oo_.exo_simul(2) ~= 1.2
+  error('initval_file problem with exogenous variable');
+end
+if oo_.endo_simul(2, 2) ~= 13
+  error('initval_file option problem with endogenous variable');
+end
+perfect_foresight_solver;
+if ~oo_.deterministic_simulation.status
+   error('Perfect foresight simulation failed')
+end
+
+oo_.exo_simul = [];
+oo_.endo_simul = [];
+
+initval_file(filename = ramst_initval_file_data_col_vec_mat);
+
+perfect_foresight_setup(periods=200);
+if oo_.exo_simul(2) ~= 1.2
+  error('initval_file problem with exogenous variable');
+end
+if oo_.endo_simul(2, 2) ~= 13
+  error('initval_file problem with endogenous variable');
+end
+perfect_foresight_solver;
+if ~oo_.deterministic_simulation.status
+   error('Perfect foresight simulation failed')
+end
+
+if ispc()
+    initval_file(filename = ramst_initval_file_excel);
+    perfect_foresight_setup(periods=200);
+    perfect_foresight_solver;
+    if ~oo_.deterministic_simulation.status
+       error('Perfect foresight simulation failed');
+    end;
+end
+
+initval_file(datafile = 'ramst_data.m');
+perfect_foresight_setup(periods = 200);
+perfect_foresight_solver;
+if ~oo_.deterministic_simulation.status
+   error('Perfect foresight simulation failed')
+end
+oo_.exo_simul = [];
+oo_.endo_simul = [];
+
+initval_file(datafile = 'ramst_data.mat');
+perfect_foresight_setup(periods = 200);
+perfect_foresight_solver;
+if ~oo_.deterministic_simulation.status
+   error('Perfect foresight simulation failed')
+end
+oo_.exo_simul = [];
+oo_.endo_simul = [];
+
+initval_file(datafile = 'ramst_data.csv');
+perfect_foresight_setup(periods = 200);
+perfect_foresight_solver;
+if ~oo_.deterministic_simulation.status
+   error('Perfect foresight simulation failed')
+end
+oo_.exo_simul = [];
+oo_.endo_simul = [];
+
+initval_file(datafile = 'ramst_data.xlsx');
+perfect_foresight_setup(periods = 200);
+perfect_foresight_solver;
+if ~oo_.deterministic_simulation.status
+   error('Perfect foresight simulation failed')
+end
+oo_.exo_simul = [];
+oo_.endo_simul = [];
+
+if ispc;
+  initval_file(datafile = 'ramst_data.cls');
+  perfect_foresight_setup(periods = 200);
+  perfect_foresight_solver;
+  oo_.exo_simul = [];
+  oo_.endo_simul = [];
+
+  if ~oo_.deterministic_simulation.status
+    error('Perfect foresight simulation failed')
+  end;
+end;
+
diff --git a/tests/histval_initval_file/ramst_initval_file_data.m b/tests/histval_initval_file/ramst_initval_file_data.m
new file mode 100644
index 000000000..cc1914ca1
--- /dev/null
+++ b/tests/histval_initval_file/ramst_initval_file_data.m
@@ -0,0 +1,14 @@
+x = vertcat([ 1; 1.2 ], repmat(1, 200, 1));
+k = repmat(13, 202, 1);
+c = repmat(1.5, 202, 1);
+save('ramst_initval_file_data_col_vec_mat.mat','c','k','x');
+
+if ispc()
+    xlswrite('ramst_initval_file_excel',[x k c],1,'A2');
+    xlswrite('ramst_initval_file_excel',{'x' 'k' 'c'},1,'A1');
+end
+
+c=c';
+k=k';
+x=x';
+save('ramst_initval_file_data_row_vec_mat.mat','c','k','x');
diff --git a/tests/histval_initval_file/sim_exo_lead_lag.mod b/tests/histval_initval_file/sim_exo_lead_lag.mod
new file mode 100644
index 000000000..5a06c1c5c
--- /dev/null
+++ b/tests/histval_initval_file/sim_exo_lead_lag.mod
@@ -0,0 +1,43 @@
+// Uses autonomous system from sim_base.mod, but adds separate system where exogenous variables have several leads and lags
+// Lags and leads on exogenous variables are substituted out by auxiliary variables
+
+var c cmav k z_backward z_forward;
+varexo x;
+
+parameters alph gam delt bet aa;
+alph=0.5;
+gam=0.5;
+delt=0.02;
+bet=0.05;
+aa=0.5;
+
+model;
+  c + k - aa*x*k(-1)^alph - (1-delt)*k(-1);
+  c^(-gam) - (1+bet)^(-1)*(aa*alph*x(+1)*k^(alph-1) + 1 - delt)*c(+1)^(-gam);
+  z_backward=0.1*1+0.9*z_backward(-1) + (x(-4) - 1);
+  z_forward=0.2*1+0.8*z_forward(+1) + (x(+4) - 1);
+  cmav = 0.2*(c(-2) + c(-1) + c + c(+1) + c(+2));
+end;
+
+initval;
+  c = 1.2;
+  cmav = 1.2;
+  k = 12;
+  x = 1; //set x(0), x(-1), x(-2), x(-3)
+  z_backward = 1;
+  z_forward = 1;
+end;
+
+shocks;
+var x; //sets x(+2)
+periods 2;
+values 0.9;
+end;
+
+simul(periods=200);
+
+if ~oo_.deterministic_simulation.status
+   error('Perfect foresight simulation failed')
+end
+
+
diff --git a/tests/histval_initval_file/sim_exo_lead_lag_histvalf.mod b/tests/histval_initval_file/sim_exo_lead_lag_histvalf.mod
new file mode 100644
index 000000000..3d9ae12bb
--- /dev/null
+++ b/tests/histval_initval_file/sim_exo_lead_lag_histvalf.mod
@@ -0,0 +1,76 @@
+// Uses autonomous system from sim_base.mod, but adds separate system where exogenous variables have several leads and lags
+// Lags and leads on exogenous variables are substituted out by auxiliary variables
+
+var c cmav k z_backward z_forward;
+varexo x;
+
+parameters alph gam delt bet aa;
+alph=0.5;
+gam=0.5;
+delt=0.02;
+bet=0.05;
+aa=0.5;
+
+model;
+  c + k - aa*x*k(-1)^alph - (1-delt)*k(-1);
+  c^(-gam) - (1+bet)^(-1)*(aa*alph*x(+1)*k^(alph-1) + 1 - delt)*c(+1)^(-gam);
+  z_backward=0.1*1+0.9*z_backward(-1) + (x(-4) - 1);
+  z_forward=0.2*1+0.8*z_forward(+1) + (x(+4) - 1);
+  cmav = 0.2*(c(-2) + c(-1) + c + c(+1) + c(+2));
+end;
+
+initval;
+  x = 1;
+end;
+
+steady_state_model;
+  k = ((bet + delt)/(aa*alph*x))^(1/(alph - 1));
+  c = aa*x*k^alph - delt*k;
+  z_backward = x;
+  z_forward = x;
+  cmav = c;
+end;
+
+steady;
+
+shocks;
+  var x;
+  periods 2;
+  values 0.9;
+end;
+
+perfect_foresight_setup(periods=200);
+perfect_foresight_solver;
+
+reference = oo_.endo_simul;
+
+data1 = repmat([oo_.steady_state' 1], 4, 1);
+ds = dseries(data1, '1Y', [M_.endo_names; M_.exo_names]);
+
+histval_file(series = ds);
+perfect_foresight_setup(periods=200);
+perfect_foresight_solver;
+
+if ~oo_.deterministic_simulation.status
+   error('Perfect foresight simulation failed');
+end
+
+if max(max(abs(reference(1:5,5:end) - oo_.endo_simul(1:5,5:end)))) > 1e-8
+    error('Simulation with leads and lags doesn''t match the one with auxiliary variables')
+end
+
+data1 = repmat([oo_.steady_state' 1], 6, 1);
+ds1 = dseries(data1, '1Y', [M_.endo_names; M_.exo_names]);
+
+histval_file(series = ds1, first_obs = 3, last_obs = 6, nobs = 4);
+
+perfect_foresight_setup(periods=200);
+perfect_foresight_solver;
+
+if ~oo_.deterministic_simulation.status
+   error('Perfect foresight simulation failed');
+end
+
+if max(max(abs(reference(1:5,:) - oo_.endo_simul(1:5,:)))) > 1e-8
+    error('Simulation with leads and lags doesn''t match the one with auxiliary variables')
+end
diff --git a/tests/histval_initval_file/sim_exo_lead_lag_initvalf.mod b/tests/histval_initval_file/sim_exo_lead_lag_initvalf.mod
new file mode 100644
index 000000000..e74b6c9b2
--- /dev/null
+++ b/tests/histval_initval_file/sim_exo_lead_lag_initvalf.mod
@@ -0,0 +1,56 @@
+// Uses autonomous system from sim_base.mod, but adds separate system where exogenous variables have several leads and lags
+// Lags and leads on exogenous variables are substituted out by auxiliary variables
+
+data1 = repmat([1.2, 1.2, 12, 1, 1, 1], 208, 1);
+data1(6, 6) = 0.9; //shock to x in period 2
+ds = dseries(data1, '1Y', {'c', 'cmav', 'k', 'z_backward', 'z_forward', 'x'});
+
+var c cmav k z_backward z_forward;
+varexo x;
+
+parameters alph gam delt bet aa;
+alph=0.5;
+gam=0.5;
+delt=0.02;
+bet=0.05;
+aa=0.5;
+
+model;
+  c + k - aa*x*k(-1)^alph - (1-delt)*k(-1);
+  c^(-gam) - (1+bet)^(-1)*(aa*alph*x(+1)*k^(alph-1) + 1 - delt)*c(+1)^(-gam);
+  z_backward=0.1*1+0.9*z_backward(-1) + (x(-4) - 1);
+  z_forward=0.2*1+0.8*z_forward(+1) + (x(+4) - 1);
+  cmav = 0.2*(c(-2) + c(-1) + c + c(+1) + c(+2));
+end;
+
+initval_file(series = ds);
+
+perfect_foresight_setup(periods=200);
+perfect_foresight_solver(maxit=100);
+
+if ~oo_.deterministic_simulation.status
+   error('Perfect foresight simulation failed');
+end
+
+base_results=load('sim_exo_lead_lag_results.mat');
+if max(max(abs(base_results.oo_.endo_simul(1:5,:) - oo_.endo_simul(1:5,:)))) > 1e-8
+    error('Simulation with leads and lags doesn''t match the one with auxiliary variables')
+end
+
+data1 = repmat([1.2, 1.2, 12, 1, 1, 1], 212, 1);
+data1(8, 6) = 0.9; //shock to x in period 2
+ds1 = dseries(data1, '1Y', {'c', 'cmav', 'k', 'z_backward', 'z_forward', 'x'});
+
+initval_file(series = ds1, first_obs = 3, last_obs = 210, nobs = 208);
+
+perfect_foresight_setup(periods=200);
+perfect_foresight_solver(maxit=100);
+
+if ~oo_.deterministic_simulation.status
+   error('Perfect foresight simulation failed');
+end
+
+base_results=load('sim_exo_lead_lag_results.mat');
+if max(max(abs(base_results.oo_.endo_simul(1:5,:) - oo_.endo_simul(1:5,:)))) > 1e-8
+    error('Simulation with leads and lags doesn''t match the one with auxiliary variables')
+end
diff --git a/tests/histval_initval_file/sim_histvalf_stoch_simul.mod b/tests/histval_initval_file/sim_histvalf_stoch_simul.mod
new file mode 100644
index 000000000..924a5f002
--- /dev/null
+++ b/tests/histval_initval_file/sim_histvalf_stoch_simul.mod
@@ -0,0 +1,86 @@
+// Uses autonomous system from sim_base.mod, but adds separate system where exogenous variables have several leads and lags
+// Lags and leads on exogenous variables are substituted out by auxiliary variables
+
+var c cmav k z_backward z_forward;
+varexo x;
+
+parameters alph gam delt bet aa;
+alph=0.5;
+gam=0.5;
+delt=0.02;
+bet=0.05;
+aa=0.5;
+
+model;
+  c + k - aa*x*k(-1)^alph - (1-delt)*k(-1);
+  c^(-gam) - (1+bet)^(-1)*(aa*alph*x(+1)*k^(alph-1) + 1 - delt)*c(+1)^(-gam);
+  z_backward=0.1*1+0.9*z_backward(-1) + (x(-4) - 1);
+  z_forward=0.2*1+0.8*z_forward(+1) + (x(+4) - 1);
+  cmav = 0.2*(c(-2) + c(-1) + c + c(+1) + c(+2));
+end;
+
+initval;
+  x = 1;
+end;
+
+steady_state_model;
+  k = ((bet + delt)/(aa*alph*x))^(1/(alph - 1));
+  c = aa*x*k^alph - delt*k;
+  z_backward = x;
+  z_forward = x;
+  cmav = c;
+end;
+
+steady;
+
+shocks;
+  var x;
+  stderr 0.01;
+end;
+
+s = rng;
+stoch_simul(periods=20, drop=0, irf=0);
+
+reference = oo_.endo_simul;
+
+data1 = repmat([oo_.steady_state' 1], 4, 1);
+ds = dseries(data1, '1Y', [M_.endo_names; M_.exo_names]);
+
+histval_file(series = ds);
+
+rng(s);
+stoch_simul(periods=20, drop=0, irf=0);
+
+if max(max(abs(reference(1:5,5:end) - oo_.endo_simul(1:5,5:end)))) > 1e-8
+    error('Simulation with leads and lags doesn''t match the reference')
+end
+
+data1 = repmat([oo_.steady_state' 1], 6, 1);
+ds1 = dseries(data1, '1Y', [M_.endo_names; M_.exo_names]);
+
+histval_file(series = ds1, first_obs = 6, last_obs = 6, nobs = 1);
+
+rng(s);
+stoch_simul(periods=20, drop=0, irf=0);
+
+if max(max(abs(reference(1:5,:) - oo_.endo_simul(1:5,:)))) > 1e-8
+    error('Simulation with leads and lags doesn''t match the reference')
+end
+
+histval_file(series = ds1, first_simulation_period = 7);
+
+rng(s);
+stoch_simul(periods=20, drop=0, irf=0);
+
+if max(max(abs(reference(1:5,:) - oo_.endo_simul(1:5,:)))) > 1e-8
+    error('Simulation with leads and lags doesn''t match the reference')
+end
+
+histval_file(series = ds1, first_simulation_period = 7Y);
+
+rng(s);
+stoch_simul(periods=20, drop=0, irf=0);
+
+if max(max(abs(reference(1:5,:) - oo_.endo_simul(1:5,:)))) > 1e-8
+    error('Simulation with leads and lags doesn''t match the reference')
+end
diff --git a/tests/histval_initval_file_unit_tests.m b/tests/histval_initval_file_unit_tests.m
new file mode 100644
index 000000000..4b81f71de
--- /dev/null
+++ b/tests/histval_initval_file_unit_tests.m
@@ -0,0 +1,220 @@
+top_test_dir = getenv('TOP_TEST_DIR');
+addpath([top_test_dir filesep '..' filesep 'matlab/']);
+dynare_config;
+
+cd('histval_initval_file');
+num_tests = 0;
+failed_tests = {};
+
+ds = dseries(randn(10,4));
+
+M = struct();
+M.fname = '';
+M.endo_nbr = 3;
+M.orig_endo_nbr = 3;
+M.endo_names = {'Variable_1','Variable_2','Variable_3'};
+M.exo_nbr = 1;
+M.exo_names = {'Variable_4'};
+M.exo_det_nbr = 0;
+
+caller = 'INITVAL';
+
+options = struct();
+options.series = ds;
+ds1 = histvalf_initvalf(caller, M, options);
+
+failed_tests = my_assert(failed_tests, all(all(ds1 == ds)), 'basic test');
+num_tests = num_tests + 1;
+
+options = struct();
+options.series = ds;
+options.first_obs = 2;
+ds1 = histvalf_initvalf(caller, M, options);
+failed_tests = my_assert(failed_tests, ds1.init == dates('2Y'), ...
+                         'init test 1');
+num_tests = num_tests + 1;
+
+options = struct();
+options.series = ds;
+options.first_obs = 2;
+options.last_obs = 9;
+ds1 = histvalf_initvalf(caller, M, options);
+failed_tests = my_assert(failed_tests, ds1.init == dates('2Y'), ...
+                         'first_obs last_obs test 1');
+failed_tests = my_assert(failed_tests, ds1.last == dates('9Y'), ...
+                         'first_obs last_obs test 2');
+num_tests = num_tests + 2;
+
+options = struct();
+options.series = ds;
+options.last_obs = 9;
+ds1 = histvalf_initvalf(caller, M, options);
+failed_tests = my_assert(failed_tests, ds1.init == dates('1Y'), ...
+                         'last_obs test 1');
+failed_tests = my_assert(failed_tests, ds1.last == dates('9Y'), ...
+                         'last_obs test 2');
+num_tests = num_tests + 2;
+
+options = struct();
+options.series = ds;
+options.first_obs = 2;
+options.last_obs = 9;
+options.nobs = 8;
+ds1 = histvalf_initvalf(caller, M, options);
+failed_tests = my_assert(failed_tests, ds1.init == dates('2Y'), ...
+                         'first_obs, last_obs, nobs test 1');
+failed_tests = my_assert(failed_tests, ds1.last == dates('9Y'), ...
+                         'first_obs, last_obs, nobs test 2');
+num_tests = num_tests + 2;
+
+options = struct();
+options.series = ds;
+options.last_obs = 9;
+options.nobs = 8;
+ds1 = histvalf_initvalf(caller, M, options);
+failed_tests = my_assert(failed_tests, ds1.init == dates('2Y'), ...
+                         'last_obs, nobs test 1');
+failed_tests = my_assert(failed_tests, ds1.last == dates('9Y'), ...
+                         'last_obs, nobs test 2');
+num_tests = num_tests + 2;
+
+options = struct();
+options.series = ds;
+options.first_obs = 2;
+options.last_obs = 9;
+options.nobs = 7;
+
+try
+    ds1 = histvalf_initvalf(caller, M, options);
+    error('This test didn''t catch the error')
+catch me
+    if strcmp(me.message, ['INITVAL_FILE: FIST_OBS, LAST_OBS and NOBS contain', ...
+                           ' inconsistent information. Use only two of these', ...
+                           ' options.']) == false
+        failed_tests = cat(1, failed_tests, 'Wrong nobs error message' );
+    end
+end
+num_tests = num_tests + 1;
+
+options = struct();
+options.series = ds;
+options.first_obs = -1;
+
+try
+    ds1 = histvalf_initvalf(caller, M, options);
+    error('This test didn''t catch the error')
+catch me
+    if strcmp(me.message, [caller, '_FILE: the first requested period is', ...
+                           ' before available data.']) == false
+        failed_tests = cat(1, failed_tests, ...
+                           'Wrong first period error message');
+    end
+end
+num_tests = num_tests + 1;
+
+options = struct();
+options.series = ds;
+options.last_obs = 11;
+
+try
+    ds1 = histvalf_initvalf(caller, M, options);
+    error('This test didn''t catch the error')
+catch me
+    if strcmp(me.message, [caller, '_FILE: the last requested period is', ...
+                           ' after available data.']) == false
+        failed_tests = cat(1, failed_tests, ...
+                           'Wrong last period error message');
+    end
+end
+num_tests = num_tests + 1;
+
+fh = fopen('data.m', 'w');
+init__ = 'INIT__ = ''1Y'';';
+fprintf(fh, [init__ '\n']);
+eval(init__);
+names__ = 'NAMES__ = {''x'', ''y''};';
+fprintf(fh, [names__ '\n']);
+eval(names__);
+tex__ = 'TEX__ = {''x'', ''y''};';
+fprintf(fh, [tex__ '\n']);
+eval(tex__);
+x = randn(10, 1);
+fprintf(fh, 'x = [');
+fprintf(fh, '%f ', x);
+fprintf(fh, '];\n');
+y = randn(10, 1);
+fprintf(fh, 'y = [');
+fprintf(fh, '%f ', y);
+fprintf(fh, '];\n');
+fclose(fh);
+
+M.endo_nbr = 1;
+M.orig_endo_nbr = 1;
+M.endo_names = {'y'};
+M.exo_nbr = 1;
+M.exo_names = {'x'};
+M.exo_det_nbr = 0;
+
+options = struct();
+options.datafile = 'data.m';
+series = histvalf_initvalf('INITVAL_FILE', M, options);
+failed_tests = my_assert(failed_tests, series.init == dates('1Y'), ...
+                         '*.m file first_obs test');
+failed_tests = my_assert(failed_tests, series.nobs == 10, ...
+                         '*.m file nobs test');
+
+save('data.mat', 'INIT__', 'NAMES__', 'TEX__', 'x', 'y');
+options = struct();
+options.datafile = 'data.mat';
+series = histvalf_initvalf('INITVAL_FILE', M, options);
+failed_tests = my_assert(failed_tests, series.init == dates('1Y'), ...
+                         '*.mat file first_obs test');
+failed_tests = my_assert(failed_tests, series.nobs == 10, ...
+                         '*.mat file nobs test');
+
+fh = fopen('data.csv', 'w');
+fprintf(fh, 'x,y\n');
+for i = 1:size(x,1)
+    fprintf(fh, '%f,%f\n', x(i), y(i));
+end
+fclose(fh);
+
+if ~verLessThan('matlab', '8.2')
+    writetable(table(x,y), 'data.xlsx')
+    options = struct();
+    options.datafile = 'data.xlsx';
+    series = histvalf_initvalf('INITVAL_FILE', M, options);
+    failed_tests = my_assert(failed_tests, series.init == dates('1Y'), ...
+                             '*.xlsx file first_obs test');
+    failed_tests = my_assert(failed_tests, series.nobs == 10, ...
+                             '*.xlsx file nobs test');
+    num_tests = num_tests + 2;
+
+    if ispc
+        writetable(table(x,y), 'data.xls')
+        options = struct();
+        options.datafile = 'data.xls';
+        series = histvalf_initvalf('INITVAL_FILE', M, options);
+        failed_tests = my_assert(failed_tests, series.init == dates('1Y'), ...
+                             '*.xls file first_obs test');
+        failed_tests = my_assert(failed_tests, series.nobs == 10, ...        
+                             '*.xls file nobs test');
+        num_tests = num_tests + 2;
+    end
+end
+
+cd(getenv('TOP_TEST_DIR'));
+fid = fopen('histval_initval_file_unit_tests.m.trs', 'w+');
+num_failed_tests = length(failed_tests)
+if num_failed_tests > 0
+  fprintf(fid,':test-result: FAIL\n');
+  fprintf(fid,':number-tests: %d\n', num_tests);
+  fprintf(fid,':number-failed-tests: %d\n', num_failed_tests);
+  fprintf(fid,':list-of-failed-tests: %s\n', failed_tests{:});
+else
+  fprintf(fid,':test-result: PASS\n');
+  fprintf(fid,':number-tests: %d\n', num_tests);
+  fprintf(fid,':number-failed-tests: 0\n');
+end
+fclose(fid);
+exit;
diff --git a/tests/initval_file/ramst_initval_file.mod b/tests/initval_file/ramst_initval_file.mod
deleted file mode 100644
index 2cc8f91ed..000000000
--- a/tests/initval_file/ramst_initval_file.mod
+++ /dev/null
@@ -1,55 +0,0 @@
-/* Test for the initval_file() command. This file needs ramst_initval_file_data.m. It should give results similar to those of ramst.mod */
-
-var c k;
-varexo x;
-
-parameters alph gam delt bet aa;
-alph=0.5;
-gam=0.5;
-delt=0.02;
-bet=0.05;
-aa=0.5;
-
-
-model;
-c + k - aa*x*k(-1)^alph - (1-delt)*k(-1);
-c^(-gam) - (1+bet)^(-1)*(aa*alph*x(+1)*k^(alph-1) + 1 - delt)*c(+1)^(-gam);
-end;
-
-initval;
-x = 1;
-k = ((delt+bet)/(1.0*aa*alph))^(1/(alph-1));
-c = aa*k^alph-delt*k;
-end;
-
-steady;
-
-initval_file(filename = ramst_initval_file_data_row_vec_mat);
-if oo_.exo_simul(2) ~= 1.2
-  error('initval_file problem with exogenous variable');
-end
-if oo_.endo_simul(2, 2) ~= 13
-  error('initval_file option problem with endogenous variable');
-end
-perfect_foresight_setup(periods=200);
-perfect_foresight_solver;
-
-oo_.exo_simul = [];
-oo_.endo_simul = [];
-
-initval_file(filename = ramst_initval_file_data_col_vec_mat);
-if oo_.exo_simul(2) ~= 1.2
-  error('initval_file problem with exogenous variable');
-end
-if oo_.endo_simul(2, 2) ~= 13
-  error('initval_file problem with endogenous variable');
-end
-
-perfect_foresight_setup(periods=200);
-perfect_foresight_solver;
-
-if ispc()
-    initval_file(filename = ramst_initval_file_excel);
-    perfect_foresight_setup(periods=200);
-    perfect_foresight_solver;
-end
diff --git a/tests/smoother2histval/fs2000_simul.mod b/tests/smoother2histval/fs2000_simul.mod
index eb65974e1..a682c76c0 100644
--- a/tests/smoother2histval/fs2000_simul.mod
+++ b/tests/smoother2histval/fs2000_simul.mod
@@ -66,6 +66,10 @@ results_estimation=load('fs2000_smooth_results');
 M_.params=results_estimation.M_.params;
 steady;
 
+OO = load('fs2000_smooth_results.mat');
+M_.params = OO.M_.params;
+
 histval_file(filename = 'fs2000_histval.mat');
 
-simul(periods = 30);
+perfect_foresight_setup(periods = 100);
+perfect_foresight_solver;
diff --git a/tests/smoother2histval/fs2000_smooth.mod b/tests/smoother2histval/fs2000_smooth.mod
index 89447eb57..65bb15b68 100644
--- a/tests/smoother2histval/fs2000_smooth.mod
+++ b/tests/smoother2histval/fs2000_smooth.mod
@@ -82,6 +82,6 @@ varobs gp_obs gy_obs;
 
 options_.solve_tolf = 1e-12;
 
-estimation(order=1,datafile=fsdat_simul,nobs=192,mh_replic=1500,mh_nblocks=1,mh_jscale=0.8,smoother,consider_all_endogenous);
+estimation(order=1,datafile=fsdat_simul,mh_replic=1500,mh_nblocks=1,mh_jscale=0.8,smoother,consider_all_endogenous);
 
 smoother2histval(period = 5, outfile = 'fs2000_histval.mat');