*! runmlwin_mcmcdiag.ado, Chris Charlton and George Leckie, 08Sep2023 program runmlwin_mcmcdiag, rclass byable(recall) if _caller() >= 12 version 12.0 if _caller() <= 9 version 9.0 #delimit ; syntax varlist [if] [in] [, RLQUANT1(real 0.025) RLQUANT2(real 0.975) RLR(real 0.005) RLS(real 0.95) BDK(int 2) BDAlpha(real 0.05) KERNALPoints(int 1000) MCSEPoints(int 1000) ACFPoints(int 100) PACFPoints(int 10) DELTA(real 0.001) THINning(int 1) KSMOOth(real 1) MCSELl(real 0.5) MCSEUl(real 20) POSIT(int 0) PMean(real -9999) PVar(real -9999) Level(cilevel) Detail Eform * ]; #delimit cr tempname N quietly count `if' `in' scalar `N' = r(N) if "`detail'"=="" { // 1 2 3 4 5 6 7 //1234567890123456789012345678901234567890123456789012345678901234567890123456789 //-------------+------------------------------------------------------------- // Parameter | N Mean Std. Dev. ESS [95% Cred. Interval] //-------------+------------------------------------------------------------- // parameter | ####### ######### ######### ####### ######### ######### display _n(1) local c = length("`level'") display /// as txt _col(4) "Parameter" /// as txt _col(14) "{c |}" /// as txt _col(22) "N" /// as txt _col(30) "Mean" /// as txt _col(37) "Std. Dev." /// as txt _col(51) "ESS" /// as txt _col(`=58 -`c'') `"[`=strsubdp("`level'")'% Cred. Interval]"' display as txt "{hline 13}{c +}{hline 61}" } foreach parameter of local varlist { tempname mean scalar `mean' = . tempname sd scalar `sd' = . tempname mode scalar `mode' = . tempname quantiles local numquants = 11 matrix `quantiles' = J(`numquants', 2, .) mat `quantiles'[1,1] = .005 // 0.5% mat `quantiles'[2,1] = .025 // 2.5% mat `quantiles'[3,1] = .05 // 5% mat `quantiles'[4,1] = .25 // 2.5% mat `quantiles'[5,1] = .5 // 50% mat `quantiles'[6,1] = .75 // 95% mat `quantiles'[7,1] = .95 // 95% mat `quantiles'[8,1] = .975 // 97.5% mat `quantiles'[9,1] = .995 // 99.5% mat `quantiles'[10,1] = .01*((100 - `level')/2) // lb mat `quantiles'[11,1] = .01*(100 - (100 - `level')/2) // ub tempname ESS scalar `ESS' = . tempname RL1 // Raftery-Lewis Nhat1 scalar `RL1' = . tempname RL2 // Raftery-Lewis Nhat2 scalar `RL2' = . tempname BD // Brooks-Draper Nhat scalar `BD' = . tempname MEANMCSE // Brooks-Draper Nhat scalar `MEANMCSE' = . tempname KD if `kernalpoints' > 0 { mata: st_matrix("`KD'", J(`kernalpoints',2,.)) } tempname PRIOR if `kernalpoints' > 0 { mata: st_matrix("`PRIOR'", J(`kernalpoints',2,.)) } tempname MCSE if `mcsepoints' > 0 { mata: st_matrix("`MCSE'", J(`mcsepoints',2,.)) } tempname ACF if `acfpoints' > 0 { mata: st_matrix("`ACF'", J(`acfpoints',2,.)) } tempname PACF if `pacfpoints' > 0 { mata: st_matrix("`PACF'", J(`pacfpoints',2,.)) } // Do calculations plugin call runmlwin_mcmcdiagnostics `parameter' `if' `in', "`mean'" "`sd'" "`mode'" "`quantiles'" "`ESS'" "`RL1'" "`RL2'" "`BD'" "`KD'" "`MCSE'" "`ACF'" "`PACF'" "`MEANMCSE'" "`PRIOR'" local lb = `quantiles'[10,2] local ub = `quantiles'[11,2] matrix `quantiles' = `quantiles'[1..9,.] if "`detail'"~="" { // 1 2 3 4 5 6 7 //1234567890123456789012345678901234567890123456789012345678901234567890123456789 //-------------+------------------------------------------------------------- // Parameter | N Mean Std. Dev. ESS [95% Cred. Interval] //-------------+------------------------------------------------------------- // parameter | ####### ######### ######### ###### ######### ######### display _n(1) as res "`parameter'" display "-------------------------------------------------------------" display as text "Raftery-Lewis (quantile): Nhat = (" as result `RL1' as text "," as result `RL2' as text ")" display as text "when q = (`rlquant1', `rlquant2'), r = `rlr' and s = `rls'" display as text "Brooks-Draper (mean): Nhat = " as result `BD' display as text "when k = `bdk' and alpha = `bdalpha'" display as text "Column: `parameter', posterior mean = " as result `mean' as text "(" as result `MEANMCSE' as text ") SD = " as result `sd' as text " mode = " as result `mode' display as text "Quantiles:" forvalues i = 1/`numquants' { display as text `quantiles'[`i', 1] as result _col(10) `quantiles'[`i', 2] } display as text "`=_N * `thinning'' actual iterations. Effective Sample Size (ESS) = " as result `ESS' } else { local parname `=abbrev("`parameter'", 12)' local p = 13 - length("`parname'") display /// as txt _col(`p') "`parname'" /// as txt _col(14) "{c |}" /// as res _col(16) %7.0g `N' /// as res _col(25) %9.0g `mean' /// as res _col(36) %9.0g `sd' /// as res _col(48) %6.0g `ESS' /// as res _col(56) %9.0g `lb' /// as res _col(67) %9.0g `ub' local r = `r' + 1 } local compval 0 if "`eform'" ~= "" { local compval 1 } if `mean' > `compval' { quietly count if `parameter' < `compval' } else { quietly count if `parameter' > `compval' } local pvalmean = r(N) / `N' if `kernalpoints' > `compval' { if `mode' > 0 { quietly count if `parameter' < `compval' } else { quietly count if `parameter' > `compval' } local pvalmode = r(N) / `N' } if `quantiles'[5, 2] > `compval' { // NOTE: This assumes this is always 50% quietly count if `parameter' < `compval' } else { quietly count if `parameter' > `compval' } local pvalmedian = r(N) / `N' return clear return scalar N = `N' return scalar mean = `mean' return scalar sd = `sd' return scalar lb = `lb' return scalar ub = `ub' return matrix quantiles = `quantiles' return scalar RL1 = `RL1' return scalar RL2 = `RL2' return scalar pvalmean = `pvalmean' return scalar pvalmedian = `pvalmedian' if `kernalpoints' > 0 { return matrix KD = `KD' if ("`pmean'" != "-9999" & "`pvar'" != "-9999") return matrix prior = `PRIOR' return scalar mode = `mode' return scalar pvalmode = `pvalmode' } if `acfpoints' > 0 { return matrix ACF = `ACF' return scalar ESS = `ESS' return scalar BD = `BD' if `mcsepoints' > 0 { return matrix MCSE = `MCSE' return scalar meanmcse = `MEANMCSE' } if `pacfpoints' > 0 { return matrix PACF = `PACF' } } } end program runmlwin_mcmcdiagnostics, plugin ****************************************************************************** exit