diff --git a/doc/dynare.texi b/doc/dynare.texi
index 90e02d11208b3eecc7c697c4c297ec2ce6a8697b..6ea52531a8f3567160ff43506d7f01a81737fb97 100644
--- a/doc/dynare.texi
+++ b/doc/dynare.texi
@@ -6419,7 +6419,7 @@ parameters in the @code{params}-command and be entered in the
 
 @item
 @math{e} are the exogenous stochastic shocks, specified in the
-@code{var_exo}-command;
+@code{varexo}-command;
 
 @item
 @math{W} is the weighting matrix;
diff --git a/matlab/stochastic_solvers.m b/matlab/stochastic_solvers.m
index 8f22a025ce62717203399927590206bba79dfc41..64b9b2c810933f2239c54e60f163334703f582c9 100644
--- a/matlab/stochastic_solvers.m
+++ b/matlab/stochastic_solvers.m
@@ -241,9 +241,19 @@ else
         dr = dyn_second_order_solver(jacobia_,hessian1,dr,M_,...
                                      options_.threads.kronecker.A_times_B_kronecker_C,...
                                      options_.threads.kronecker.sparse_hessian_times_B_kronecker_C);
+                                 
+        % reordering second order derivatives, used for deterministic
+        % variables below
+        k1 = nonzeros(M_.lead_lag_incidence(:,order_var)');
+        kk = [k1; length(k1)+(1:M_.exo_nbr+M_.exo_det_nbr)'];
+        nk = size(kk,1);
+        kk1 = reshape([1:nk^2],nk,nk);
+        kk1 = kk1(kk,kk);
+        hessian1 = hessian1(:,kk1(:));
     end
 end
 
+
 %exogenous deterministic variables
 if M_.exo_det_nbr > 0
     gx = dr.gx;
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 271da5c1565de168a6d7b012ca17dcaf0d3b6fe1..1dbce930e1c00f107f11f678075f0ed8f69e106b 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -180,6 +180,7 @@ MODFILES = \
 	shock_decomposition/fs2000_est_varlist.mod \
 	stochastic_purely_forward/stochastic_purely_forward.mod \
 	stochastic_purely_forward/stochastic_purely_forward_with_static.mod \
+	forecast/Hansen_exo_det_forecast.mod \
 	gradient/fs2000_numgrad_13.mod \
 	gradient/fs2000_numgrad_15.mod \
 	gradient/fs2000_numgrad_2.mod \
diff --git a/tests/forecast/Hansen_exo_det_forecast.mod b/tests/forecast/Hansen_exo_det_forecast.mod
new file mode 100644
index 0000000000000000000000000000000000000000..22744908a5c8534855ef981ea5e2ce0f3c658677
--- /dev/null
+++ b/tests/forecast/Hansen_exo_det_forecast.mod
@@ -0,0 +1,95 @@
+/*
+* The Hansen model following McCandless, George T. (2008): The ABCs of RBCs, Hardvard University Press, Chapter 6
+*
+* This mod-file tests the correctness of forecasting with exogenous deterministic variables.
+* A forecast starting at the steady state in t=0, where the only shock is a perfectly 
+* anticipated shock in t=8, is equal to the IRF to a 7 period anticipated news shock.
+* Note the timing difference due to the fact that in forecasting. the agent starts at the 
+* steady state at time 0 and has the first endogenous reaction period at t=1 so that the shock at 
+* t=8 is only 7 period anticipated
+
+* This implementation was written by Johannes Pfeifer. Please note that the
+* following copyright notice only applies to this Dynare implementation of the
+* model.
+*/
+
+/*
+ * Copyright (C) 2014 Johannes Pfeifer
+ *
+ * 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/>.
+ */
+
+
+
+var c r y h k a;
+varexo eps_a_surp eps_a_antic;
+varexo_det e_det;
+
+parameters beta delta theta rho_a A_disutil;
+
+beta = 0.99;
+delta = 0.025;
+theta = 0.36;
+rho_a = 0.35;
+A_disutil = 2;
+
+model;
+1/c = beta*((1/c(+1))*(r(+1) +(1-delta)));
+(1-h)*(1-theta)*(y/h) = A_disutil*c;
+c = y +(1-delta)*k(-1) - k;
+y = a*k(-1)^(theta)*h^(1-theta);
+r = theta*(y/k(-1));
+log(a)=rho_a*log(a(-1))+eps_a_surp+eps_a_antic(-7)+e_det;
+end;
+
+steady_state_model;
+a = 1;
+h = (1+(A_disutil/(1-theta))*(1 - (beta*delta*theta)/(1-beta*(1-delta))))^(-1);
+k = h*(theta*a/(1/beta -(1-delta)))^(1/(1-theta));
+y = a*k^(theta)*h^(1-theta);
+c = y-delta*k;
+r =  1/beta - (1-delta);
+end;
+
+steady;
+
+shocks;
+var eps_a_surp; stderr 0;
+var eps_a_antic; stderr 0.01;
+var e_det;
+periods 8;
+values 0.01;
+end;
+
+stoch_simul(irf=40,nomoments, order=1);
+
+% make_ex_;
+
+forecast(periods=40);
+figure
+for ii=1:M_.orig_endo_nbr
+    subplot(3,3,ii)    
+    var_name=deblank(M_.endo_names(ii,:));
+    var_index=strmatch(var_name,deblank(M_.endo_names),'exact');
+    plot(1:40,oo_.dr.ys(var_index)+oo_.irfs.([var_name,'_eps_a_antic']),'b',1:40,... %
+            oo_.forecast.Mean.(var_name),'r--')
+    title(var_name)
+    difference(ii,:)=oo_.dr.ys(var_index)+oo_.irfs.([var_name,'_eps_a_antic'])-oo_.forecast.Mean.(var_name)';
+end
+
+if max(max(abs(difference)))>1e-10
+   error('Forecasts with exogenous deterministic variable is wrong') 
+end
\ No newline at end of file