*
*  @bayestst series  start end
*
*  Tests series for a unit root using the Bayesian procedure outlined in
*  Sims, "Bayesian Skepticism on Unit Root Econometrics", J. of Economic
*  Dynamics and Control, 1988 no. 2/3.
*
*  This generalizes Sims' formula to allow for intercept and trend.  These
*  however, are, in effect, given flat priors, which is not necessarily
*  realistic: the priors on constant, trend and rho should be correlated
*  since the three coefficients combined determine the level and trend of
*  the series.
*
*  The options are the following:
*
*     LAGS = number of total AR lags to be estimated
*     TREND/[NOTREND] : allow for a trending series
*     ALPHA=prior probability on stationary rhos [.8]
*     LIMIT=stationary prior concentrated on (LIMIT,1)
*
*  Output from the procedure is
*    The squared t-statistic for the unit root from the regression
*    The "Schwarz" asymptotic Bayesian limit
*    The small-sample limit (depends upon ALPHA and LIMIT)
*    The marginal ALPHA - with ALPHA set at this value (given LIMIT), the
*     posterior odds ratio is even.  A small value indicates that only a
*     very strong prior on the unit root will overcome the data evidence
*     against it.
*
*  Written by Tom Doan, November, 1989
*
procedure bayestst series start end
type series series
type integer start end
*
local series sdiff strend
local integer startl endl
local real teststat margalph alphtest
local real sigmarho compare1 compare2
*
option integer lags    1
option switch  trend   0
option real    alpha  .8
option real    limit  .5
*
inquire(series=series) startl<<start endl<<end
*
set sdiff startl+1 endl = series-series{1}
set strend = t
if lags>1 {
  if trend
     linreg(noprint) series startl+lags endl
     # series{1} constant strend sdiff{1 to lags-1}
  else
     linreg(noprint) series startl+lags endl
     # series{1} constant sdiff{1 to lags-1}
  }
else {
  if trend
     linreg(noprint) series startl+1 endl
     # series{1} constant strend
  else
     linreg(noprint) series startl+1 endl
     # series{1} constant
  }
*
compute sigmarho=%seesq*%xx(1,1)
compute teststat=(%beta(1)-1.0)/sqrt(sigmarho)
compute compare1=-log(sigmarho)
compute compare2=-log(sigmarho)+2*log((1-alpha)/alpha)-log(2*%pi)+2*log(1-limit)
compute alphtest=teststat**2+log(sigmarho)+log(2*%pi)-2*log(1-limit)
compute margalph=1./(1+exp(.5*alphtest))
display 'Bayesian Unit Root Test'
display ' Squared t     Schwarz Limit     Small Sample Limit   Marginal Alpha'
display #####.### teststat**2 @20 #####.### compare1 $
        @40 #####.### compare2 @60 #.#### margalph
end

