diff --git a/matlab/histvalf.m b/matlab/histvalf.m
index b77968c54..59c98cbee 100644
--- a/matlab/histvalf.m
+++ b/matlab/histvalf.m
@@ -25,6 +25,16 @@ 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);
+    else
+        oo_.exo_simul = repmat(ex0_',M_.maximum_lag,1);
+    end
+end
+
 S = load(fname);
 
 outvars = fieldnames(S);
@@ -35,26 +45,43 @@ for i = 1:length(outvars)
         ov = ov_(1:end-1);
         j = strmatch(ov, M_.endo_names, 'exact');
         if isempty(j)
-            error(['smoother2histval: output variable ' ov ' does not exist.'])
+            warning(['smoother2histval: output variable ' ov ' does not exist.'])
         end
     else
-        % Lagged endogenous, search through aux vars
-        z = strsplit(ov_, '_');
-        ov = z{1};
-        lead_lag = str2num(z{2});
+        % 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 = str2num(ov_((undidx+1):end));
         j = [];
         for i = 1:length(M_.aux_vars)
-            if M_.aux_vars(i).type ~= 1 
+            if M_.aux_vars(i).type ~= 1 && M_.aux_vars(i).type ~= 3
                 continue
             end
-            orig_var = deblank(M_.endo_names(M_.aux_vars(i).orig_index, :));
+            if M_.aux_vars(i).type == 1
+                % Endogenous
+                orig_var = deblank(M_.endo_names(M_.aux_vars(i).orig_index, :));
+            else
+                % Exogenous
+                orig_var = deblank(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) = getfield(S, ov_);
+            end
             continue
         end
     end
     M_.endo_histval(j, :) = getfield(S, ov_);
 end
+
diff --git a/matlab/smoother2histval.m b/matlab/smoother2histval.m
index 87e50c821..bfdcab373 100644
--- a/matlab/smoother2histval.m
+++ b/matlab/smoother2histval.m
@@ -1,11 +1,11 @@
 function smoother2histval(opts)
-% This function takes values from oo_.SmoothedVariables and copies them into
-% M_.histval.
+% This function takes values from oo_.SmoothedVariables (and possibly
+% oo_.SmoothedShocks) and copies them into M_.histval.
 % 
 % Optional fields in 'opts' structure:
 %    infile:      An optional *_results MAT file created by Dynare.
-%                 If present, oo_.SmoothedVariables is read from there.
-%                 Otherwise, it is read from the global workspace.
+%                 If present, oo_.Smoothed{Variables,Shocks} are read from 
+%                 there. Otherwise, they are read from the global workspace.
 %    invars:      An optional char or cell array listing variables to read in
 %                 oo_.SmoothedVariables. If absent, all the endogenous
 %                 variables present in oo_.SmoothedVariables are used.
@@ -47,17 +47,19 @@ if ~isfield(opts, 'infile')
     if ~isfield(oo_, 'SmoothedVariables')
         error('Could not find smoothed variables; did you set the "smoother" option?')
     end
-    smoothedvals = oo_.SmoothedVariables;
+    smoothedvars = oo_.SmoothedVariables;
+    smoothedshocks = oo_.SmoothedShocks;
 else
     S = load(opts.infile);
     if ~isfield(S, 'oo_') || ~isfield(S.oo_, 'SmoothedVariables')
         error('Could not find smoothed variables in file; is this a Dynare results file, and did you set the "smoother" option when producing it?')
     end
-    smoothedvals = S.oo_.SmoothedVariables;
+    smoothedvars = S.oo_.SmoothedVariables;
+    smoothedshocks = S.oo_.SmoothedShocks;
 end
 
 % Hack to determine if oo_.SmoothedVariables was computed after a Metropolis
-if isstruct(getfield(smoothedvals, fieldnames(smoothedvals){1}))
+if isstruct(getfield(smoothedvars, fieldnames(smoothedvars){1}))
     post_metropolis = 1;
 else
     post_metropolis = 0;
@@ -66,7 +68,8 @@ end
 % If post-Metropolis, select the parameter set
 if isempty(options_.parameter_set)
     if post_metropolis
-        smoothedvals = smoothedvals.Mean;
+        smoothedvars = smoothedvars.Mean;
+        smoothedshocks = smoothedshocks.Mean;
     end
 else
     switch options_.parameter_set
@@ -82,19 +85,21 @@ else
         if ~post_metropolis
             error('Option parameter_set=posterior_mean is not consistent with computed smoothed values.')
         end
-        smoothedvals = smoothedvals.Mean;
+        smoothedvars = smoothedvars.Mean;
+        smoothedshocks = smoothedshocks.Mean;
       case 'posterior_median'
         if ~post_metropolis
             error('Option parameter_set=posterior_median is not consistent with computed smoothed values.')
         end
-        smoothedvals = smoothedvals.Median;
+        smoothedvars = smoothedvars.Median;
+        smoothedshocks = smoothedshocks.Median;
       otherwise
         error([ 'Option parameter_set=' options_.parameter_set ' unsupported.' ])
     end
 end
 
 % Determine number of periods
-n = size(getfield(smoothedvals, fieldnames(smoothedvals){1}));
+n = size(getfield(smoothedvars, fieldnames(smoothedvars){1}));
 
 if n < M_.maximum_endo_lag
     error('Not enough observations to create initial conditions')
@@ -106,7 +111,7 @@ if isfield(opts, 'invars')
         invars = cellstr(invars);
     end
 else
-    invars = fieldnames(smoothedvals);
+    invars = [fieldnames(smoothedvars); fieldnames(smoothedshocks)];
 end
 
 if isfield(opts, 'period')
@@ -142,9 +147,13 @@ else
     o = struct();
 end
 
-% Handle all variables to be copied
+% Handle all endogenous variables to be copied
 for i = 1:length(invars)
-    s = getfield(smoothedvals, invars{i});
+    if isempty(strmatch(invars{i}, M_.endo_names))
+        % Skip exogenous
+        continue
+    end
+    s = getfield(smoothedvars, invars{i});
     v = s((period-M_.maximum_endo_lag+1):period);
     if ~isfield(opts, 'outfile')
         j = strmatch(outvars{i}, M_.endo_names, 'exact');
@@ -159,15 +168,25 @@ for i = 1:length(invars)
     end
 end
 
-% Handle auxiliary variables for lags
+% Handle auxiliary variables for lags (both on endogenous and exogenous)
 for i = 1:length(M_.aux_vars)
-    if M_.aux_vars(i).type ~= 1
+    if M_.aux_vars(i).type ~= 1 && M_.aux_vars(i).type ~= 3
         continue
     end
-    orig_var = deblank(M_.endo_names(M_.aux_vars(i).orig_index, :));
+    if M_.aux_vars(i).type == 1
+        % Endogenous
+        orig_var = deblank(M_.endo_names(M_.aux_vars(i).orig_index, :));
+    else
+        % Exogenous
+        orig_var = deblank(M_.exo_names(M_.aux_vars(i).orig_index, :));
+    end
     [m, k] = ismember(orig_var, outvars);
     if m
-        s = getfield(smoothedvals, invars{k});
+        if ~isempty(strmatch(invars{k}, M_.endo_names))
+            s = getfield(smoothedvars, invars{k});
+        else
+            s = getfield(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')
diff --git a/tests/smoother2histval/fs2000_simul.mod b/tests/smoother2histval/fs2000_simul.mod
index 0a7c018f4..3dcfc5f60 100644
--- a/tests/smoother2histval/fs2000_simul.mod
+++ b/tests/smoother2histval/fs2000_simul.mod
@@ -26,7 +26,7 @@ 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));
+y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a(-1)));
 gy_obs = dA*y/y(-2);
 gp_obs = (P/P(-1))*m(-1)/dA;
 end;
diff --git a/tests/smoother2histval/fs2000_smooth.mod b/tests/smoother2histval/fs2000_smooth.mod
index 26d3e7249..5f2931f28 100644
--- a/tests/smoother2histval/fs2000_smooth.mod
+++ b/tests/smoother2histval/fs2000_smooth.mod
@@ -26,7 +26,7 @@ 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));
+y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a(-1)));
 gy_obs = dA*y/y(-2);
 gp_obs = (P/P(-1))*m(-1)/dA;
 end;
@@ -73,6 +73,6 @@ varobs gp_obs gy_obs;
 
 options_.solve_tolf = 1e-12;
 
-estimation(order=1,datafile=fsdat_simul,nobs=192,loglinear,mh_replic=1500,mh_nblocks=1,mh_jscale=0.8,smoother);
+estimation(order=1,datafile=fsdat_simul,nobs=192,mh_replic=1500,mh_nblocks=1,mh_jscale=0.8,smoother,consider_all_endogenous);
 
 smoother2histval(period = 5, outfile = 'fs2000_histval.mat');