*
* @JohMLE( options)  start end
* # list of endogenous variables
*
* Computes Johansen lambda tests for cointegrating rank, and the ML
* estimator for a single cointegrating vector.
*
* Options:
* LAGS=#of lags in the VAR
* DETERM=NONE/[CONSTANT]/TREND - deterministic variables in the VAR (not in the CV)
* [PRINT]/NOPRINT
* CV=ML estimator for a single cointegrating vector [not used]
* VECTORS=matrix of (column) eigenvectors
* DUALVECTORS=matrix of eigenvectors of the dual problem
*
* Revision Schedule:
*  11/03 Written by Tom Doan, Estima
*  05/05 Revised to add output, return the ML estimator for a
*     single cointegrating vector.
*  01/06 Revised to improve output, add log likelihood to display, add VECTORS and
*     DUALVECTORS options, have procedure set %SIGMA
*  10/06 Flip TRACE column so it decrements from low rank to high
*  10/06 Add eigenvalues option
*  10/06 Have procedure define %s11, %s00 and %s01 - the components
*        of the covariance matrix with short run components swept out.
*
procedure JohMLE start end
type integer start end
*
option switch  print 1
option integer lags  1
option choice  determ 2 none constant trend
option vector  *cv
option rect    *vectors
option rect    *dualvectors
option vector  *eigenvalues
*
declare symm %%s00 %%s11
declare rect %%s01
*
local integer numvar i startl endl
local vect[int] v reglist
local vect[series] dv
local series trend
local rect s
local symm sss
local rect sf
local vect eigval
local rect eigvec
local vect a1
local vect[labels] vlabels
local real runlogl baselogl lambdamax
local string detdesc vs
*
* Pull in list of endogenous variables
*
enter(varying) v
compute numvar=%rows(v)
*
* Create differences of them
*
dim dv(numvar) vlabels(numvar)
do i=1,numvar
   set dv(i) = v(i){0}-v(i){1}
   compute vlabels(i)=%l(v(i))
end do i
*
set trend = t
*
* Generate list of required variables. Start with the
* current differences and lagged levels
*
dim reglist(0)
do i=1,numvar
   compute reglist=%rladdone(reglist,dv(i))
end do i
do i=1,numvar
   compute reglist=%rladdlag(reglist,v(i),1)
end do i
*
* Now add any required lags of the differences
*
if lags>1
   do i=1,numvar
      compute reglist=%rladdlaglist(reglist,dv(i),%seq(1,lags-1))
   end do i
*
* Add in the deterministic components
*
if determ>=2
   compute reglist=%rladdone(reglist,constant)
if determ>=3
   compute reglist=%rladdone(reglist,trend)
*
* Compute the cross product matrix and scale it to a covariance matrix
*
inquire(reglist) startl<<start endl<<end
# reglist
cmom startl endl
# reglist
compute s=%cmom*(1.0/%nobs)
*
* Sweep out the additional lags and the deterministic components.
*
compute s=%sweeplist(s,%seq(2*numvar+1,%ncmom))
*
compute %%s00=%xsubmat(s,       1,  numvar,       1,  numvar)
compute %%s11=%xsubmat(s,numvar+1,2*numvar,numvar+1,2*numvar)
compute %%s01=%xsubmat(s,       1,  numvar,numvar+1,2*numvar)
*
* Transform the problem so <<eigen>> can handle it, then
* backtransform the eigenvectors
*
compute sf=inv(%decomp(%%s11))
compute sss=sf*tr(%%s01)*inv(%%s00)*%%s01*tr(sf)
eigen sss eigval eigvec
compute eigvec=tr(sf)*eigvec
*
* Compute the log likelihood of the model w/o cointegration, that
* is, the VAR in first differences
*
compute %logl=%logconcdensity(%%s00,%nobs)
if determ==1
   compute detdesc=' No Deterministics'
else
if determ==2
   compute detdesc=' with Constant'
else
   compute detdesc=' with Constant + Trend'

report(action=define)
report(atrow=1,atcol=1,span) 'Likelihood Based Analysis of Cointegration'
*
* %logl is the log likelihood of a VAR with no cointegration
* baselogl starts as the log likelihood gap between the fully restricted
* and fully unrestricted VAR's, then adjusts in the loop below to show
* the LR trace statistics.
*
compute vs=''
compute baselogl=0
do i=1,numvar
   compute vs=vs+' '+%l(v(i))
   compute baselogl=baselogl-.5*%nobs*log(1-eigval(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) 'Data Points '+%nobs+' Lags '+lags+detdesc
report(atrow=6,atcol=1,span) 'Unrestricted eigenvalues and -T log(1-lambda)'
report(atrow=7,atcol=1) 'Rank' 'EigVal' 'Lambda-max' 'Trace' 'LogL'
compute runlogl=%logl
report(row=new,atcol=1) 0 '' '' '' runlogl
do i=1,numvar
   compute lambdamax=-%nobs*log(1-eigval(i))
   compute runlogl=runlogl+.5*lambdamax
   report(row=new,atcol=1) i eigval(i) lambdamax 2.0*baselogl runlogl
   compute baselogl=baselogl-.5*lambdamax
end do i
report(action=format,picture='*.###')
*
* Pull the first eigenvector
*
compute a1=%xcol(eigvec,1)
*
report(row=new)
report(row=new,atcol=1,span) 'Cointegrating Vector for Largest Eigenvalue'
report(row=new,atcol=1) vlabels
report(row=new,atcol=1) a1
if print
   report(action=show)
*
if %defined(cv)
   compute cv=a1
if %defined(vectors)
   compute vectors=eigvec
if %defined(eigenvalues)
   compute eigenvalues=eigval
*
compute %sigma=%%s00
if %defined(dualvectors) {
   compute sf=inv(%decomp(%%s00))
   compute sss=sf*%%s01*inv(%%s11)*tr(%%s01)*tr(sf)
   eigen sss eigval dualvectors
   compute dualvectors=tr(sf)*dualvectors
}
end


