*
* @PerronNgMTests  series  start  end
* computes one or more of the Perron-Ng "M" unit root tests.
*
* Parameters:
*   series     series to analyze. This should already be detrended.
*   start end  range of <<series>> to use (*not* range over which test is run)
*              By default, the defined range of <<series>>
*
* Options:
*   LAGS=number of additional lags [0]
*
*   GLS/[NOGLS]
*   DET=NONE/[CONSTANT]/TREND
*   CBAR=(positive) value for GLS rho=1-cbar/T. Default depends upon the model.
*   The DET and CBAR options are used in detrending the data using local-to-
*   unity detrending.
*
* These options allow you to retrieve any of the five variants
*
*   MZA=(output) modified Z-alpha statistic
*   MSB=(output) modified SB statistic
*   MZT=(output) modified Z-t statistic
*   ADF=(output) standard DF t-statistic
*   PT=(output)  point optimal Pt test
*
* References:
*   Perron and Ng(1996), "Useful Modifications to Some Unit Root Tests with
*    Dependent Errors and Their Local Asymptotic Properties", Review of Economic
*    Studies, vol 63, 435-463.
*
* Revision Schedule:
*   09/2009 Written by Tom Doan, Estima
*
source glsdetrend.src
*
procedure PerronNgMTests series start end
type series  series
type integer start end
*
option integer lags  0
option switch  gls   0
option choice  det   1  none constant trend
option real    cbar
*
option real    *mza
option real    *msb
option real    *mzt
option real    *adf
option real    *pt
*
local integer startl endl
local series  sdiff  xseries
local real    fiddle
local real    sumysq ssq ysq ssra ssr1 alphabar
*
if .not.%defined(series) {
   disp "Syntax: @PerronNgMTests(options)  SERIES start end"
   return
}

inquire(series=series) startl<<start endl<<end

if gls {
   @glsdetrend(det=det,cbar=cbar) series startl endl xseries
   compute alphabar=%rho,ssra=%rss
}
else
   set xseries startl endl = series

set sdiff  startl+1 endl = xseries-xseries{1}
linreg(noprint) sdiff startl+lags+1 endl
# xseries{1} sdiff{1 to lags}
if %defined(adf)
   compute adf=%tstats(1)
summarize(noprint)
# sdiff{1 to lags}
*
compute fiddle=(1.0-%sumlc)
compute ssq   =%rss/(%nobs*fiddle^2)
*
* We use the %NOBS from this instruction in the formulas
*
sstats startl+1 endl xseries{1}^2>>sumysq
compute sumysq=sumysq/%nobs^2
compute ysq   =xseries(endl)^2
*
if gls.and.%defined(pt) {
   @glsdetrend(det=det,cbar=0.0) series startl endl xseries
   compute ssr1=%rss
   compute pt=(ssra-alphabar*ssr1)/ssq
}

if %defined(mza)
   compute mza=(ysq/%nobs-ssq)/(2.0*sumysq)
if %defined(msb)
   compute msb=sqrt(sumysq/ssq)
if %defined(mzt)
   compute mzt=(ysq/%nobs-ssq)/sqrt(4.0*ssq*sumysq)
end
