fs2000.mod: provide actual replication

Closes https://git.dynare.org/Dynare/dynare/-/issues/1905

examples/fs2000.mod

@@ -3,20 +3,20 @@
  * in the paper) described in Frank Schorfheide (2000): "Loss function-based 
  * evaluation of DSGE models", Journal of Applied Econometrics, 15(6), 645-670.
  *
- * The data are in file "fsdat_simul.m", and have been artificially generated.
- * They are therefore different from the original dataset used by Schorfheide.
+ * The data are taken from the replication package at 
+ * http://dx.doi.org/10.15456/jae.2022314.0708799949
  *
  * The prior distribution follows the one originally specified in Schorfheide's
- * paper, except for parameter rho. In the paper, the elicited beta prior for rho
+ * paper. Note that the elicited beta prior for rho in the paper
  * implies an asymptote and corresponding prior mode at 0. It is generally
- * recommended to avoid this extreme type of prior. Some optimizers, for instance
- * mode_compute=12 (Mathworks' particleswarm algorithm) may find a posterior mode
- * with rho equal to zero. We lowered the value of the prior standard deviation
- * (changing .223 to .100) to remove the asymptote.
+ * recommended to avoid this extreme type of prior. 
+ *
+ * Because the data are already logged and we use the loglinear option to conduct 
+ * a full log-linearization, we need to use the logdata option.
  *
  * The equations are taken from J. Nason and T. Cogley (1994): "Testing the
  * implications of long-run neutrality for monetary business cycle models",
- * Journal of Applied Econometrics, 9, S37-S70.
+ * Journal of Applied Econometrics, 9, S37-S70, NC in the following.
  * Note that there is an initial minus sign missing in equation (A1), p. S63.
  *
  * This implementation was originally written by Michel Juillard. Please note that the
@@ -25,7 +25,7 @@
  */
 
 /*
- * Copyright © 2004-2017 Dynare Team
+ * Copyright © 2004-2023 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -43,33 +43,71 @@
  * along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
  */
 
-var m P c e W R k d n l gy_obs gp_obs y dA;
-varexo e_a e_m;
+var m       ${m}$           (long_name='money growth')
+    P       ${P}$           (long_name='Price level')
+    c       ${c}$           (long_name='consumption')
+    e       ${e}$           (long_name='capital stock')
+    W       ${W}$           (long_name='Wage rate')
+    R       ${R}$           (long_name='interest rate')
+    k       ${k}$           (long_name='capital stock')
+    d       ${d}$           (long_name='dividends')
+    n       ${n}$           (long_name='labor')
+    l       ${l}$           (long_name='loans')
+    gy_obs  ${\Delta \ln GDP}$  (long_name='detrended capital stock')
+    gp_obs  ${\Delta \ln P}$    (long_name='detrended capital stock')
+    y       ${y}$           (long_name='detrended output')
+    dA      ${\Delta A}$    (long_name='TFP growth')
+    ;
+varexo e_a  ${\epsilon_A}$      (long_name='TFP shock')
+    e_m     ${\epsilon_M}$      (long_name='Money growth shock')
+    ;
 
-parameters alp bet gam mst rho psi del;
+parameters alp  ${\alpha}$      (long_name='capital share')
+    bet         ${\beta}$       (long_name='discount factor')
+    gam         ${\gamma}$      (long_name='long-run TFP growth')
+    mst         ${m^*}$         (long_name='long-run money growth')
+    rho         ${\rho}$        (long_name='autocorrelation money growth')
+    phi         ${\phi}$        (long_name='labor weight in consumption')
+    del         ${\delta}$      (long_name='depreciation rate')
+    ;
 
+% roughly picked values to allow simulating the model before estimation
 alp = 0.33;
 bet = 0.99;
 gam = 0.003;
 mst = 1.011;
 rho = 0.7;
-psi = 0.787;
+phi = 0.787;
 del = 0.02;
 
 model;
+[name='NC before eq. (1), TFP growth equation']
 dA = exp(gam+e_a);
+[name='NC eq. (2), money growth rate']
 log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
+[name='NC eq. (A1), Euler equation']
 -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;
+[name='NC below eq. (A1), firm borrowing constraint']
 W = l/n;
--(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
+[name='NC eq. (A2), intratemporal labour market condition']
+-(phi/(1-phi))*(c*P/(1-n))+l/n = 0;
+[name='NC below eq. (A2), credit market clearing']
 R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
+[name='NC eq. (A3), credit market optimality']
 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;
+[name='NC eq. (18), aggregate resource constraint']
 c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
+[name='NC eq. (19), money market condition']
 P*c = m;
+[name='NC eq. (20), credit market equilibrium condition']
 m-1+d = l;
+[name='Definition TFP shock']
 e = exp(e_a);
+[name='Implied by NC eq. (18), production function']
 y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
+[name='Observation equation GDP growth']
 gy_obs = dA*y/y(-1);
+[name='Observation equation price level']
 gp_obs = (P/P(-1))*m(-1)/dA;
 end;
 
@@ -84,12 +122,12 @@ steady_state_model;
   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 );
+  nust = phi*mst^2/( (1-alp)*(1-phi)*bet*gst^alp*khst^alp );
   n  = xist/(nust+xist);
   P  = xist + nust;
   k  = khst*n;
 
-  l  = psi*mst*n/( (1-psi)*(1-n) );
+  l  = phi*mst*n/( (1-phi)*(1-n) );
   c  = mst/P;
   d  = l - mst + 1;
   y  = k^alp*n^(1-alp)*gst^alp;
@@ -104,17 +142,18 @@ steady_state_model;
   gy_obs = dA;
 end;
 
-steady;
 
+steady;
 check;
 
+% Table 1 of Schorfheide (2000)
 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.100;
-psi, beta_pdf, 0.65, 0.05;
+rho, beta_pdf, 0.129, 0.223;
+phi, 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;
@@ -122,14 +161,8 @@ end;
 
 varobs gp_obs gy_obs;
 
-estimation(order=1, datafile=fsdat_simul, nobs=192, loglinear, mh_replic=2000, mh_nblocks=2, mh_jscale=0.8, mode_check);
+estimation(order=1, datafile=fs2000_data, loglinear,logdata, mh_replic=2000, mh_nblocks=2, mh_jscale=0.8, mode_check);
 
-
-/*
- * The following lines were used to generate the data file. If you want to
- * generate another random data file, comment the "estimation" line and uncomment
- * the following lines.
- */
-
-//stoch_simul(periods=200, order=1);
-//datatomfile('fsdat_simul', {'gy_obs', 'gp_obs'});
+%uncomment the following lines to generate LaTeX-code of the model equations
+%write_latex_original_model(write_equation_tags);
+%collect_latex_files;

examples/fs2000_data.m

@@ -0,0 +1,215 @@
+%This file is a direct Matlab implementation of the loaddata.g and data.prn files
+%of Schorfheide, Frank (2000): Loss function-based evaluation of DSGE models 
+%(replication data). Version: 1. Journal of Applied Econometrics. Dataset. 
+%http://dx.doi.org/10.15456/jae.2022314.0708799949
+% Copyright: 2000-2022 Frank Schorfheide
+% Copyright: 2023 Dynare Team
+
+% License: CC BY 4.0
+% (https://creativecommons.org/licenses/by/4.0/legalcode)
+
+% Time series, extracted 05/04/00
+% columms are quarterly data from 1949:IV to 1997:IV
+% 1: GDPD = GROSS DOMESTIC PRODUCT:IMPLICIT PRICE DEFLATOR (INDEX,92=100)(T7.1)
+% 2: GDPQ = GROSS DOMESTIC PRODUCT
+% 3: GPOP = POPULATION, NIPA basis (THOUS.,NSA)
+
+data_q=[18.02 1474.5 150.2
+17.94 1538.2 150.9
+18.01 1584.5 151.4
+18.42 1644.1 152
+18.73 1678.6 152.7
+19.46 1693.1 153.3
+19.55 1724   153.9
+19.56 1758.2 154.7
+19.79 1760.6 155.4
+19.77 1779.2 156
+19.82 1778.8 156.6
+20.03 1790.9 157.3
+20.12 1846   158
+20.1  1882.6 158.6
+20.14 1897.3 159.2
+20.22 1887.4 160
+20.27 1858.2 160.7
+20.34 1849.9 161.4
+20.39 1848.5 162
+20.42 1868.9 162.8
+20.47 1905.6 163.6
+20.56 1959.6 164.3
+20.62 1994.4 164.9
+20.78 2020.1 165.7
+21    2030.5 166.5
+21.2  2023.6 167.2
+21.33 2037.7 167.9
+21.62 2033.4 168.7
+21.71 2066.2 169.5
+22.01 2077.5 170.2
+22.15 2071.9 170.9
+22.27 2094   171.7
+22.29 2070.8 172.5
+22.56 2012.6 173.1
+22.64 2024.7 173.8
+22.77 2072.3 174.5
+22.88 2120.6 175.3
+22.92 2165   176.045
+22.91 2223.3  176.727
+22.94 2221.4  177.481
+23.03 2230.95 178.268
+23.13 2279.22 179.694
+23.22 2265.48 180.335
+23.32 2268.29 181.094
+23.4  2238.57 181.915
+23.45 2251.68 182.634
+23.51 2292.02 183.337
+23.56 2332.61 184.103
+23.63 2381.01 184.894
+23.75 2422.59 185.553
+23.81 2448.01 186.203
+23.87 2471.86 186.926
+23.94 2476.67 187.68
+24    2508.7  188.299
+24.07 2538.05 188.906
+24.12 2586.26 189.631
+24.29 2604.62 190.362
+24.35 2666.69 190.954
+24.41 2697.54 191.56
+24.52 2729.63 192.256
+24.64 2739.75 192.938
+24.77 2808.88 193.467
+24.88 2846.34 193.994
+25.01 2898.79 194.647
+25.17 2970.48 195.279
+25.32 3042.35 195.763
+25.53 3055.53 196.277
+25.79 3076.51 196.877
+26.02 3102.36 197.481
+26.14 3127.15 197.967
+26.31 3129.53 198.455
+26.6  3154.19 199.012
+26.9  3177.98 199.572
+27.21 3236.18 199.995
+27.49 3292.07 200.452
+27.75 3316.11 200.997
+28.12 3331.22 201.538
+28.39 3381.86 201.955
+28.73 3390.23 202.419
+29.14 3409.65 202.986
+29.51 3392.6  203.584
+29.94 3386.49 204.086
+30.36 3391.61 204.721
+30.61 3422.95 205.419
+31.02 3389.36 206.13
+31.5  3481.4  206.763
+31.93 3500.95 207.362
+32.27 3523.8  208
+32.54 3533.79 208.642
+33.02 3604.73 209.142
+33.2  3687.9  209.637
+33.49 3726.18 210.181
+33.95 3790.44 210.737
+34.36 3892.22 211.192
+34.94 3919.01 211.663
+35.61 3907.08 212.191
+36.29 3947.11 212.708
+37.01 3908.15 213.144
+37.79 3922.57 213.602
+38.96 3879.98 214.147
+40.13 3854.13 214.7
+41.05 3800.93 215.135
+41.66 3835.21 215.652
+42.41 3907.02 216.289
+43.19 3952.48 216.848
+43.69 4044.59 217.314
+44.15 4072.19 217.776
+44.77 4088.49 218.338
+45.57 4126.39 218.917
+46.32 4176.28 219.427
+47.07 4260.08 219.956
+47.66 4329.46 220.573
+48.63 4328.33 221.201
+49.42 4345.51 221.719
+50.41 4510.73 222.281
+51.27 4552.14 222.933
+52.35 4603.65 223.583
+53.51 4605.65 224.152
+54.65 4615.64 224.737
+55.82 4644.93 225.418
+56.92 4656.23 226.117
+58.18 4678.96 226.754
+59.55 4566.62 227.389
+61.01 4562.25 228.07
+62.59 4651.86 228.689
+64.15 4739.16 229.155
+65.37 4696.82 229.674
+66.65 4753.02 230.301
+67.87 4693.76 230.903
+68.86 4615.89 231.395
+69.72 4634.88 231.906
+70.66 4612.08 232.498
+71.44 4618.26 233.074
+72.08 4662.97 233.546
+72.83 4763.57 234.028
+73.48 4849    234.603
+74.19 4939.23 235.153
+75.02 5053.56 235.605
+75.58 5132.87 236.082
+76.25 5170.34 236.657
+76.81 5203.68 237.232
+77.63 5257.26 237.673
+78.25 5283.73 238.176
+78.76 5359.6  238.789
+79.45 5393.57 239.387
+79.81 5460.83 239.861
+80.22 5466.95 240.368
+80.84 5496.29 240.962
+81.45 5526.77 241.539
+82.09 5561.8  242.009
+82.68 5618    242.52
+83.33 5667.39 243.12
+84.09 5750.57 243.721
+84.67 5785.29 244.208
+85.56 5844.05 244.716
+86.66 5878.7  245.354
+87.44 5952.83 245.966
+88.45 6010.96 246.46
+89.39 6055.61 247.017
+90.13 6087.96 247.698
+90.88 6093.51 248.374
+92    6152.59 248.928
+93.18 6171.57 249.564
+94.14 6142.1  250.299
+95.11 6078.96 251.031
+96.27 6047.49 251.65
+97    6074.66 252.295
+97.7  6090.14 253.033
+98.31 6105.25 253.743
+99.13 6175.69 254.338
+99.79 6214.22 255.032
+100.17 6260.74 255.815
+100.88 6327.12 256.543
+101.84 6327.93 257.151
+102.35 6359.9  257.785
+102.83 6393.5  258.516
+103.51 6476.86 259.191
+104.13 6524.5  259.738
+104.71 6600.31 260.351
+105.39 6629.47 261.04
+106.09 6688.61 261.692
+106.75 6717.46 262.236
+107.24 6724.2  262.847
+107.75 6779.53 263.527
+108.29 6825.8  264.169
+108.91 6882    264.681
+109.24 6983.91 265.258
+109.74 7020    265.887
+110.23 7093.12 266.491
+111    7166.68 266.987
+111.43 7236.5  267.545
+111.76 7311.24 268.171
+112.08 7364.63 268.815];
+
+%Compute growth rates: from 1950:I to 1997:IV
+gy_obs=1000*data_q(:,2)./data_q(:,3); %real GDP per capita
+gy_obs=diff(log(gy_obs)); 
+
+gp_obs = diff(log(data_q(:,1))); %GDP deflator inflation

license.txt

@@ -245,6 +245,11 @@ Copyright: 2005-2010 Pascal Getreuer
            2017 Dynare Team
 License: BSD-2-clause
 
+Files: examples/fs2000_data.m
+Copyright: 2000-2022 Frank Schorfheide
+Copyright: 2023 Dynare Team
+License: CC-BY-SA-4.0
+
 Files: doc/*.rst doc/*.tex doc/*.svg doc/*.pdf doc/*.bib
 Copyright: 1996-2022 Dynare Team
 License: GFDL-NIV-1.3+

matlab/distributions/inverse_gamma_specification.m

@@ -85,7 +85,7 @@ s = [];
 nu = [];
 
 sigma = sqrt(sigma2);
-mu2 = mu*mu;
+mu2 = (mu-lb)*(mu-lb);
 
 if type == 2       % Inverse Gamma 2
     nu   = 2*(2+mu2/sigma2);
@@ -132,10 +132,10 @@ elseif type == 1   % Inverse Gamma 1
         end
         s = (sigma2+mu2)*(nu-2);
         if check_solution_flag
-            if abs(log(mu)-log(sqrt(s/2))-gammaln((nu-1)/2)+gammaln(nu/2))>1e-7
+            if abs(log(mu-lb)-log(sqrt(s/2))-gammaln((nu-1)/2)+gammaln(nu/2))>1e-7
                 error('inverse_gamma_specification:: Failed in solving for the hyperparameters!');
             end
-            if abs(sigma-sqrt(s/(nu-2)-mu*mu))>1e-7
+            if abs(sigma-sqrt(s/(nu-2)-(mu-lb)*(mu-lb)))>1e-7
                 error('inverse_gamma_specification:: Failed in solving for the hyperparameters!');
             end
         end

matlab/dynare_estimation_init.m

@@ -498,18 +498,18 @@ if options_.heteroskedastic_filter
             'for heteroskedastic_filter'])
     end
 
-    M_.heteroskedastic_shocks.Qvalue = NaN(M_.exo_nbr,options_.nobs);
-    M_.heteroskedastic_shocks.Qscale = NaN(M_.exo_nbr,options_.nobs);
+    M_.heteroskedastic_shocks.Qvalue = NaN(M_.exo_nbr,options_.nobs+1);
+    M_.heteroskedastic_shocks.Qscale = NaN(M_.exo_nbr,options_.nobs+1);
 
     for k=1:length(M_.heteroskedastic_shocks.Qvalue_orig)
         v = M_.heteroskedastic_shocks.Qvalue_orig(k);
-        temp_periods=v.periods(v.periods<options_.nobs+options_.first_obs);
+        temp_periods=v.periods(v.periods<options_.nobs+options_.first_obs+1);
         temp_periods=temp_periods(temp_periods>=options_.first_obs);
         M_.heteroskedastic_shocks.Qvalue(v.exo_id, temp_periods-(options_.first_obs-1)) = v.value^2;
     end
     for k=1:length(M_.heteroskedastic_shocks.Qscale_orig)
         v = M_.heteroskedastic_shocks.Qscale_orig(k);
-        temp_periods=v.periods(v.periods<options_.nobs+options_.first_obs);
+        temp_periods=v.periods(v.periods<options_.nobs+options_.first_obs+1);
         temp_periods=temp_periods(temp_periods>=options_.first_obs);
         M_.heteroskedastic_shocks.Qscale(v.exo_id, temp_periods-(options_.first_obs-1)) = v.scale^2;
     end

matlab/kalman/get_Qvec_heteroskedastic_filter.m

@@ -27,7 +27,7 @@ function Qvec=get_Qvec_heteroskedastic_filter(Q,smpl,M_)
 
 isqdiag = all(all(abs(Q-diag(diag(Q)))<1e-14)); % ie, the covariance matrix is diagonal...
 Qvec=repmat(Q,[1 1 smpl+1]);
-for k=1:smpl
+for k=1:smpl+1
     inx = ~isnan(M_.heteroskedastic_shocks.Qvalue(:,k));
     if any(inx)
         if isqdiag

matlab/model_diagnostics.m

@@ -213,7 +213,7 @@ for b=1:nb
             if n_rel  > 1
                 disp(['Relation ' int2str(i)])
             end
-            disp('Colinear variables:')
+            disp('Collinear variables:')
             for j=1:10
                 k = find(abs(ncol(:,i)) > 10^-j);
                 if max(abs(jacob(:,k)*ncol(k,i))) < 1e-6
@@ -236,7 +236,7 @@ for b=1:nb
             if n_rel  > 1
                 disp(['Relation ' int2str(i)])
             end
-            disp('Colinear equations')
+            disp('Collinear equations')
             for j=1:10
                 k = find(abs(neq(:,i)) > 10^-j);
                 if max(abs(jacob(k,:)'*neq(k,i))) < 1e-6

matlab/perfect-foresight-models/sim1.m

@@ -1,4 +1,5 @@
 function [endogenousvariables, success, err, iter] = sim1(endogenousvariables, exogenousvariables, steadystate, M_, options_)
+% [endogenousvariables, success, err, iter] = sim1(endogenousvariables, exogenousvariables, steadystate, M_, options_)
 % Performs deterministic simulations with lead or lag of one period, using
 % a basic Newton solver on sparse matrices.
 % Uses perfect_foresight_problem DLL to construct the stacked problem.
@@ -90,6 +91,7 @@ for iter = 1:options_.simul.maxit
                 end
             end            
         end
+        check_Jacobian_for_singularity(full(A),M_.endo_names,options_);
     end
     if options_.endogenous_terminal_period && iter > 1
         for it = 1:periods
@@ -281,3 +283,64 @@ if any(isinf(dyy))
     fprintf('%s, ', endo_names{:}),
     fprintf('\n'),
 end
+
+
+function check_Jacobian_for_singularity(jacob,endo_names,options_)
+n_vars_jacob=size(jacob,2);
+try
+    if (~isoctave && matlab_ver_less_than('9.12')) || isempty(options_.jacobian_tolerance)
+        rank_jacob = rank(jacob); %can sometimes fail
+    else
+        rank_jacob = rank(jacob,options_.jacobian_tolerance); %can sometimes fail
+    end
+catch
+    rank_jacob=size(jacob,1);
+end
+if rank_jacob < size(jacob,1)
+    disp(['sim1:  The Jacobian of the dynamic model is ' ...
+        'singular'])
+    disp(['sim1:  there is ' num2str(n_vars_jacob-rank_jacob) ...
+        ' collinear relationships between the variables and the equations'])
+    if (~isoctave && matlab_ver_less_than('9.12')) || isempty(options_.jacobian_tolerance)
+        ncol = null(jacob);
+    else
+        ncol = null(jacob,options_.jacobian_tolerance); %can sometimes fail
+    end
+    n_rel = size(ncol,2);
+    for i = 1:n_rel
+        if n_rel  > 1
+            disp(['Relation ' int2str(i)])
+        end
+        disp('Collinear variables:')
+        for j=1:10
+            k = find(abs(ncol(:,i)) > 10^-j);
+            if max(abs(jacob(:,k)*ncol(k,i))) < 1e-6
+                break
+            end
+        end
+        fprintf('%s\n',endo_names{mod(k,length(endo_names))})
+    end
+    if (~isoctave && matlab_ver_less_than('9.12')) || isempty(options_.jacobian_tolerance)
+        neq = null(jacob'); %can sometimes fail
+    else
+        neq = null(jacob',options_.jacobian_tolerance); %can sometimes fail
+    end
+    n_rel = size(neq,2);
+    for i = 1:n_rel
+        if n_rel  > 1
+            disp(['Relation ' int2str(i)])
+        end
+        disp('Collinear equations')
+        for j=1:10
+            k = find(abs(neq(:,i)) > 10^-j);
+            if max(abs(jacob(k,:)'*neq(k,i))) < 1e-6
+                break
+            end
+        end
+        equation=mod(k,length(endo_names));
+        period=ceil(k/length(endo_names));
+        for ii=1:length(equation)
+            fprintf('Equation %5u, period %5u\n',equation(ii),period(ii))
+        end
+    end
+end

matlab/posterior_sampler_core.m

@@ -141,49 +141,38 @@ for curr_block = fblck:nblck
         options_=set_dynare_seed_local_options(options_,options_.DynareRandomStreams.seed+curr_block);
     end
     mh_recover_flag=0;
-    if (options_.load_mh_file~=0) && (fline(curr_block)>1) && OpenOldFile(curr_block) %load previous draws and likelihood
-        load([BaseName '_mh' int2str(NewFile(curr_block)) '_blck' int2str(curr_block) '.mat'])
-        x2 = [x2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,npar)];
-        logpo2 = [logpo2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,1)];
+    if options_.mh_recover && exist([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat'],'file')==2 && OpenOldFile(curr_block)
+        % this should be done whatever value of load_mh_file
+        load([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat']);
+        draw_iter = size(neval_this_chain,2)+1;
+        draw_index_current_file = draw_iter+fline(curr_block)-1;
+        feval_this_chain = sum(sum(neval_this_chain));
+        feval_this_file = sum(sum(neval_this_chain));
+        if feval_this_chain>draw_index_current_file-fline(curr_block)
+            % non Metropolis type of sampler
+            accepted_draws_this_chain = draw_index_current_file-fline(curr_block);
+            accepted_draws_this_file = draw_index_current_file-fline(curr_block);
+        else
+            accepted_draws_this_chain = 0;
+            accepted_draws_this_file = 0;
+        end
+        mh_recover_flag=1;
+        set_dynare_random_generator_state(LastSeeds.(['file' int2str(NewFile(curr_block))]).Unifor, LastSeeds.(['file' int2str(NewFile(curr_block))]).Normal);
+        last_draw(curr_block,:)=x2(draw_index_current_file-1,:);
+        last_posterior(curr_block)=logpo2(draw_index_current_file-1);
         OpenOldFile(curr_block) = 0;
     else
-        if options_.mh_recover && exist([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat'],'file')==2
-            load([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat']);
-            draw_iter = size(neval_this_chain,2)+1;
-            draw_index_current_file = draw_iter;
-            feval_this_chain = sum(sum(neval_this_chain));
-            feval_this_file = sum(sum(neval_this_chain));
-            if feval_this_chain>draw_iter-1
-                % non Metropolis type of sampler
-                accepted_draws_this_chain = draw_iter-1;
-                accepted_draws_this_file = draw_iter-1;
-            else
-                accepted_draws_this_chain = 0;
-                accepted_draws_this_file = 0;
-            end
-            mh_recover_flag=1;
-            set_dynare_random_generator_state(LastSeeds.(['file' int2str(NewFile(curr_block))]).Unifor, LastSeeds.(['file' int2str(NewFile(curr_block))]).Normal);
-            last_draw(curr_block,:)=x2(draw_iter-1,:);
-            last_posterior(curr_block)=logpo2(draw_iter-1);
-
+        if (options_.load_mh_file~=0) && (fline(curr_block)>1) && OpenOldFile(curr_block) %load previous draws and likelihood
+            load([BaseName '_mh' int2str(NewFile(curr_block)) '_blck' int2str(curr_block) '.mat'])
+            x2 = [x2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,npar)];
+            logpo2 = [logpo2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,1)];
+            OpenOldFile(curr_block) = 0;
         else
 
             x2 = zeros(InitSizeArray(curr_block),npar);
             logpo2 = zeros(InitSizeArray(curr_block),1);
         end
     end
-    %Prepare waiting bars
-    if whoiam
-        refresh_rate = sampler_options.parallel_bar_refresh_rate;
-        bar_title = sampler_options.parallel_bar_title;
-        prc0=(curr_block-fblck)/(nblck-fblck+1)*(isoctave || options_.console_mode);
-        hh_fig = dyn_waitbar({prc0,whoiam,options_.parallel(ThisMatlab)},[bar_title ' (' int2str(curr_block) '/' int2str(options_.mh_nblck) ')...']);
-    else
-        refresh_rate = sampler_options.serial_bar_refresh_rate;
-        bar_title = sampler_options.serial_bar_title;
-        hh_fig = dyn_waitbar(0,[bar_title ' (' int2str(curr_block) '/' int2str(options_.mh_nblck) ')...']);
-        set(hh_fig,'Name',bar_title);
-    end
     if mh_recover_flag==0
         accepted_draws_this_chain = 0;
         accepted_draws_this_file = 0;
@@ -192,6 +181,18 @@ for curr_block = fblck:nblck
         draw_iter = 1;
         draw_index_current_file = fline(curr_block); %get location of first draw in current block
     end
+    %Prepare waiting bars
+    if whoiam
+        refresh_rate = sampler_options.parallel_bar_refresh_rate;
+        bar_title = sampler_options.parallel_bar_title;
+        prc0=(curr_block-fblck)/(nblck-fblck+1)*(isoctave || options_.console_mode)+(draw_iter-1)/nruns(curr_block);
+        hh_fig = dyn_waitbar({prc0,whoiam,options_.parallel(ThisMatlab)},[bar_title ' (' int2str(curr_block) '/' int2str(options_.mh_nblck) ')...']);
+    else
+        refresh_rate = sampler_options.serial_bar_refresh_rate;
+        bar_title = sampler_options.serial_bar_title;
+        hh_fig = dyn_waitbar(0,[bar_title ' (' int2str(curr_block) '/' int2str(options_.mh_nblck) ')...']);
+        set(hh_fig,'Name',bar_title);
+    end
     sampler_options.curr_block = curr_block;
     while draw_iter <= nruns(curr_block)
 

matlab/posterior_sampler_initialization.m

@@ -170,6 +170,12 @@ if ~options_.load_mh_file && ~options_.mh_recover
         ilogpo2 = zeros(NumberOfBlocks,1);
         ilogpo2(:,1) = record0.LastLogPost;
         ix2(:,IA) = record0.LastParameters(:,IB);
+        for j=1:NumberOfBlocks
+            if not(all(ix2(j,:)' >= mh_bounds.lb) && all(ix2(j,:)' <= mh_bounds.ub))
+                new_estimated_parameters = logical(new_estimated_parameters + (ix2(j,:)' < mh_bounds.lb));
+                new_estimated_parameters = logical(new_estimated_parameters + (ix2(j,:)' > mh_bounds.ub));
+            end
+        end
     else
         new_estimated_parameters = true(1,npar);
     end
@@ -458,7 +464,7 @@ elseif options_.mh_recover
     AllMhFiles = dir([BaseName '_mh*_blck*.mat']);
     TotalNumberOfMhFiles = length(AllMhFiles)-length(dir([BaseName '_mh_tmp*_blck*.mat']));
     % Quit if no crashed mcmc chain can be found as there are as many files as expected
-    if (TotalNumberOfMhFiles==ExpectedNumberOfMhFiles)
+    if (ExpectedNumberOfMhFilesPerBlock>LastFileNumberInThePreviousMh) && (TotalNumberOfMhFiles==ExpectedNumberOfMhFiles)
         if isnumeric(options_.parallel)
             fprintf('%s: It appears that you don''t need to use the mh_recover option!\n',dispString);
             fprintf('                  You have to edit the mod file and remove the mh_recover option\n');
@@ -469,15 +475,23 @@ elseif options_.mh_recover
     % 2. Something needs to be done; find out what
     % Count the number of saved mh files per block.
     NumberOfMhFilesPerBlock = zeros(NumberOfBlocks,1);
+    is_chain_complete = true(NumberOfBlocks,1);
     for b = 1:NumberOfBlocks
         BlckMhFiles = dir([BaseName '_mh*_blck' int2str(b) '.mat']);
         NumberOfMhFilesPerBlock(b) = length(BlckMhFiles)-length(dir([BaseName '_mh_tmp*_blck' int2str(b) '.mat']));
+        if (ExpectedNumberOfMhFilesPerBlock==LastFileNumberInThePreviousMh)
+            % here we need to check for the number of lines
+            tmpdata = load([BaseName '_mh' int2str(LastFileNumberInThePreviousMh) '_blck' int2str(b) '.mat'],'logpo2');
+            if length(tmpdata.logpo2) == LastLineNumberInThePreviousMh
+                is_chain_complete(b) = false;
+            end
+        end
     end
     % Find FirstBlock (First block), an integer targeting the crashed mcmc chain.
     FirstBlock = 1; %initialize
     FBlock = zeros(NumberOfBlocks,1);
     while FirstBlock <= NumberOfBlocks
-        if  NumberOfMhFilesPerBlock(FirstBlock) < ExpectedNumberOfMhFilesPerBlock
+        if  (NumberOfMhFilesPerBlock(FirstBlock) < ExpectedNumberOfMhFilesPerBlock) || not(is_chain_complete(FirstBlock))
             fprintf('%s: Chain %u is not complete!\n', dispString,FirstBlock);
             FBlock(FirstBlock)=1;
         else