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