Compare commits
11 Commits
a7f5fd571d
...
7cdb09991f
Author | SHA1 | Date |
---|---|---|
Johannes Pfeifer | 7cdb09991f | |
Sébastien Villemot | 441ef7e102 | |
Marco Ratto | f102a992aa | |
Marco Ratto | 53b57da8ba | |
Marco Ratto | aad5c36081 | |
Marco Ratto | de152a3de3 | |
Marco Ratto | 8f73564634 | |
Sébastien Villemot | 858b534c22 | |
Johannes Pfeifer | e17bf15042 | |
Sébastien Villemot | ea28fcb4b4 | |
Johannes Pfeifer | 42fc1ec40a |
|
@ -3,20 +3,20 @@
|
||||||
* in the paper) described in Frank Schorfheide (2000): "Loss function-based
|
* in the paper) described in Frank Schorfheide (2000): "Loss function-based
|
||||||
* evaluation of DSGE models", Journal of Applied Econometrics, 15(6), 645-670.
|
* 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.
|
* The data are taken from the replication package at
|
||||||
* They are therefore different from the original dataset used by Schorfheide.
|
* http://dx.doi.org/10.15456/jae.2022314.0708799949
|
||||||
*
|
*
|
||||||
* The prior distribution follows the one originally specified in Schorfheide's
|
* 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
|
* implies an asymptote and corresponding prior mode at 0. It is generally
|
||||||
* recommended to avoid this extreme type of prior. Some optimizers, for instance
|
* recommended to avoid this extreme type of prior.
|
||||||
* 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
|
* Because the data are already logged and we use the loglinear option to conduct
|
||||||
* (changing .223 to .100) to remove the asymptote.
|
* a full log-linearization, we need to use the logdata option.
|
||||||
*
|
*
|
||||||
* The equations are taken from J. Nason and T. Cogley (1994): "Testing the
|
* The equations are taken from J. Nason and T. Cogley (1994): "Testing the
|
||||||
* implications of long-run neutrality for monetary business cycle models",
|
* 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.
|
* 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
|
* 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.
|
* This file is part of Dynare.
|
||||||
*
|
*
|
||||||
|
@ -43,33 +43,71 @@
|
||||||
* along with Dynare. If not, see <https://www.gnu.org/licenses/>.
|
* 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;
|
var m ${m}$ (long_name='money growth')
|
||||||
varexo e_a e_m;
|
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;
|
alp = 0.33;
|
||||||
bet = 0.99;
|
bet = 0.99;
|
||||||
gam = 0.003;
|
gam = 0.003;
|
||||||
mst = 1.011;
|
mst = 1.011;
|
||||||
rho = 0.7;
|
rho = 0.7;
|
||||||
psi = 0.787;
|
phi = 0.787;
|
||||||
del = 0.02;
|
del = 0.02;
|
||||||
|
|
||||||
model;
|
model;
|
||||||
|
[name='NC before eq. (1), TFP growth equation']
|
||||||
dA = exp(gam+e_a);
|
dA = exp(gam+e_a);
|
||||||
|
[name='NC eq. (2), money growth rate']
|
||||||
log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
|
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;
|
-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;
|
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;
|
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;
|
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);
|
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;
|
P*c = m;
|
||||||
|
[name='NC eq. (20), credit market equilibrium condition']
|
||||||
m-1+d = l;
|
m-1+d = l;
|
||||||
|
[name='Definition TFP shock']
|
||||||
e = exp(e_a);
|
e = exp(e_a);
|
||||||
|
[name='Implied by NC eq. (18), production function']
|
||||||
y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
|
y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
|
||||||
|
[name='Observation equation GDP growth']
|
||||||
gy_obs = dA*y/y(-1);
|
gy_obs = dA*y/y(-1);
|
||||||
|
[name='Observation equation price level']
|
||||||
gp_obs = (P/P(-1))*m(-1)/dA;
|
gp_obs = (P/P(-1))*m(-1)/dA;
|
||||||
end;
|
end;
|
||||||
|
|
||||||
|
@ -84,12 +122,12 @@ steady_state_model;
|
||||||
m = mst;
|
m = mst;
|
||||||
khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1));
|
khst = ( (1-gst*bet*(1-del)) / (alp*gst^alp*bet) )^(1/(alp-1));
|
||||||
xist = ( ((khst*gst)^alp - (1-gst*(1-del))*khst)/mst )^(-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);
|
n = xist/(nust+xist);
|
||||||
P = xist + nust;
|
P = xist + nust;
|
||||||
k = khst*n;
|
k = khst*n;
|
||||||
|
|
||||||
l = psi*mst*n/( (1-psi)*(1-n) );
|
l = phi*mst*n/( (1-phi)*(1-n) );
|
||||||
c = mst/P;
|
c = mst/P;
|
||||||
d = l - mst + 1;
|
d = l - mst + 1;
|
||||||
y = k^alp*n^(1-alp)*gst^alp;
|
y = k^alp*n^(1-alp)*gst^alp;
|
||||||
|
@ -104,17 +142,18 @@ steady_state_model;
|
||||||
gy_obs = dA;
|
gy_obs = dA;
|
||||||
end;
|
end;
|
||||||
|
|
||||||
steady;
|
|
||||||
|
|
||||||
|
steady;
|
||||||
check;
|
check;
|
||||||
|
|
||||||
|
% Table 1 of Schorfheide (2000)
|
||||||
estimated_params;
|
estimated_params;
|
||||||
alp, beta_pdf, 0.356, 0.02;
|
alp, beta_pdf, 0.356, 0.02;
|
||||||
bet, beta_pdf, 0.993, 0.002;
|
bet, beta_pdf, 0.993, 0.002;
|
||||||
gam, normal_pdf, 0.0085, 0.003;
|
gam, normal_pdf, 0.0085, 0.003;
|
||||||
mst, normal_pdf, 1.0002, 0.007;
|
mst, normal_pdf, 1.0002, 0.007;
|
||||||
rho, beta_pdf, 0.129, 0.100;
|
rho, beta_pdf, 0.129, 0.223;
|
||||||
psi, beta_pdf, 0.65, 0.05;
|
phi, beta_pdf, 0.65, 0.05;
|
||||||
del, beta_pdf, 0.01, 0.005;
|
del, beta_pdf, 0.01, 0.005;
|
||||||
stderr e_a, inv_gamma_pdf, 0.035449, inf;
|
stderr e_a, inv_gamma_pdf, 0.035449, inf;
|
||||||
stderr e_m, inv_gamma_pdf, 0.008862, inf;
|
stderr e_m, inv_gamma_pdf, 0.008862, inf;
|
||||||
|
@ -122,14 +161,8 @@ end;
|
||||||
|
|
||||||
varobs gp_obs gy_obs;
|
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);
|
||||||
|
|
||||||
|
%uncomment the following lines to generate LaTeX-code of the model equations
|
||||||
/*
|
%write_latex_original_model(write_equation_tags);
|
||||||
* The following lines were used to generate the data file. If you want to
|
%collect_latex_files;
|
||||||
* 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'});
|
|
|
@ -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
|
|
@ -245,6 +245,11 @@ Copyright: 2005-2010 Pascal Getreuer
|
||||||
2017 Dynare Team
|
2017 Dynare Team
|
||||||
License: BSD-2-clause
|
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
|
Files: doc/*.rst doc/*.tex doc/*.svg doc/*.pdf doc/*.bib
|
||||||
Copyright: 1996-2022 Dynare Team
|
Copyright: 1996-2022 Dynare Team
|
||||||
License: GFDL-NIV-1.3+
|
License: GFDL-NIV-1.3+
|
||||||
|
|
|
@ -85,7 +85,7 @@ s = [];
|
||||||
nu = [];
|
nu = [];
|
||||||
|
|
||||||
sigma = sqrt(sigma2);
|
sigma = sqrt(sigma2);
|
||||||
mu2 = mu*mu;
|
mu2 = (mu-lb)*(mu-lb);
|
||||||
|
|
||||||
if type == 2 % Inverse Gamma 2
|
if type == 2 % Inverse Gamma 2
|
||||||
nu = 2*(2+mu2/sigma2);
|
nu = 2*(2+mu2/sigma2);
|
||||||
|
@ -132,10 +132,10 @@ elseif type == 1 % Inverse Gamma 1
|
||||||
end
|
end
|
||||||
s = (sigma2+mu2)*(nu-2);
|
s = (sigma2+mu2)*(nu-2);
|
||||||
if check_solution_flag
|
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!');
|
error('inverse_gamma_specification:: Failed in solving for the hyperparameters!');
|
||||||
end
|
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!');
|
error('inverse_gamma_specification:: Failed in solving for the hyperparameters!');
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
|
|
@ -498,18 +498,18 @@ if options_.heteroskedastic_filter
|
||||||
'for heteroskedastic_filter'])
|
'for heteroskedastic_filter'])
|
||||||
end
|
end
|
||||||
|
|
||||||
M_.heteroskedastic_shocks.Qvalue = 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);
|
M_.heteroskedastic_shocks.Qscale = NaN(M_.exo_nbr,options_.nobs+1);
|
||||||
|
|
||||||
for k=1:length(M_.heteroskedastic_shocks.Qvalue_orig)
|
for k=1:length(M_.heteroskedastic_shocks.Qvalue_orig)
|
||||||
v = M_.heteroskedastic_shocks.Qvalue_orig(k);
|
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);
|
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;
|
M_.heteroskedastic_shocks.Qvalue(v.exo_id, temp_periods-(options_.first_obs-1)) = v.value^2;
|
||||||
end
|
end
|
||||||
for k=1:length(M_.heteroskedastic_shocks.Qscale_orig)
|
for k=1:length(M_.heteroskedastic_shocks.Qscale_orig)
|
||||||
v = M_.heteroskedastic_shocks.Qscale_orig(k);
|
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);
|
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;
|
M_.heteroskedastic_shocks.Qscale(v.exo_id, temp_periods-(options_.first_obs-1)) = v.scale^2;
|
||||||
end
|
end
|
||||||
|
|
|
@ -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...
|
isqdiag = all(all(abs(Q-diag(diag(Q)))<1e-14)); % ie, the covariance matrix is diagonal...
|
||||||
Qvec=repmat(Q,[1 1 smpl+1]);
|
Qvec=repmat(Q,[1 1 smpl+1]);
|
||||||
for k=1:smpl
|
for k=1:smpl+1
|
||||||
inx = ~isnan(M_.heteroskedastic_shocks.Qvalue(:,k));
|
inx = ~isnan(M_.heteroskedastic_shocks.Qvalue(:,k));
|
||||||
if any(inx)
|
if any(inx)
|
||||||
if isqdiag
|
if isqdiag
|
||||||
|
|
|
@ -213,7 +213,7 @@ for b=1:nb
|
||||||
if n_rel > 1
|
if n_rel > 1
|
||||||
disp(['Relation ' int2str(i)])
|
disp(['Relation ' int2str(i)])
|
||||||
end
|
end
|
||||||
disp('Colinear variables:')
|
disp('Collinear variables:')
|
||||||
for j=1:10
|
for j=1:10
|
||||||
k = find(abs(ncol(:,i)) > 10^-j);
|
k = find(abs(ncol(:,i)) > 10^-j);
|
||||||
if max(abs(jacob(:,k)*ncol(k,i))) < 1e-6
|
if max(abs(jacob(:,k)*ncol(k,i))) < 1e-6
|
||||||
|
@ -236,7 +236,7 @@ for b=1:nb
|
||||||
if n_rel > 1
|
if n_rel > 1
|
||||||
disp(['Relation ' int2str(i)])
|
disp(['Relation ' int2str(i)])
|
||||||
end
|
end
|
||||||
disp('Colinear equations')
|
disp('Collinear equations')
|
||||||
for j=1:10
|
for j=1:10
|
||||||
k = find(abs(neq(:,i)) > 10^-j);
|
k = find(abs(neq(:,i)) > 10^-j);
|
||||||
if max(abs(jacob(k,:)'*neq(k,i))) < 1e-6
|
if max(abs(jacob(k,:)'*neq(k,i))) < 1e-6
|
||||||
|
|
|
@ -1,4 +1,5 @@
|
||||||
function [endogenousvariables, success, err, iter] = sim1(endogenousvariables, exogenousvariables, steadystate, M_, options_)
|
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
|
% Performs deterministic simulations with lead or lag of one period, using
|
||||||
% a basic Newton solver on sparse matrices.
|
% a basic Newton solver on sparse matrices.
|
||||||
% Uses perfect_foresight_problem DLL to construct the stacked problem.
|
% Uses perfect_foresight_problem DLL to construct the stacked problem.
|
||||||
|
@ -90,6 +91,7 @@ for iter = 1:options_.simul.maxit
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
end
|
end
|
||||||
|
check_Jacobian_for_singularity(full(A),M_.endo_names,options_);
|
||||||
end
|
end
|
||||||
if options_.endogenous_terminal_period && iter > 1
|
if options_.endogenous_terminal_period && iter > 1
|
||||||
for it = 1:periods
|
for it = 1:periods
|
||||||
|
@ -281,3 +283,64 @@ if any(isinf(dyy))
|
||||||
fprintf('%s, ', endo_names{:}),
|
fprintf('%s, ', endo_names{:}),
|
||||||
fprintf('\n'),
|
fprintf('\n'),
|
||||||
end
|
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
|
|
@ -141,49 +141,38 @@ for curr_block = fblck:nblck
|
||||||
options_=set_dynare_seed_local_options(options_,options_.DynareRandomStreams.seed+curr_block);
|
options_=set_dynare_seed_local_options(options_,options_.DynareRandomStreams.seed+curr_block);
|
||||||
end
|
end
|
||||||
mh_recover_flag=0;
|
mh_recover_flag=0;
|
||||||
if (options_.load_mh_file~=0) && (fline(curr_block)>1) && OpenOldFile(curr_block) %load previous draws and likelihood
|
if options_.mh_recover && exist([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat'],'file')==2 && OpenOldFile(curr_block)
|
||||||
load([BaseName '_mh' int2str(NewFile(curr_block)) '_blck' int2str(curr_block) '.mat'])
|
% this should be done whatever value of load_mh_file
|
||||||
x2 = [x2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,npar)];
|
load([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat']);
|
||||||
logpo2 = [logpo2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,1)];
|
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;
|
OpenOldFile(curr_block) = 0;
|
||||||
else
|
else
|
||||||
if options_.mh_recover && exist([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat'],'file')==2
|
if (options_.load_mh_file~=0) && (fline(curr_block)>1) && OpenOldFile(curr_block) %load previous draws and likelihood
|
||||||
load([BaseName '_mh_tmp_blck' int2str(curr_block) '.mat']);
|
load([BaseName '_mh' int2str(NewFile(curr_block)) '_blck' int2str(curr_block) '.mat'])
|
||||||
draw_iter = size(neval_this_chain,2)+1;
|
x2 = [x2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,npar)];
|
||||||
draw_index_current_file = draw_iter;
|
logpo2 = [logpo2;zeros(InitSizeArray(curr_block)-fline(curr_block)+1,1)];
|
||||||
feval_this_chain = sum(sum(neval_this_chain));
|
OpenOldFile(curr_block) = 0;
|
||||||
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);
|
|
||||||
|
|
||||||
else
|
else
|
||||||
|
|
||||||
x2 = zeros(InitSizeArray(curr_block),npar);
|
x2 = zeros(InitSizeArray(curr_block),npar);
|
||||||
logpo2 = zeros(InitSizeArray(curr_block),1);
|
logpo2 = zeros(InitSizeArray(curr_block),1);
|
||||||
end
|
end
|
||||||
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
|
if mh_recover_flag==0
|
||||||
accepted_draws_this_chain = 0;
|
accepted_draws_this_chain = 0;
|
||||||
accepted_draws_this_file = 0;
|
accepted_draws_this_file = 0;
|
||||||
|
@ -192,6 +181,18 @@ for curr_block = fblck:nblck
|
||||||
draw_iter = 1;
|
draw_iter = 1;
|
||||||
draw_index_current_file = fline(curr_block); %get location of first draw in current block
|
draw_index_current_file = fline(curr_block); %get location of first draw in current block
|
||||||
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)+(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;
|
sampler_options.curr_block = curr_block;
|
||||||
while draw_iter <= nruns(curr_block)
|
while draw_iter <= nruns(curr_block)
|
||||||
|
|
||||||
|
|
|
@ -170,6 +170,12 @@ if ~options_.load_mh_file && ~options_.mh_recover
|
||||||
ilogpo2 = zeros(NumberOfBlocks,1);
|
ilogpo2 = zeros(NumberOfBlocks,1);
|
||||||
ilogpo2(:,1) = record0.LastLogPost;
|
ilogpo2(:,1) = record0.LastLogPost;
|
||||||
ix2(:,IA) = record0.LastParameters(:,IB);
|
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
|
else
|
||||||
new_estimated_parameters = true(1,npar);
|
new_estimated_parameters = true(1,npar);
|
||||||
end
|
end
|
||||||
|
@ -458,7 +464,7 @@ elseif options_.mh_recover
|
||||||
AllMhFiles = dir([BaseName '_mh*_blck*.mat']);
|
AllMhFiles = dir([BaseName '_mh*_blck*.mat']);
|
||||||
TotalNumberOfMhFiles = length(AllMhFiles)-length(dir([BaseName '_mh_tmp*_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
|
% 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)
|
if isnumeric(options_.parallel)
|
||||||
fprintf('%s: It appears that you don''t need to use the mh_recover option!\n',dispString);
|
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');
|
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
|
% 2. Something needs to be done; find out what
|
||||||
% Count the number of saved mh files per block.
|
% Count the number of saved mh files per block.
|
||||||
NumberOfMhFilesPerBlock = zeros(NumberOfBlocks,1);
|
NumberOfMhFilesPerBlock = zeros(NumberOfBlocks,1);
|
||||||
|
is_chain_complete = true(NumberOfBlocks,1);
|
||||||
for b = 1:NumberOfBlocks
|
for b = 1:NumberOfBlocks
|
||||||
BlckMhFiles = dir([BaseName '_mh*_blck' int2str(b) '.mat']);
|
BlckMhFiles = dir([BaseName '_mh*_blck' int2str(b) '.mat']);
|
||||||
NumberOfMhFilesPerBlock(b) = length(BlckMhFiles)-length(dir([BaseName '_mh_tmp*_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
|
end
|
||||||
% Find FirstBlock (First block), an integer targeting the crashed mcmc chain.
|
% Find FirstBlock (First block), an integer targeting the crashed mcmc chain.
|
||||||
FirstBlock = 1; %initialize
|
FirstBlock = 1; %initialize
|
||||||
FBlock = zeros(NumberOfBlocks,1);
|
FBlock = zeros(NumberOfBlocks,1);
|
||||||
while FirstBlock <= NumberOfBlocks
|
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);
|
fprintf('%s: Chain %u is not complete!\n', dispString,FirstBlock);
|
||||||
FBlock(FirstBlock)=1;
|
FBlock(FirstBlock)=1;
|
||||||
else
|
else
|
||||||
|
|
Loading…
Reference in New Issue