*
* @MCMCPostProc(NDRAWS=# of draws,other options) stats
* Post processor for Markov Chain Monte Carlo statistics
*
* Parameters:
*  stats     SERIES[VECT] of saved statistics
*
* Options:
*  NDRAWS=(input, REQUIRED) # of draws (length of stats)
*  MEAN=(output) [VECTOR] averages of each component
*  STDERRS=(output) [VECTORS] standard errors of each component
*  CV=(output) [SYMMETRIC] estimated covariance matrix
*  NSE=(output) [VECTOR] numerical standard errors of each component
*
*  PERCENTILES=(input) ||desired percentiles|| [not used]
*  QUANTILES=(output) VECT[VECTOR] with percentiles
*
*  CD=[VECTOR] Geweke CD measures
*  BW=[VECTOR] Between-within measures
*  BLOCKS=number of blocks for BW [10]
*
* Revision Schedule:
*  05/2009 Written by Tom Doan, Estima
*  07/2009 Revised to cap the number of lags on the Bartlett window
*  05/2014 Add BW option
*  05/2017 Change lags for calcuation of NSE's and CD measures.
*
procedure MCMCPostProc stats
type series[vect] stats
option integer    ndraws
option vector     *mean
option vector     *stderrs
option symm       *cv
option vector     *nse
option vector     percentiles
option vec[vec]   *quantiles
option vector     *bw
option vector     *cd
option integer    blocks  10
*
local integer i j n lags
local symm    vpost
local vect    bpost
local vect    tquant
local series  u
local symm    t_s0
local integer l1 l2 u1 u2
local real    b1 b2 nv1 nv2
local real    tss wss
local integer tdf wdf
*
if .not.%defined(stats).or..not.%defined(ndraws) {
   disp "###@MCMCPostProcessor(NDRAWS=# of draws,other options)  SERIES[VECT] of statistics"
   return
}
compute n=%rows(stats(1))
dim vpost(n,n) bpost(n)
*
* Compute the means
*
do i=1,n
   sstats(mean) 1 ndraws stats(t)(i)>>bpost(i)
end do i
*
* Given the means, compute the covariance matrix
*
compute vpost=%zeros(n,n)
do t=1,ndraws
   compute vpost=vpost+%outerxx(stats(t)-bpost)
end do t
compute vpost=vpost/ndraws
*
* If requested, save the vector of means and the covariance matrix
*
if %defined(mean)
   compute mean=bpost
if %defined(cv)
   compute cv=vpost
*
* If requested, save the standard errors
*
if %defined(stderrs)
   compute stderrs=%sqrt(%xdiag(vpost))
*
* This is the time-consuming calculation. Leave if we don't need it.
*
if .not.(%defined(nse).or.%defined(cd).or.%defined(bw).or.%defined(percentiles))
   return

if %defined(nse)
   dim nse(n)
if %defined(cd)
   dim cd(n)
if %defined(bw)
   dim bw(n)
if %defined(quantiles)
   dim quantiles(n)

do i=1,n
   *
   * Compute the LR variance over the full set of draws. Use that to
   * compute the numerical standard error.
   *
   set u 1 ndraws = stats(t)(i)-bpost(i)
   compute lags=%imin(fix(2.0*ndraws^(.4)),400)
   mcov(lags=lags,lwindow=bartlett,matrix=t_s0) 1 ndraws u
   # constant
   if %defined(nse)
      compute nse(i)=sqrt(t_s0(1,1)/ndraws^2)
   *
   if .not.%defined(cd)
      goto skipcd
   *
   * Calculate CD
   *
   * Compute the sample mean and LR variance over the early subsample
   *
   compute l1=1,u1=fix(.1*ndraws)
   compute u2=ndraws,l2=u2-fix(.4*ndraws)+1
   *
   sstats(mean) l1 u1 stats(t)(i)>>b1
   set u l1 u1 = stats(t)(i)-b1
   mcov(lags=fix(lags*.1^.4),lwindow=bartlett,matrix=t_s0) l1 u1 u
   # constant
   compute nv1=t_s0(1,1)/%nobs^2
   *
   * And over the late subsample
   *
   sstats(mean) l2 u2 stats(t)(i)>>b2
   set u l2 u2 = stats(t)(i)-b2
   mcov(lags=fix(lags*.4^.4),lwindow=bartlett,matrix=t_s0) l2 u2 u
   # constant
   compute nv2=t_s0(1,1)/%nobs^2
   *
   * Convert into Z-score on the difference
   *
   compute cd(i)=(b1-b2)/sqrt(nv1+nv2)
   *

:skipcd
   *
   if .not.%defined(bw)
      goto skipbw

   stats(noprint) u 1 ndraws
   compute tss=%variance*(%nobs-1),tdf=%nobs-1
   compute wss=0.0,wdf=0
   do j=1,blocks
      compute l1=(j-1)*ndraws/blocks+1,l2=j*ndraws/blocks
      stats(noprint) u l1 l2
      compute wss=wss+%variance*(%nobs-1)
      compute wdf=wdf+(%nobs-1)
   end do j
   compute bw(i)=((tss-wss)/(tdf-wdf))/(wss/wdf)

:skipbw
   compute tquant=%fractiles(u,percentiles)
   if %defined(quantiles)
      compute quantiles(i)=tquant+bpost(i)
end do i
end MCMCPostProcessor
