diff --git a/matlab/dynare_estimation_1.m b/matlab/dynare_estimation_1.m
index c084a211c68b0b81c1b5de69a00adc04f848fa3c..4f4fa37f330c1dba55640ab1078bcae4094f8bde 100644
--- a/matlab/dynare_estimation_1.m
+++ b/matlab/dynare_estimation_1.m
@@ -106,6 +106,20 @@ if ~isempty(options_.mode_file) && ~options_.mh_posterior_mode_estimation
load(options_.mode_file);
end
+%% load optimal_mh_scale parameter if previous run was with
+%% mode_compute=6
+mh_scale_fname = [M_.fname '_optimal_mh_scale_parameter.mat'];
+if exist(mh_scale_fname)
+ if options_.mode_compute == 0
+ tmp = load(mh_scale_fname,'Scale');
+ bayestopt_.mh_jscale = tmp.Scale;
+ clear tmp;
+ else
+ % remove the file if mode_compute ~= 0
+ delete('mh_scale_fname')
+ end
+end
+
if ~isempty(estim_params_)
set_parameters(xparam1);
end
@@ -170,11 +184,6 @@ if isequal(options_.mode_compute,0) && isempty(options_.mode_file) && options_.m
return
end
-if isequal(options_.mode_compute,6)
- % Erase previously computed optimal mh scale parameter.
- delete([M_.fname '_optimal_mh_scale_parameter.mat'])
-end
-
%% Estimation of the posterior mode or likelihood mode
if ~isequal(options_.mode_compute,0) && ~options_.mh_posterior_mode_estimation
diff --git a/matlab/random_walk_metropolis_hastings_core.m b/matlab/random_walk_metropolis_hastings_core.m
index b5ad698ecb0306a60d2cb0e4257eae0df495c358..69062f66c735746eb316f0d69b9b9f360e057342 100644
--- a/matlab/random_walk_metropolis_hastings_core.m
+++ b/matlab/random_walk_metropolis_hastings_core.m
@@ -112,23 +112,12 @@ if strcmpi(ProposalFun,'rand_multivariate_normal')
elseif strcmpi(ProposalFun,'rand_multivariate_student')
n = options_.student_degrees_of_freedom;
end
-% load([MhDirectoryName '/' ModelName '_mh_history.mat'],'record');
+
%%%%
%%%% NOW i run the (nblck-fblck+1) metropolis-hastings chains
%%%%
-
-if any(isnan(bayestopt_.jscale)) || (options_.mode_compute==0 && ~isempty(options_.mode_file))
- if exist([ModelName '_optimal_mh_scale_parameter.mat'])% This file is created by mode_compute=6.
- load([ModelName '_optimal_mh_scale_parameter'])
- proposal_covariance_Cholesky_decomposition = d*Scale;
- else
- error('mh:: Something is wrong. I can''t figure out the value of the scale parameter.')
- end
-else
- proposal_covariance_Cholesky_decomposition = d*diag(bayestopt_.jscale);
-end
-
+proposal_covariance_Cholesky_decomposition = d*diag(bayestopt_.jscale);
jloop=0;
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 9d3cd67296bf630ba7cdef44bbd1243f544722d4..9e3c6d614f3cf1a637c8ee40d39e532e13487776 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -139,8 +139,11 @@ MODFILES = \
deterministic_simulations/rbc_det3.mod \
deterministic_simulations/rbc_det4.mod \
deterministic_simulations/rbc_det5.mod \
- fs2000/fs2000_sd.mod
-
+ fs2000/fs2000_sd.mod \
+ estimation/fs2000_mc4.mod \
+ estimation/fs2000_mc4_mf.mod \
+ estimation/fs2000_mc6.mod \
+ estimation/fs2000_mc6_mf.mod
EXTRA_DIST = \
$(MODFILES) \
homotopy/common.mod \
@@ -191,6 +194,7 @@ EXTRA_DIST = \
kalman_filter_smoother/fsdat_simul.m \
kalman_filter_smoother/fs2000a_steadystate.m \
identification/kim/kim2_steadystate.m \
+ estimation/fsdat_simul.m \
ep/mean_preserving_spread.m
TARGETS =
diff --git a/tests/estimation/fs2000.mod b/tests/estimation/fs2000.mod
new file mode 100644
index 0000000000000000000000000000000000000000..981dd35f61a47e0be29b6333f06cf73609791971
--- /dev/null
+++ b/tests/estimation/fs2000.mod
@@ -0,0 +1,75 @@
+// See fs2000.mod in the examples/ directory for details on the model
+
+var m P c e W R k d n l gy_obs gp_obs y dA;
+varexo e_a e_m;
+
+parameters alp bet gam mst rho psi del;
+
+alp = 0.33;
+bet = 0.99;
+gam = 0.003;
+mst = 1.011;
+rho = 0.7;
+psi = 0.787;
+del = 0.02;
+
+model;
+dA = exp(gam+e_a);
+log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
+-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
+W = l/n;
+-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
+R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
+1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
+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));
+gy_obs = dA*y/y(-1);
+gp_obs = (P/P(-1))*m(-1)/dA;
+end;
+
+initval;
+k = 6;
+m = mst;
+P = 2.25;
+c = 0.45;
+e = 1;
+W = 4;
+R = 1.02;
+d = 0.85;
+n = 0.19;
+l = 0.86;
+y = 0.6;
+gy_obs = exp(gam);
+gp_obs = exp(-gam);
+dA = exp(gam);
+end;
+
+shocks;
+var e_a; stderr 0.014;
+var e_m; stderr 0.005;
+end;
+
+steady;
+
+check;
+
+estimated_params;
+alp, beta_pdf, 0.356, 0.02;
+bet, beta_pdf, 0.993, 0.002;
+gam, normal_pdf, 0.0085, 0.003;
+mst, normal_pdf, 1.0002, 0.007;
+rho, beta_pdf, 0.129, 0.223;
+psi, beta_pdf, 0.65, 0.05;
+del, beta_pdf, 0.01, 0.005;
+stderr e_a, inv_gamma_pdf, 0.035449, inf;
+stderr e_m, inv_gamma_pdf, 0.008862, inf;
+end;
+
+varobs gp_obs gy_obs;
+
+options_.solve_tolf = 1e-12;
+
+estimation(order=1,datafile=fsdat_simul,nobs=192,loglinear,mh_replic=2000,mh_nblocks=2,mh_jscale=0.8);
diff --git a/tests/estimation/fs2000_mc4.mod b/tests/estimation/fs2000_mc4.mod
new file mode 100644
index 0000000000000000000000000000000000000000..8b9ebb4a1ffaf2a2b04a2375618553c68b4d0ea7
--- /dev/null
+++ b/tests/estimation/fs2000_mc4.mod
@@ -0,0 +1,73 @@
+// test mode_compute default (=4)
+// See fs2000.mod in the examples/ directory for details on the model
+
+var m P c e W R k d n l gy_obs gp_obs y dA;
+varexo e_a e_m;
+
+parameters alp bet gam mst rho psi del;
+
+alp = 0.33;
+bet = 0.99;
+gam = 0.003;
+mst = 1.011;
+rho = 0.7;
+psi = 0.787;
+del = 0.02;
+
+model;
+dA = exp(gam+e_a);
+log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
+-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
+W = l/n;
+-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
+R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
+1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
+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));
+gy_obs = dA*y/y(-1);
+gp_obs = (P/P(-1))*m(-1)/dA;
+end;
+
+steady_state_model;
+ dA = exp(gam);
+ gst = 1/dA;
+ m = mst;
+
+ khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1));
+ xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1);
+ nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp );
+ n = xist/(nust+xist);
+ P = xist + nust;
+ k = khst*n;
+
+ l = psi*mst*n/( (1-psi)*(1-n) );
+ c = mst/P;
+ d = l - mst + 1;
+ y = k^alp*n^(1-alp)*gst^alp;
+ R = mst/bet;
+
+ W = l/n;
+ e = 1;
+
+ gp_obs = m/dA;
+ gy_obs = dA;
+end;
+
+estimated_params;
+alp, beta_pdf, 0.356, 0.02;
+bet, beta_pdf, 0.993, 0.002;
+gam, normal_pdf, 0.0085, 0.003;
+mst, normal_pdf, 1.0002, 0.007;
+rho, beta_pdf, 0.129, 0.223;
+psi, beta_pdf, 0.65, 0.05;
+del, beta_pdf, 0.01, 0.005;
+stderr e_a, inv_gamma_pdf, 0.035449, inf;
+stderr e_m, inv_gamma_pdf, 0.008862, inf;
+end;
+
+varobs gp_obs gy_obs;
+
+estimation(order=1,datafile=fsdat_simul,nobs=192,loglinear,mh_replic=0);
diff --git a/tests/estimation/fs2000_mc4_mf.mod b/tests/estimation/fs2000_mc4_mf.mod
new file mode 100644
index 0000000000000000000000000000000000000000..6c0a9428e5a9f06f10dcaeb0beb70fd222064b0f
--- /dev/null
+++ b/tests/estimation/fs2000_mc4_mf.mod
@@ -0,0 +1,66 @@
+// test mode_compute=0 with mode_file computed with default (mode_compute=4)
+// See fs2000.mod in the examples/ directory for details on the model
+
+var m P c e W R k d n l gy_obs gp_obs y dA;
+varexo e_a e_m;
+
+parameters alp bet gam mst rho psi del;
+
+
+model;
+dA = exp(gam+e_a);
+log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
+-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
+W = l/n;
+-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
+R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
+1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
+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));
+gy_obs = dA*y/y(-1);
+gp_obs = (P/P(-1))*m(-1)/dA;
+end;
+
+steady_state_model;
+ dA = exp(gam);
+ gst = 1/dA;
+ m = mst;
+
+ khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1));
+ xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1);
+ nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp );
+ n = xist/(nust+xist);
+ P = xist + nust;
+ k = khst*n;
+
+ l = psi*mst*n/( (1-psi)*(1-n) );
+ c = mst/P;
+ d = l - mst + 1;
+ y = k^alp*n^(1-alp)*gst^alp;
+ R = mst/bet;
+
+ W = l/n;
+ e = 1;
+
+ gp_obs = m/dA;
+ gy_obs = dA;
+end;
+
+estimated_params;
+alp, beta_pdf, 0.356, 0.02;
+bet, beta_pdf, 0.993, 0.002;
+gam, normal_pdf, 0.0085, 0.003;
+mst, normal_pdf, 1.0002, 0.007;
+rho, beta_pdf, 0.129, 0.223;
+psi, beta_pdf, 0.65, 0.05;
+del, beta_pdf, 0.01, 0.005;
+stderr e_a, inv_gamma_pdf, 0.035449, inf;
+stderr e_m, inv_gamma_pdf, 0.008862, inf;
+end;
+
+varobs gp_obs gy_obs;
+
+estimation(order=1,datafile=fsdat_simul,nobs=192,loglinear,mode_compute=0,mode_file=fs2000_mc4_mode,mh_replic=2000,mh_jscale=0.8);
diff --git a/tests/estimation/fs2000_mc6.mod b/tests/estimation/fs2000_mc6.mod
new file mode 100644
index 0000000000000000000000000000000000000000..0819affdfc2fd27fdf34bc5b7dcbc1b62de340ea
--- /dev/null
+++ b/tests/estimation/fs2000_mc6.mod
@@ -0,0 +1,73 @@
+// test mode_compute=6
+// See fs2000.mod in the examples/ directory for details on the model
+
+var m P c e W R k d n l gy_obs gp_obs y dA;
+varexo e_a e_m;
+
+parameters alp bet gam mst rho psi del;
+
+alp = 0.33;
+bet = 0.99;
+gam = 0.003;
+mst = 1.011;
+rho = 0.7;
+psi = 0.787;
+del = 0.02;
+
+model;
+dA = exp(gam+e_a);
+log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
+-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
+W = l/n;
+-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
+R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
+1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
+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));
+gy_obs = dA*y/y(-1);
+gp_obs = (P/P(-1))*m(-1)/dA;
+end;
+
+steady_state_model;
+ dA = exp(gam);
+ gst = 1/dA;
+ m = mst;
+
+ khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1));
+ xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1);
+ nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp );
+ n = xist/(nust+xist);
+ P = xist + nust;
+ k = khst*n;
+
+ l = psi*mst*n/( (1-psi)*(1-n) );
+ c = mst/P;
+ d = l - mst + 1;
+ y = k^alp*n^(1-alp)*gst^alp;
+ R = mst/bet;
+
+ W = l/n;
+ e = 1;
+
+ gp_obs = m/dA;
+ gy_obs = dA;
+end;
+
+estimated_params;
+alp, beta_pdf, 0.356, 0.02;
+bet, beta_pdf, 0.993, 0.002;
+gam, normal_pdf, 0.0085, 0.003;
+mst, normal_pdf, 1.0002, 0.007;
+rho, beta_pdf, 0.129, 0.223;
+psi, beta_pdf, 0.65, 0.05;
+del, beta_pdf, 0.01, 0.005;
+stderr e_a, inv_gamma_pdf, 0.035449, inf;
+stderr e_m, inv_gamma_pdf, 0.008862, inf;
+end;
+
+varobs gp_obs gy_obs;
+
+estimation(order=1,datafile=fsdat_simul,nobs=192,loglinear,mh_replic=0,mode_compute=6);
diff --git a/tests/estimation/fs2000_mc6_mf.mod b/tests/estimation/fs2000_mc6_mf.mod
new file mode 100644
index 0000000000000000000000000000000000000000..577abb4a6db16fdd801e89b89b329c5a9a5befa2
--- /dev/null
+++ b/tests/estimation/fs2000_mc6_mf.mod
@@ -0,0 +1,66 @@
+// test mode_compute=0 with mode_file computed with mode_compute=6
+// See fs2000.mod in the examples/ directory for details on the model
+
+var m P c e W R k d n l gy_obs gp_obs y dA;
+varexo e_a e_m;
+
+parameters alp bet gam mst rho psi del;
+
+
+model;
+dA = exp(gam+e_a);
+log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
+-P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c(+2)*P(+2)*m(+1))=0;
+W = l/n;
+-(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
+R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
+1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
+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));
+gy_obs = dA*y/y(-1);
+gp_obs = (P/P(-1))*m(-1)/dA;
+end;
+
+steady_state_model;
+ dA = exp(gam);
+ gst = 1/dA;
+ m = mst;
+
+ khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1));
+ xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-1);
+ nust = psi*mst^2/( (1-alp)*(1-psi)*bet*gst^alp*khst^alp );
+ n = xist/(nust+xist);
+ P = xist + nust;
+ k = khst*n;
+
+ l = psi*mst*n/( (1-psi)*(1-n) );
+ c = mst/P;
+ d = l - mst + 1;
+ y = k^alp*n^(1-alp)*gst^alp;
+ R = mst/bet;
+
+ W = l/n;
+ e = 1;
+
+ gp_obs = m/dA;
+ gy_obs = dA;
+end;
+
+estimated_params;
+alp, beta_pdf, 0.356, 0.02;
+bet, beta_pdf, 0.993, 0.002;
+gam, normal_pdf, 0.0085, 0.003;
+mst, normal_pdf, 1.0002, 0.007;
+rho, beta_pdf, 0.129, 0.223;
+psi, beta_pdf, 0.65, 0.05;
+del, beta_pdf, 0.01, 0.005;
+stderr e_a, inv_gamma_pdf, 0.035449, inf;
+stderr e_m, inv_gamma_pdf, 0.008862, inf;
+end;
+
+varobs gp_obs gy_obs;
+
+estimation(order=1,datafile=fsdat_simul,nobs=192,loglinear,mode_compute=0,mode_file=fs2000_mc6_mode,mh_replic=2000,mh_jscale=0.8);