diff --git a/tests/Makefile.am b/tests/Makefile.am
index ea1725e7c..924772063 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -214,7 +214,8 @@ MODFILES = \
 	ep/rbc2.mod \
 	ep/rbcii.mod \
 	ep/linearmodel.mod \
-	deterministic_simulations/deterministic_model_purely_forward.mod \
+	deterministic_simulations/purely_forward/ar1.mod \
+        deterministic_simulations/purely_forward/nk.mod \
 	deterministic_simulations/purely_backward/ar1.mod \
 	deterministic_simulations/rbc_det1.mod \
 	deterministic_simulations/rbc_det2.mod \
diff --git a/tests/deterministic_simulations/deterministic_model_purely_forward.mod b/tests/deterministic_simulations/deterministic_model_purely_forward.mod
deleted file mode 100644
index c13e5e055..000000000
--- a/tests/deterministic_simulations/deterministic_model_purely_forward.mod
+++ /dev/null
@@ -1,54 +0,0 @@
-var y i pi rbar ;
-
-varexo r tauw taus taua gn;
-
-parameters khia khiw khis phipi phiy taubs taubw tauba w sigma psi kappa alpha mu beta teta;
-
-teta = 12.7721;
-sigma = 1.1599;
-beta = 0.9970;
-alpha = 0.7747;
-mu = 0.9030;
-taubs = 0.05;
-taubw = 0.02;
-tauba = 0;
-w = 1.5692;
-phipi = 1.5;
-phiy = 0.5/4;
-
-khia = (1-beta)/(1-tauba);
-khiw = 1/(1-taubw);
-khis = 1/(1+taubs);
-psi = 1/(sigma + w);
-kappa = (1-alpha)*(1-alpha*beta)*(sigma+w)/(alpha*(1+w*teta));
-
-model(linear);
-y = y(+1)-sigma*(i-pi(+1)-r)+(gn-gn(+1))+(sigma)^-1*khis*(taus(+1)-taus)+sigma*khia*taua;
-
-pi=kappa*y+kappa*psi*(khiw*tauw+khis*taus-sigma*gn)+beta*pi(+1);
-
-i=max(0,r+phipi*pi+phiy*y);
-
-rbar = -((kappa*phipi+(1-beta*mu)*phiy)*sigma^-1*khia*taus)/((1-mu+sigma^-1*phiy)*(1-beta*mu)+kappa*sigma^-1*(phipi-mu))
-- (((1-mu)*kappa*psi*phipi+sigma^-1*mu*kappa*psi*phiy)*khiw*tauw)/((1-mu+sigma^-1*phiy)*(1-beta*mu)+kappa*sigma*(phipi-mu))
--(kappa*sigma*(1-mu)*(sigma^-1-psi)*phipi+((1-mu)*(1-beta*mu)-kappa*psi*mu)*phiy)*(gn-sigma^-1*khis*taus)/((1-mu-sigma^-1*phiy)*(1-beta*mu)+kappa*sigma^-1*(phipi-mu));
-
-end;
-
-initval;
-y=0;
-i=-log(beta);
-pi=0;
-rbar = 0;
-end;
-
-steady;
-check;
-
-shocks;
-var r;
-periods 1:9;
-values -0.0104;
-end;
-
-simul(periods=2100);
diff --git a/tests/deterministic_simulations/purely_forward/ar1.mod b/tests/deterministic_simulations/purely_forward/ar1.mod
new file mode 100644
index 000000000..0942e8466
--- /dev/null
+++ b/tests/deterministic_simulations/purely_forward/ar1.mod
@@ -0,0 +1,33 @@
+var y;
+
+varexo eps;
+
+model;
+    y = sqrt(y(1))*exp(eps);
+end;
+
+initval;
+    y = 1;
+    eps = 0;
+end;
+
+steady;
+
+check;
+
+
+shocks;
+    var eps;
+    periods 1 2;
+    values 1 -1;
+end;
+
+simul(periods=5);
+
+expected_y = ones(1, 6);
+expected_y(2) = exp(-1);
+expected_y(1) = sqrt(exp(-1))*exp(1);
+
+if max(abs(oo_.endo_simul-expected_y))>options_.dynatol.x
+    error('Wrong solution!')
+end
diff --git a/tests/deterministic_simulations/purely_forward/nk.mod b/tests/deterministic_simulations/purely_forward/nk.mod
new file mode 100644
index 000000000..5d5dcfa8b
--- /dev/null
+++ b/tests/deterministic_simulations/purely_forward/nk.mod
@@ -0,0 +1,76 @@
+var y i pi ;
+
+varexo e;
+
+parameters a1 a2 a3 a4 a5;
+
+a1 = -.5;
+a2 =  .1;
+a3 =  .9;
+a4 = 1.5;
+a5 = 0.5;
+
+model;
+
+y = y(1)+a1*(i-pi(1))+e;
+
+pi = a2*y + a3*pi(1);
+
+i = max(0, a4*pi+a5*y);
+
+end;
+
+steady_state_model;
+y=0;
+i=0;
+pi=0;
+end;
+
+steady;
+check;
+
+shocks;
+var e;
+periods 1 2;
+values  .3  -0.1;
+end;
+
+simul(periods=5);
+
+% Initialize the analytical solution for the endogenous variables.
+expected_y = zeros(1, 6);
+expected_pi = zeros(1, 6);
+expected_i = zeros(1, 6);
+
+% Set period 2.
+tmp = inv(eye(3)-[0 0 a1; a2 0 0; a5 a4 0])*[oo_.exo_simul(2); 0; 0];
+
+if tmp(3)>=0
+    expected_y(2) = tmp(1);
+    expected_pi(2) = tmp(2);
+    expected_i(2) = tmp(3);
+else
+    expected_y(2) = oo_.exo_simul(2);
+    expected_pi(2) = expected_y(2)*a2;
+    expected_i(2) = 0;
+end
+
+% Set period 1.
+tmp = inv(eye(3)-[0 0 a1; a2 0 0; a5 a4 0])*[oo_.exo_simul(1)+expected_y(2)-a1*expected_pi(2); a3*expected_pi(2); 0];
+
+if tmp(3)>=0
+    expected_y(1) = tmp(1);
+    expected_pi(1) = tmp(2);
+    expected_i(1) = tmp(3);
+else
+    expected_y(1) = oo_.exo_simul(1)+expected_y(2)-a1*expected_pi(2);
+    expected_pi(1) = expected_y(2)*a2+a3*expected_pi(2);
+    expected_i(1) = 0;
+end
+
+% Compare the paths returned by sim1_purely_forward routine and the analytical solution.
+if max(abs(oo_.endo_simul(1,:)-expected_y))>options_.dynatol.x || ...
+     max(abs(oo_.endo_simul(2,:)-expected_i))>options_.dynatol.x || ...
+     max(abs(oo_.endo_simul(3,:)-expected_pi))>options_.dynatol.x
+    error('Wrong solution!')
+end