*
* @SWDOLS(options)  start  end
* # list of endogenous variables
*
* Estimates cointegrating vectors using Stock and Watson's dynamic OLS.
*
* Parameters:
*   start end   Range over which to run the regression (by default, the maximum
*               range permitted by the variables and their leads and lags)
*
* Options:
*   rank=number of cointegrating vectors to estimate [1]
*     This estimates "rank" equations, with the first "rank" endogenous
*     variables used as the LHS variables of these, while the remaining
*     variables are the RHS variables.
*
*   lags=number of leads and lags of the differences of the RHS variables
*     to include (in addition to the current level and the current
*     difference) [1]
*
*   det=[constant]/trend
*     Deterministic components to include
*
*   ar=number of lags in the AR(VAR if rank>1) used in computing the
*     spectral density matrix at 0 frequency.
*
*   title=title for report ["Dynamic OLS Estimates"]
*
*   ect=(output) VECT[EQUATION] of equations describing the cointegrating vectors
*
* Example:
*
* @swdols(rank=2)
* # c in y
*
* will estimate cointegrating vectors of the forms c = betac * y and
* in = betai * y. (Rank=2, so the first two endogenous variables are used
* on the LHS and the remaining variable is on the RHS).
*
* Description:
*
* The dynamic regressions are estimated by OLS equation by equation with
* the RHS endogenous variable(s), a constant and current, leads and lags
* of the differenced RHS endogenous variable(s) as explanatory
* variables. The covariance matrix is computed as Su(0) Kroneker
* inv(X'X) where Su(0) is the spectral density at zero frequency of the
* residuals, and inv(X'X) is the inverse cross product of the explanatory
* variables. The output shows only the estimated cointegrating
* relationships, omitting the constant and lags and leads of the
* differenced variables.
*
* Variables Defined:
* %BETA and %XX are defined as they would be for a SUR estimate on the
* cointegrating vectors, that is, they are for a stacked system of
* equations. In the example above, %BETA(1) is the betac, and %BETA(2)
* is betai with %XX(1,1) the estimated variance of the betac
* coefficient. After using @SWDOLS, You can use TEST, RESTRICT and
* MRESTRICT to test hypotheses on the cointegrating vectors. For
* instance, in the example above, you could test betar=1 and betac=1
* (jointly) by
*
* TEST(TITLE="Test of Balanced Growth")
* # 1 2
* # 1.0 1.0
*
* Reference:
*   Stock and Watson(1993), "A Simple Estimator of Cointegrating Vectors in Higher
*    Order Integrated Systems", Econometrica, vol 61, 783-820.
*
* Revision Schedule:
*   09/2005 Written by Tom Doan, Estima
*   05/2011 Added TITLE and DET options
*   10/2014 Add ECT option
*
procedure SWDOLS start end
type integer start end
*
option integer       rank 1
option integer       lags 1
option integer       ar   1
option string        title
option choice        det  1  constant trend
option vect[equ]     *ect
*
local vect[int]     vlist
local integer       i j n hrows rhs
local integer       startl endl
local vect[series]  dy u
local vect[integer] reglist reglistu
local equation      tempeq baseeq
local rect          cmomx cmomu
local rect[integer] temptbl copytbl
local rect          sw
local string        ltitle
local report        dreport
local series        trend
local label         label
*
* Get the list of endogenous variables
*
enter(varying) vlist
compute n=%rows(vlist)
*
* Figure out the defined range of the endogenous variables
*
inquire(reglist) startl endl
# vlist
*
* Allow for lags to get the target regression range
*
compute startl=startl+lags+1,endl=endl-lags
*
* Take the input range if given, but clip it to the range
* permitted by the data.
*
if %defined(start)
   compute startl=%imax(start,startl)
if %defined(end)
   compute endl  =%imin(end  ,endl)
*
* Generate the differences of the variables which will be on the RHS
*
dim dy(n-rank) u(rank)
do i=1,n-rank
   diff vlist(rank+i) startl-lags endl+lags dy(i)
end do i
*
* Build the regressor list, starting with the endogenous variables, then
* adding the constant and trend (if needed), and the lags and leads of
* the differenced RHS endogenous variables.
*
compute reglist=vlist
compute reglist=%rladdone(reglist,constant)
if det==2 {
   set trend = t
   compute reglist=%rladdone(reglist,trend)
}
*
* Save the equation form using only the endogenous variables and the
* deterministics.
*
equation baseeq
# reglist
*
do i=1,n-rank
   compute reglist=%rladdlaglist(reglist,dy(i),%seq(-lags,lags))
end do i
*
* Compute the cross product matrix, sweep out all RHS variables
*
cmom startl endl
# reglist
compute cmomx=%sweeplist(%cmom,%seq(rank+1,%ncmom))
*
* This takes the regressor list just used and removes the LHS endogenous
* variables. (It's a bit clumsy because the submatrix functions don't
* work on integer matrices).
*
equation tempeq
# reglist
compute temptbl=%eqntable(tempeq)
dim copytbl(2,%cols(temptbl)-rank)
ewise copytbl(i,j)=temptbl(i,j+rank)
equation tempeq
# %rlfromtable(copytbl)
*
* Generate the residuals from the DOLS regression. %eqnxvector(tempeq,t)
* pulls the RHS variables at time t, and column i, rows rank+1 to %ncmom
* has the regression coefficients.
*
do i=1,rank
   set u(i) startl endl = vlist(i){0}-$
      %dot(%eqnxvector(tempeq,t),%xsubmat(cmomx,rank+1,%ncmom,i,i))
end do i
*
* Build a regressor list for the VAR on the residuals. We have to do
* this a lag at a time because we need to pull out the matrices of lag
* coefficients in order to recolor the spectral density.
*
dim reglistu(0)
do i=1,rank
   compute reglistu=%rladdone(reglistu,u(i))
end do i
do j=1,ar
   do i=1,rank
      compute reglistu=%rladdlag(reglistu,u(i),j)
   end do i
end do j
*
* Compute the cross product matrix, scale through by the degrees of
* freedom and sweep out the lags of the residuals.
*
cmom startl+ar endl
# reglistu
compute cmomu=%cmom/(%nobs-%eqnsize(tempeq))
compute cmomu=%sweeplist(cmomu,%seq(rank+1,%ncmom))
*
* Compute the Su(0) matrix
*
compute sw=%identity(rank)
do j=1,ar
   compute sw=sw+%xsubmat(cmomu,1,rank,j*rank+1,(j+1)*rank)
end do j
compute sw=%mqform(%xsubmat(cmomu,1,rank,1,rank),tr(inv(sw)))
*
if %defined(title)
   compute ltitle=title
else
   compute ltitle="Dynamic OLS Estimates"
report(use=dreport,action=define,title=ltitle)
report(use=dreport,atrow=1,atcol=1,span) ltitle
report(use=dreport,atrow=2,atcol=1,span) "Estimated from "+%datelabel(startl)+" to "+%datelabel(endl)
report(use=dreport,atrow=3,atcol=1,span) "Lags/Leads Used in Dynamic Regressions: "+lags
report(use=dreport,atrow=4,atcol=1,span) "AR Lags in Computing S(0): "+ar
if det==2
   report(use=dreport,atrow=5,atcol=1,span) "Using Constant + Trend as deterministic components"
else
   report(use=dreport,atrow=5,atcol=1,span) "Using Constant as deterministic component"
*
compute hrows=5
*
* LHS Endogenous variables across the top
*
do i=1,rank
   report(use=dreport,atrow=hrows+1,atcol=i+1,align=center) %l(vlist(i))
end do i
*
* RHS variables. These are the leftover endogenous variables.
*
compute rhs=n-rank
do j=1,rhs
   compute label=%l(vlist(rank+j))
   report(use=dreport,atrow=hrows+2*j,atcol=1,align=left) label
   do i=1,rank
      report(use=dreport,atrow=hrows+2*j,atcol=i+1) cmomx(rank+j,i)
      report(use=dreport,atrow=hrows+2*j+1,atcol=i+1,special=parens) sqrt(sw(i,i)*cmomx(rank+j,rank+j))
   end do i
end do j
report(use=dreport,action=format,align=decimal,width=8,atrow=hrows+2,atcol=2)
report(use=dreport,action=show)
*
* Create a fake regression so hypothesis tests can be done on the output
*
dim reglist(0)
do i=1,rank
   do j=rank+1,n
      compute reglist=%rladdone(reglist,vlist(j))
   end do j
end do i
*
linreg(create,noscale,coeffs=%vec(%xsubmat(cmomx,rank+1,n,1,rank)),$
  covmat=%kroneker(sw,%xsubmat(cmomx,rank+1,n,rank+1,n)),noprint) vlist(1) startl endl
# reglist
*
* Define ECT's
*
if %defined(ect) {
   dim ect(rank)
   do i=1,rank
      equation(coeffs=%zeros(rank,1)~~%xsubmat(cmomx,rank+1,%eqnsize(baseeq),i,i)) $
         ect(i) vlist(i)
      # %rlfromtable(%eqntable(baseeq))
   end do i
}
end
