1function [raftery_lewis] = raftery_lewis(runs,q,r,s) 2% function raftery_lewis = raftery_lewis(runs,q,r,s) 3% Computes the convergence diagnostics of Raftery and Lewis (1992), i.e. the 4% number of draws needed in MCMC to estimate the posterior cdf of the q-quantile 5% within an accuracy r with probability s 6% 7% Inputs: 8% - draws [n_draws by n_var] double matrix of draws from the sampler 9% - q [scalar] quantile of the quantity of interest 10% - r [scalar] level of desired precision 11% - s [scalar] probability associated with r 12% 13% Output: 14% raftery_lewis [structure] containing the fields: 15% - M_burn [n_draws by 1] number of draws required for burn-in 16% - N_prec [n_draws by 1] number of draws required to achieve desired precision r 17% - k_thin [n_draws by 1] thinning required to get 1st order MC 18% - k_ind [n_draws by 1] thinning required to get independence 19% - I_stat [n_draws by 1] I-statistic of Raftery/Lewis (1992b) 20% measures increase in required 21% iterations due to dependence in chain 22% - N_min [scalar] # draws if the chain is white noise 23% - N_total [n_draws by 1] nburn + nprec 24% 25 26% --------------------------------------------------------------------- 27% NOTES: Example values of q, r, s: 28% 0.025, 0.005, 0.95 (for a long-tailed distribution) 29% 0.025, 0.0125, 0.95 (for a short-tailed distribution); 30% 31% - The result is quite sensitive to r, being proportional to the 32% inverse of r^2. 33% - For epsilon (closeness of probabilities to equilibrium values), 34% Raftery/Lewis use 0.001 and argue that the results 35% are quite robust to changes in this value 36% 37% --------------------------------------------------------------------- 38% REFERENCES: 39% Raftery, Adrien E./Lewis, Steven (1992a): "How many iterations in the Gibbs sampler?" 40% in: Bernardo/Berger/Dawid/Smith (eds.): Bayesian Statistics, Vol. 4, Clarendon Press: Oxford, 41% pp. 763-773. 42% Raftery, Adrien E./Lewis, Steven (1992b): "Comment: One long run with diagnostics: 43% Implementation strategies for Markov chain Monte Carlo." Statistical Science, 44% 7(4), pp. 493-497. 45% 46% ---------------------------------------------------- 47 48% Copyright (C) 2016 Benjamin Born and Johannes Pfeifer 49% Copyright (C) 2016-2017 Dynare Team 50% 51% This file is part of Dynare. 52% 53% Dynare is free software: you can redistribute it and/or modify 54% it under the terms of the GNU General Public License as published by 55% the Free Software Foundation, either version 3 of the License, or 56% (at your option) any later version. 57% 58% Dynare is distributed in the hope that it will be useful, 59% but WITHOUT ANY WARRANTY; without even the implied warranty of 60% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 61% GNU General Public License for more details. 62% 63% You should have received a copy of the GNU General Public License 64% along with Dynare. If not, see <http://www.gnu.org/licenses/>. 65 66 67 68[n_runs, n_vars] = size(runs); 69 70raftery_lewis.M_burn=NaN(n_vars,1); 71raftery_lewis.N_prec=NaN(n_vars,1); 72raftery_lewis.k_thin=NaN(n_vars,1); 73raftery_lewis.k_ind=NaN(n_vars,1); 74raftery_lewis.I_stat=NaN(n_vars,1); 75raftery_lewis.N_total=NaN(n_vars,1); 76 77 78thinned_chain = zeros(n_runs,1); 79%quantities that can be precomputed as they are independent of variable 80Phi = norminv((s+1)/2); %note the missing ^{-1} at the Phi in equation top page 5, see RL (1995) 81raftery_lewis.N_min = fix(Phi^2*(1-q)*q/r^2+1); 82 83for nv = 1:n_vars % big loop over variables 84 if q > 0 && q < 1 85 work = (runs(:,nv) <= quantile(runs(:,nv),q)); 86 else 87 error('Quantile must be between 0 and 1'); 88 end 89 90 k_thin_current_var = 1; 91 bic = 1; 92 epss = 0.001; 93 % Find thinning factor for which first-order Markov Chain is preferred to second-order one 94 while(bic > 0) 95 thinned_chain=work(1:k_thin_current_var:n_runs,1); 96 [g2, bic] = first_vs_second_order_MC_test(thinned_chain); 97 k_thin_current_var = k_thin_current_var+1; 98 end 99 100 k_thin_current_var = k_thin_current_var-1; %undo last step 101 102 %compute transition probabilities 103 transition_matrix = zeros(2,2); 104 for i1 = 2:size(thinned_chain,1) 105 transition_matrix(thinned_chain(i1-1)+1,thinned_chain(i1)+1) = transition_matrix(thinned_chain(i1-1)+1,thinned_chain(i1)+1)+1; 106 end 107 alpha = transition_matrix(1,2)/(transition_matrix(1,1)+transition_matrix(1,2)); %prob of going from 1 to 2 108 beta = transition_matrix(2,1)/(transition_matrix(2,1)+transition_matrix(2,2)); %prob of going from 2 to 1 109 110 kmind=k_thin_current_var; 111 [g2, bic]=independence_chain_test(thinned_chain); 112 113 while(bic > 0) 114 thinned_chain=work(1:kmind:n_runs,1); 115 [g2, bic] = independence_chain_test(thinned_chain); 116 kmind = kmind+1; 117 end 118 119 m_star = log((alpha + beta)*epss/max(alpha,beta))/log(abs(1 - alpha - beta)); %equation bottom page 4 120 raftery_lewis.M_burn(nv) = fix((m_star+1)*k_thin_current_var); 121 n_star = (2 - (alpha + beta))*alpha*beta*(Phi^2)/((alpha + beta)^3 * r^2); %equation top page 5 122 raftery_lewis.N_prec(nv) = fix(n_star+1)*k_thin_current_var; 123 raftery_lewis.I_stat(nv) = (raftery_lewis.M_burn(nv) + raftery_lewis.N_prec(nv))/raftery_lewis.N_min; 124 raftery_lewis.k_ind(nv) = max(fix(raftery_lewis.I_stat(nv)+1),kmind); 125 raftery_lewis.k_thin(nv) = k_thin_current_var; 126 raftery_lewis.N_total(nv)= raftery_lewis.M_burn(nv)+raftery_lewis.N_prec(nv); 127end 128 129end 130 131function [g2, bic] = first_vs_second_order_MC_test(d) 132%conducts a test of first vs. second order Markov Chain via BIC criterion 133n_obs=size(d,1); 134g2 = 0; 135tran=zeros(2,2,2); 136for t_iter=3:n_obs % count state transitions 137 tran(d(t_iter-2,1)+1,d(t_iter-1,1)+1,d(t_iter,1)+1)=tran(d(t_iter-2,1)+1,d(t_iter-1,1)+1,d(t_iter,1)+1)+1; 138end 139% Compute the log likelihood ratio statistic for second-order MC vs first-order MC. G2 statistic of Bishop, Fienberg and Holland (1975) 140for ind_1 = 1:2 141 for ind_2 = 1:2 142 for ind_3 = 1:2 143 if tran(ind_1,ind_2,ind_3) ~= 0 144 fitted = (tran(ind_1,ind_2,1) + tran(ind_1,ind_2,2))*(tran(1,ind_2,ind_3) + tran(2,ind_2,ind_3))/... 145 (tran(1,ind_2,1) + tran(1,ind_2,2) + tran(2,ind_2,1) + tran(2,ind_2,2)); 146 focus = tran(ind_1,ind_2,ind_3); 147 g2 = g2 + log(focus/fitted)*focus; 148 end 149 end % end of for i3 150 end % end of for i2 151end % end of for i1 152g2 = g2*2; 153bic = g2 - log(n_obs-2)*2; 154 155end 156 157 158function [g2, bic] = independence_chain_test(d) 159%conducts a test of independence Chain via BIC criterion 160n_obs=size(d,1); 161trans = zeros(2,2); 162for ind_1 = 2:n_obs 163 trans(d(ind_1-1)+1,d(ind_1)+1)=trans(d(ind_1-1)+1,d(ind_1)+1)+1; 164end 165dcm1 = n_obs - 1; 166g2 = 0; 167% Compute the log likelihood ratio statistic for second-order MC vs first-order MC. G2 statistic of Bishop, Fienberg and Holland (1975) 168for ind_1 = 1:2 169 for ind_2 = 1:2 170 if trans(ind_1,ind_2) ~= 0 171 fitted = ((trans(ind_1,1) + trans(ind_1,2))*(trans(1,ind_2) + trans(2,ind_2)))/dcm1; 172 focus = trans(ind_1,ind_2); 173 g2 = g2 + log(focus/fitted)*focus; 174 end 175 end 176end 177g2 = g2*2; 178bic = g2 - log(dcm1); 179end 180