diff --git a/test/models/herbst_schorfheide/nk_est.mod b/test/models/herbst_schorfheide/nk_est.mod
new file mode 100644
index 00000000..43336059
--- /dev/null
+++ b/test/models/herbst_schorfheide/nk_est.mod
@@ -0,0 +1,88 @@
+var y R pi z g YGR INFL INT;
+varexo eg ep ez;
+
+parameters tau rho_R psi1 psi2 rho_g rho_z kappa  piA gammaQ rA s_ep s_eg s_ez; 
+
+rA = 0.42;
+piA = 3.3;
+gammaQ = 0.52;
+tau = 2.83;
+kappa = 0.78;
+psi1 = 1.8;
+psi2 = 0.63;
+rho_R = 0.77;
+rho_g = 0.98;
+rho_z = 0.88;
+s_ep = 0.22;
+s_eg = 0.72;
+s_ez = 0.31;
+
+
+model;
+#  beta = 1/(1 + rA/400);
+  y = y(+1) - (1/tau)*(R - pi(+1) - z(+1)) + g - g(+1);
+  pi = beta*pi(+1) + kappa*(y - g);
+  R = rho_R*R(-1) + (1 - rho_R)*(psi1*pi + psi2*(y - g)) + s_ep*ep/100;
+  g = rho_g*g(-1) + s_eg*eg/100;
+  z = rho_z*z(-1) + s_ez*ez/100;
+  YGR = gammaQ + 100*(y - y(-1) + z);
+  INFL = piA + 400*pi;
+  INT = piA + rA + 4*gammaQ + 400*R;
+end;
+
+steady_state_model;
+  y = 0;
+  pi = 0;
+  R = 0;
+  g = 0;
+  z = 0;
+  YGR = gammaQ;
+  INFL = piA;
+  INT = piA + rA + 4*gammaQ;
+end;
+
+shocks;
+  var ep; stderr 1.0;
+  var eg; stderr 1.0;
+  var ez; stderr 1.0;
+  // calibrated measurement errors
+  var YGR;  stderr (0.20*0.579923);
+  var INFL; stderr (0.20*1.470832);
+  var INT;  stderr(0.20*2.237937);
+end;
+
+estimated_params;
+rA, gamma_pdf, 0.5, 0.2;
+piA, gamma_pdf, 7, 2;
+gammaQ, normal_pdf, 0.4, 0.2;
+tau, gamma_pdf, 2, 0.5;
+kappa, beta_pdf, 0.5, 0.2;
+psi1, gamma_pdf, 1.5, 0.25;
+psi2, gamma_pdf, 0.5, 0.25;
+rho_R, beta_pdf, 0.5, 0.2;
+rho_g, beta_pdf, 0.5, 0.2;
+rho_z, beta_pdf, 0.5, 0.2;
+s_ep, inv_gamma_pdf, 0.5013256549262002, 0.06867258771281654;
+s_eg, inv_gamma_pdf, 1.2533141373155003, 0.4292036732051032;
+s_ez, inv_gamma_pdf, 0.6266570686577502, 0.1073009183012758;
+end;
+
+estimated_params_init;
+  tau, 2.263895645223946;
+  kappa, 0.9;
+  psi1, 1.932357745633903;
+  psi2, 0.465714640873466;
+  rA, 0.333453026636814;
+  piA, 3.427852381637088;
+  gammaQ, 0.624767199308368;
+  rho_R, 0.764476492352786;
+  rho_g, 0.990359718181570;
+  rho_z, 0.917302546854851;
+  s_ep,   0.211571288023659;
+  s_eg, 0.632177304244290;
+  s_ez, 0.192160214260411;
+end;
+
+varobs YGR INFL INT;
+
+estimation(datafile='dsge1_data.csv', mh_replic=10000, mh_jscale=0.2);
diff --git a/test/models/herbst_schorfheide/nk_est_a.mod b/test/models/herbst_schorfheide/nk_est_a.mod
new file mode 100644
index 00000000..94359a3a
--- /dev/null
+++ b/test/models/herbst_schorfheide/nk_est_a.mod
@@ -0,0 +1,89 @@
+var y R pi z g YGR INFL INT;
+varexo eg ep ez;
+
+parameters tau rho_R psi1 psi2 rho_g rho_z kappa  piA gammaQ rA s_ep s_eg s_ez; 
+
+rA = 0.42;
+piA = 3.3;
+gammaQ = 0.52;
+tau = 2.83;
+kappa = 0.78;
+psi1 = 1.8;
+psi2 = 0.63;
+rho_R = 0.77;
+rho_g = 0.98;
+rho_z = 0.88;
+s_ep = 0.22;
+s_eg = 0.72;
+s_ez = 0.31;
+
+
+model;
+#  beta = 1/(1 + rA/400);
+  y = y(+1) - (1/tau)*(R - pi(+1) - z(+1)) + g - g(+1);
+  pi = beta*pi(+1) + kappa*(y - g);
+  R = rho_R*R(-1) + (1 - rho_R)*(psi1*pi + psi2*(y - g)) + s_ep*ep/100;
+  g = rho_g*g(-1) + s_eg*eg/100;
+  z = rho_z*z(-1) + s_ez*ez/100;
+  YGR = gammaQ + 100*(y - y(-1) + z);
+  INFL = piA + 400*pi;
+  INT = piA + rA + 4*gammaQ + 400*R;
+end;
+
+steady_state_model;
+  y = 0;
+  pi = 0;
+  R = 0;
+  g = 0;
+  z = 0;
+  YGR = gammaQ;
+  INFL = piA;
+  INT = piA + rA + 4*gammaQ;
+end;
+
+shocks;
+  var ep; stderr 1.0;
+  var eg; stderr 1.0;
+  var ez; stderr 1.0;
+  // calibrated measurement errors
+  var YGR;  stderr (0.20*0.579923);
+  var INFL; stderr (0.20*1.470832);
+  var INT;  stderr(0.20*2.237937);
+end;
+
+estimated_params;
+rA, gamma_pdf, 0.5, 0.2;
+piA, gamma_pdf, 7, 2;
+gammaQ, normal_pdf, 0.4, 0.2;
+tau, gamma_pdf, 2, 0.5;
+kappa, beta_pdf, 0.5, 0.2;
+psi1, gamma_pdf, 1.5, 0.25;
+psi2, gamma_pdf, 0.5, 0.25;
+rho_R, beta_pdf, 0.5, 0.2;
+rho_g, beta_pdf, 0.5, 0.2;
+rho_z, beta_pdf, 0.5, 0.2;
+s_ep, inv_gamma_pdf, 0.5013256549262002, 0.06867258771281654;
+s_eg, inv_gamma_pdf, 1.2533141373155003, 0.4292036732051032;
+s_ez, inv_gamma_pdf, 0.6266570686577502, 0.1073009183012758;
+end;
+
+estimated_params_init;
+  tau, 2.263895645223946;
+  kappa, 0.9;
+  psi1, 1.932357745633903;
+  psi2, 0.465714640873466;
+  rA, 0.333453026636814;
+  piA, 3.427852381637088;
+  gammaQ, 0.624767199308368;
+  rho_R, 0.764476492352786;
+  rho_g, 0.990359718181570;
+  rho_z, 0.917302546854851;
+  s_ep,   0.211571288023659;
+  s_eg, 0.632177304244290;
+  s_ez, 0.192160214260411;
+end;
+
+varobs YGR INFL INT;
+
+//estimation(datafile='dsge1_data.csv', mh_replic=0);
+rwmh_compute!(datafile="dsge1_data.csv", mcmc_replic=100000, mcmc_jscale=0.03, mcmc_chains=2, transformed_parameters=false, initial_values=context.work.estimated_parameters.initialvalue);
diff --git a/test/models/herbst_schorfheide/nk_recfor.mod b/test/models/herbst_schorfheide/nk_recfor.mod
new file mode 100644
index 00000000..8fffa4cc
--- /dev/null
+++ b/test/models/herbst_schorfheide/nk_recfor.mod
@@ -0,0 +1,60 @@
+var y R pi z g YGR INFL INT;
+varexo eg ep ez;
+
+parameters tau rho_R psi1 psi2 rho_g rho_z kappa  piA gammaQ rA s_ep s_eg s_ez; 
+
+rA = 0.42;
+piA = 3.3;
+gammaQ = 0.52;
+tau = 2.83;
+kappa = 0.78;
+psi1 = 1.8;
+psi2 = 0.63;
+rho_R = 0.77;
+rho_g = 0.98;
+rho_z = 0.88;
+s_ep = 0.22;
+s_eg = 0.72;
+s_ez = 0.31;
+
+
+model;
+#  beta = 1/(1 + rA/400);
+  y = y(+1) - (1/tau)*(R - pi(+1) - z(+1)) + g - g(+1);
+  pi = beta*pi(+1) + kappa*(y - g);
+  R = rho_R*R(-1) + (1 - rho_R)*(psi1*pi + psi2*(y - g)) + s_ep*ep/100;
+  g = rho_g*g(-1) + s_eg*eg/100;
+  z = rho_z*z(-1) + s_ez*ez/100;
+  YGR = gammaQ + 100*(y - y(-1) + z);
+  INFL = piA + 400*pi;
+  INT = piA + rA + 4*gammaQ + 400*R;
+end;
+
+steady_state_model;
+  y = 0;
+  pi = 0;
+  R = 0;
+  g = 0;
+  z = 0;
+  YGR = gammaQ;
+  INFL = piA;
+  INT = piA + rA + 4*gammaQ;
+end;
+
+varobs YGR INFL INT;
+
+shocks;
+  var ep; stderr 1.0;
+  var eg; stderr 1.0;
+  var ez; stderr 1.0;
+  // calibrated measurement errors
+  var YGR;  stderr (0.20*0.579923);
+  var INFL; stderr (0.20*1.470832);
+  var INT;  stderr(0.20*2.237937);
+end;
+
+stoch_simul(order=1, irf=0);
+
+recursive_forecasting!(periods=12, first_period=78, last_period=80, datafile="dsge1_data.csv");
+
+
diff --git a/test/models/herbst_schorfheide/nk_sol.mod b/test/models/herbst_schorfheide/nk_sol.mod
new file mode 100644
index 00000000..3aa43d1d
--- /dev/null
+++ b/test/models/herbst_schorfheide/nk_sol.mod
@@ -0,0 +1,50 @@
+var y R pi z g YGR INFL INT;
+varexo eg ep ez;
+
+parameters tau rho_R psi1 psi2 rho_g rho_z kappa  piA gammaQ rA s_ep s_eg s_ez; 
+
+rA = 0.42;
+piA = 3.3;
+gammaQ = 0.52;
+tau = 2.83;
+kappa = 0.78;
+psi1 = 1.8;
+psi2 = 0.63;
+rho_R = 0.77;
+rho_g = 0.98;
+rho_z = 0.88;
+s_ep = 0.22;
+s_eg = 0.72;
+s_ez = 0.31;
+
+
+model;
+#  beta = 1/(1 + rA/400);
+  y = y(+1) - (1/tau)*(R - pi(+1) - z(+1)) + g - g(+1);
+  pi = beta*pi(+1) + kappa*(y - g);
+  R = rho_R*R(-1) + (1 - rho_R)*(psi1*pi + psi2*(y - g)) + s_ep*ep/100;
+  g = rho_g*g(-1) + s_eg*eg/100;
+  z = rho_z*z(-1) + s_ez*ez/100;
+  YGR = gammaQ + 100*(y - y(-1) + z);
+  INFL = piA + 400*pi;
+  INT = piA + rA + 4*gammaQ + 400*R;
+end;
+
+steady_state_model;
+  y = 0;
+  pi = 0;
+  R = 0;
+  g = 0;
+  z = 0;
+  YGR = gammaQ;
+  INFL = piA;
+  INT = piA + rA + 4*gammaQ;
+end;
+
+shocks;
+  var ep; stderr 1.0;
+  var eg; stderr 1.0;
+  var ez; stderr 1.0;
+end;
+
+stoch_simul(order=1, irf=0);