*
* @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]
*
* Revision Schedule:
*  11/03 Written by Tom Doan, Estima
*  05/05 Revised to add output, return the ML estimator for a
*     single cointegrating vector.
*
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
*
local integer numvar i
local vect[int] v reglist
local vect[series] dv
local series trend
local rect s
local symm suu svv
local rect suv sss
local vect eigval
local rect eigvec
local vect a1
local vect[labels] vlabels
*
* 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
*
cmom start end
# 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 suu=%xsubmat(s,       1,  numvar,       1,  numvar)
compute svv=%xsubmat(s,numvar+1,2*numvar,numvar+1,2*numvar)
compute suv=%xsubmat(s,       1,  numvar,numvar+1,2*numvar)
*
compute sss=inv(svv)*tr(suv)*inv(suu)*suv
eigen sss eigval eigvec
report(action=define)
report(atrow=1,atcol=1,span) 'Likelihood Based Analysis of Cointegration'
report(atrow=2,atcol=1,span) 'Data Points' %nobs 'Lags' lags
report(atrow=4,atcol=1,span) 'Unrestricted eigenvalues and -T log(1-lambda)'
report(atrow=5,atcol=1) 'Rank' 'EigVal' 'Lambda-max'
 do i=1,numvar
   report(row=new,atcol=1) i eigval(i) -%nobs*log(1-eigval(i))
end do i
report(action=format,picture='*.###')
*
*compute sssx=inv(suu)*suv*inv(svv)*tr(suv)
*eigen sssx eigx betax
*disp eigx
*disp betax
*disp eigval
*disp eigvec
*
* Pull the first eigenvector
*
compute a1=%xcol(eigvec,1)
*
* This normalizes a1 so that tr(a1)*Svv*a1=1
*
compute a1=a1*(1.0/sqrt(%qform(svv,a1)))
*
* And this normalizes a1 to have a unit coefficient in the first element
*
compute a1=a1*(1.0/a1(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
end




