*
* @SWTrends( options )  start end
* # list of endogenous variables
*
* Test for cointegration rank via Stock-Watson common trends approach.
*
* Parameters:
*   start end   Range of data to use. By default, the maximum range permitted
*               by the series involved.
*
* Options:
*   LAGS=#of lags in the VAR
*   NULL=# of unit roots under the null
*   DET=NONE/[CONSTANT]/TREND deterministic parts to be removed before analysis
*
* Reference:
*   Stock and Watson(1988), "Testing for Common Trends", JASA, vol. 83,
*    1097-1107.
*
* Revision Schedule:
*  09/2005 Written by Tom Doan, Estima
*
procedure SWTrends start end
type integer start end
*
option switch  print 1
option integer lags  1
option choice  determ 2 none constant trend
option integer null 1
*
local integer i j k n startl endl
local vect[int] v reglist reglistx
local series trend
local rect cmomx eigvec gamma
local integer ntotal
local vector eigval
local vect[complex] cv
local rect[complex]  cvect
local vect[series] w dw dt
local real teststat
local integer notim
local string vs detdesc
*
* Pull in list of endogenous variables
*
enter(varying) v
compute n=%rows(v)
*
* Figure out the range of data to use
*
inquire(reglist) startl<<start endl<<end
# v
*
* Detrend or demean the data as necessary
*
dim dt(n)
set trend startl endl = t
do i=1,n
   linreg(entries=determ-1,noprint) v(i) startl endl dt(i)
   # constant trend
end do i
compute notim=%nobs
*
dim w(null) dw(null)
*
* The original paper indicates that the cointegrating vectors are
* estimated under the null by taking the principal components. The
* source code, however, uses a different calculation. This doesn't
* matter if the null is no cointegration (NULL=dim of x). I'm not
* convinced that this next calculation is correct, and it definitely
* doesn't work well with the way that RATS normalizes eigenvectors.
*
dim reglistx(0)
do i=1,n
   compute reglistx=%rladdone(reglistx,dt(i))
end do i
do i=1,n
   compute reglistx=%rladdlag(reglistx,dt(i),1)
end do i
cmom start end
# reglistx
compute cmomx=%sweeplist(%cmom,%seq(n+1,2*n))
eigen(cval=cv,cvect=cvect) %xsubmat(cmomx,n+1,2*n,1,n)
dim eigvec(n,n)
*
* Rearrange the eigenvectors so the largest roots are at the end
*
ewise eigvec(i,j)=%real(cvect(i,n+1-j))+%imag(cvect(i,n+1-j))
*
* This does the principal components method. Done this way, the
* eigenvectors that would be omitted (if null<dim(x)) have the smallest
* variance, and are thus (apparently) the ones corresponding to
* cointegrating vectors.
*
cmom startl endl
# v
eigen %cmom eigval eigvec
*
* Extract the linear combinations which are assumed (under the null) to
* be non-stationary.
*
do i=1,null
   set w(i)  startl endl = %dot(%xt(dt,t),%xcol(eigvec,i))
   set dw(i) startl endl = w(i)
end do i
*
* Compute a one lag VAR on the "W" series (the copy of it in DW). This
* is for the initial filtering.
*
dim reglistx(0)
do i=1,null
   compute reglistx=%rladdone(reglistx,dw(i))
end do i
do i=1,null
   compute reglistx=%rladdlag(reglistx,dw(i),1)
end do i
*
cmom start end
# reglistx
compute cmomx=%sweeplist(%cmom,%seq(null+1,2*null))
*
* Filter "w" to make "dw". (This could also have been done by just first
* differencing w; that would also give a valid test).
*
do i=1,null
   set dw(i) startl+1 endl = w(i)-%dot(%xt(w,t-1),%xsubmat(cmomx,null+1,2*null,i,i))
end do i
*
* Now set up a VAR of dw on its lags to compute the whitening filter
*
dim reglist(0)
do i=1,null
   compute reglist=%rladdone(reglist,dw(i))
end do i
*
* Now add any required lags. These are grouped by lags to make it easier
* to use the sweep operator to do the regressions.
*
do j=1,lags
   do i=1,null
      compute reglist=%rladdlag(reglist,dw(i),j)
   end do i
end do j
compute ntotal=null*(lags+1)
*
* Compute the cross product matrix
*
cmom startl+lags+1 end
# reglist
compute cmomx=%cmom/%nobs
*
* Regress out the lags
*
compute cmomx=%sweeplist(cmomx,%seq(null+1,ntotal))
*
* Filter the w series by PI(L) to make a new DW.
*
do i=1,null
   set dw(i) startl+lags endl = w(i)
   do j=1,lags
      do k=1,null
         set dw(i) startl+lags endl = dw(i)-cmomx(j*null+k,i)*w(k){j}
      end do k
   end do j
end do i
*
* Finally, compute the one lag VAR on the filtered data and
* examine the eigenvalues of the coefficient matrix.
*
cmom startl+lags+1 endl
# reglistx
compute cmomx=%sweeplist(%cmom,%seq(null+1,2*null))
compute gamma=%xsubmat(cmomx,null+1,2*null,1,null)
eigen(cvalues=cv) gamma
*
* The Stock-Watson procedure uses as "T" the total number
* of observations, not the number used in the final calculation,
* so we do that as well.
*
if determ==1
   compute detdesc=" None"
else
if determ==2
   compute detdesc=" Constant"
else
   compute detdesc=" Constant + Trend"

report(action=define)
report(atrow=1,atcol=1,span) "Stock-Watson Test of Number of Common Trends"
compute vs=""
do i=1,n
   compute vs=vs+" "+%l(v(i))
end do i
report(atrow=2,atcol=1,span) "Variables: "+vs
report(atrow=3,atcol=1,span) "Estimated from "+%datelabel(startl)+" to "+%datelabel(endl)
report(atrow=4,atcol=1,span) "Deterministic Components Removed:"+detdesc
report(atrow=5,atcol=1,span) "Lags Used in Prewhitening Filter "+lags
*
report(atrow=6,atcol=1,align=center) "Null" "Versus" "Statistic"
do i=null-1,0,-1
   compute teststat=notim*(%real(cv(i+1))-1)
   report(row=new,atcol=1) null i teststat
end do i
report(action=show)
end


