*
* @paneldols( options ) start end
* # list of variables (list dependent variable first)
*
* Procedure for multivariate group mean panel DOLS tests in heterogeneous panels.
*
* Parameters:
*   start,end  = Range of entries to use. Defaults to all available data.
*                If you set these, do so as follows:
*                 start = first observation of the first individual in the desired range
*                 end   = last observation of the last individual in the desired range
*
* Options:
*   block = time dimension [by default is taken from CALENDAR(PANELOBS)]
*           (see note below for unbalanced panels)
*   tdum/[notdum]   Use TDUM to subtract out common time effects
*   bvec =VECTOR with hypothesized slope coefficients [all zeros]. The t-statistics
*     reported are for tests against this vector.
*   print=none/short/[full]
*     Sets the level of printed output. PRINT=FULL includes the estimates for the
*     individuals, while PRINT=SHORT just does the full sample estimates.
*   average=[simple]/sqrt/precision
*     Determines how the individual estimates are combined to compute the full sample
*     estimate. AVERAGE=SIMPLE takes a simple arithmetic average. This is the default,
*     and the behavior of Pedroni's original program. AVERAGE=SQRT weights each
*     individual by the diagonal matrix formed by taking the square roots of the
*     precision matrix (inverse covariance matrix) of the estimates for that individual.
*     This matches up with the averaging done in computing the t-statistics, that is,
*     the coefficients and covariance matrix from AVERAGE=SQRT will reproduce the average
*     t-statistics. AVERAGE=PRECISION weights each individual by the precision of its
*     estimates.
*
*   ibetas  = (output) RECT of estimated individual betas.
*             Dimension number of individuals x number of RHS variables.
*   itstats = (output) RECT of estimated individual t-statistics.
*             Dimension number of individuals x number of RHS variables.
*
*   dlags = number of leads and lags on the differenced RHS variables [2]
*   lags  = number of lags for kernel estimator for variance [5]
*
* For unbalanced panels, you should enter into your data set either blanks or the value
* NA for missing values, so that the end result is that each member has uniform length T,
* and the data is stacked into a column of length N*T.
*
* References:
*
*   Pedroni, Peter (2001) "Purchasing Power Parity Tests in Cointegrated Panels,"
*   Review of Economics and Statistics, 83, 727-731.
*
* Revision Schedule:
*   01/2000 Written by Peter Pedroni
*   04/2008 Converted to procedure by Tom Doan, Estima. PRINT and AVERAGE options added.
*   11/2009 Fixed a problem when the first n individuals were skipped. Solution is to compute
*             first and last individuals and use those in DO I and DO J loops, rather than "1,nsecs".
*
procedure panelDOLS  startp  endp
type integer startp endp
*
option integer block
option switch  tdum     0
option integer dlags    2
option integer lags     5
option vector  bveco
option choice  average  1   simple sqrt   precision
option choice  print    3   none   short  full
option rect    *ibetas
option rect    *itstats
*
local integer     Nsecs
local integer     Tperiods jperiods
local integer     m
local vect[int]   datavec rhs
local vec[label]  labelvec
local vec         bvec
local integer     unbal
local real        mint maxt

local vec[series] xvec dxvec dvec rvec ytilde
local integer     start end
local series      u tvec dolsres
local vec[int]    savecal
local integer     i j k row
local vect[int]   reglist
local real        idlrvar
local vect        pdlbsum1 pdlbsum2 pdlbsum3 pdltsum bdols tdols
local symm        pdlvsum1 pdlvsum2 pdlvsum3 pdlhsum2 pdlhsum3 vdols hdols ddols
local series[vec] betas tstats
local integer     firstInd LastInd
*
* Get list of variables (LHS variable first)
*
enter(varying) datavec
compute m=%rows(datavec)-1
*
* Create list to be used for RHS (that is, datavec w/o the first entry, which will
* be used for the dependent variable).
*
dim rhs(m)
ewise rhs(i)=datavec(i+1)
dim labelvec(m+1)
ewise labelvec(i)=%l(datavec(i))

*
*  Get the block size
*

if .not.%defined(block)
   compute tperiods=2//1-1//1
else
   compute tperiods=block

*
* Get the data range
*

inquire(reglist) start<<startp end<<endp
# datavec

*
* Normalize <<start>> to begin at the first entry in an individual's record, and
* <<end>> to end at the last entry.
*

compute start=(start-1)/tperiods*tperiods+1
compute end  =((end-1)/tperiods+1)*tperiods
compute nsecs=(end-start+1)/tperiods
compute firstInd = (start/tperiods)+1
compute lastInd = firstInd+(nsecs-1)

*
* Run an OLS mainly to figure out what data points are available.
*

linreg(noprint) datavec(1) start end u
# rhs

*
* See if the data set is unbalanced. This also gets the count of the number of
* observations per individual into <<tvec>>.
*

compute savecal = %calendar()
calendar(panelobs=tperiods)
panel(icount=1,compress) u start end tvec
sstats(maximum) 1 nsecs tvec(t)>>maxt
sstats(minimum) 1 nsecs tvec(t)>>mint
compute unbal = fix(maxt<>mint)
calendar(recall=savecal)

*
* Set up the hypothesized beta vector
*

if .not.%defined(bveco)
   compute bvec=%zeros(m,1)
else
   compute bvec=bveco
*
* Dimension everything
*

dim ytilde(m+1) dvec(m+1) rvec(m+1) xvec(m) dxvec(m) tdols(m)
*
*  Copy either the raw data, or (if tdum) the deviations from time period means
*  means into ytilde.
*
do K=1,m+1
   if tdum==1
      panel(block=tperiods,entry=1.0,time=-1.0) datavec(K) start end ytilde(k)
   else
      set ytilde(k) start end = datavec(k){0}
end do k
*
* Build the regressor list for the DOLS regression, starting with the RHS
* endogenous variables (with any prior adjustments), then adding the constant and
* the lags and leads of the differenced RHS endogenous variables.
*
dim reglist(0)
do i=1,m
   compute reglist=%rladdone(reglist,xvec(i))
end do i
compute reglist=%rladdone(reglist,constant)
do i=1,m
   compute reglist=%rladdlaglist(reglist,dxvec(i),%seq(-dlags,dlags))
end do i

*
* Zero out the matrices for sums across individuals
*

compute pdlbsum1 = %zeros(m,1)
compute pdlbsum2 = %zeros(m,1)
compute pdlbsum3 = %zeros(m,1)
compute pdltsum  = %zeros(m,1)
compute pdlvsum1 = %zeros(m,m)
compute pdlvsum2 = %zeros(m,m)
compute pdlhsum2 = %zeros(m,m)
compute pdlhsum3 = %zeros(m,m)

infobox(action=define,progress,lower=1,upper=NSecs) "Computing Panel Statistics"
do J=firstInd,lastInd
   *
   * Copy the information into the first entries. SAMPLE with the SMPL option
   * is used to ensure that we get the same block of data for each series in
   * the individual's record. Difference and de-mean the data.
   *
   do K=1,M+1
      sample(smpl=%valid(u)) ytilde(k) (j-1)*tperiods+1 j*tperiods dvec(k) 1
      compute jperiods=%nobs
      diff dvec(K) 2 jperiods rvec(K)
      diff(center) rvec(K) 2 jperiods rvec(K)
   end do K
   *
   * Copy the subvector of data and differenced data for the RHS variables
   *
   do K=1,M
      set xvec(K)  1 jperiods = dvec(K+1)
      set dxvec(K) 2 jperiods = rvec(K+1)
   end do K
   *
   * Do DOLS augmented regression
   *
   linreg(noprint) dvec(1) 2 jperiods dolsres
   # reglist
   *
   * Compute the long-run variance of the residuals.
   *
   mcov(lwindow=bartlett,lags=lags,noprint) 2 jperiods
   # dolsres
   compute idlrvar = %cmom(1,1)/(jperiods-1)
   *
   * Pull out of this regression:
   *  1. The coefficients on the RHS endogenous variables
   *  2. The covariance matrix of those coefficients
   *  3. The precision matrix of those coefficients (inv of #2)
   *  4. The matrix of reciprocal square roots of the diagonal elements of the
   *     covariance matrix. That's the weighting matrix for the coefficients which
   *     will match up with the averaged t-statistics.
   *
   compute bdols=%xsubvec(%beta,1,m)
   compute vdols=%xsubmat(%xx,1,m,1,m)*idlrvar
   compute hdols=inv(vdols)
   compute ddols=inv(%sqrt(%diag(%xdiag(vdols))))
   *
   * Compute t-stat vs the hypothesized vector
   *
   ewise tdols(i)  = (bdols(i)-bvec(i))*ddols(i,i)
   *
   * Keep track of several combined measures:
   *  pdlbsum1 is the arithmetic average
   *  pdlbsum2 is weighted by the square root of the diagonals of the precision
   *  pdlbsum3 is weighted by the precision
   *  pdlvsum1 is the sum of the variance matrices
   *  pdlvsum2 is the sum of the #2 rescaled variance matrices
   *  pdlhsum2 is the sum of the #2 weights
   *  pdlhsum3 is the sum of the #3 weights
   *  pdltsum  is the sum of the t-statistics
   *
   compute pdlbsum1 = pdlbsum1  + bdols
   compute pdlbsum2 = pdlbsum2  + ddols * bdols
   compute pdlbsum3 = pdlbsum3  + hdols * bdols
   compute pdlvsum1 = pdlvsum1  + vdols
   compute pdlvsum2 = pdlvsum2  + ddols * vdols * ddols
   compute pdlhsum2 = pdlhsum2  + ddols
   compute pdlhsum3 = pdlhsum3  + hdols
   compute pdltsum  = pdltsum   + tdols
   gset betas  j j = bdols
   gset tstats j j = tdols
   infobox(action=modify,current=j)
end do J
infobox(action=remove)
*
* Standardize the estimates as appropriate
*
compute pdltsum  = (1.0/sqrt(NSecs))*pdltsum
compute pdlbsum1 = (1.0/NSecs)*pdlbsum1
compute pdlbsum2 = inv(pdlhsum2)*pdlbsum2
compute pdlbsum3 = inv(pdlhsum3)*pdlbsum3
compute pdlvsum1 = (1.0/NSecs)**2*pdlvsum1
compute pdlvsum2 = inv(pdlhsum2)*pdlvsum2*inv(pdlhsum2)
compute pdlvsum3 = inv(pdlhsum3)
*
* Define variables
*
compute %tstats=pdltsum
if average==1
   compute %beta=pdlbsum1,%xx=pdlvsum1
else
if average==2
   compute %beta=pdlbsum2,%xx=pdlvsum2
else
if average==3
   compute %beta=pdlbsum3,%xx=pdlvsum3

    *** DISPLAY INDIVIDUAL AND PANEL RESULTS ***

report(action=define)
report(atrow=1,atcol=1,span) "Mean Group Panel DOLS Estimation"
report(atrow=2,atcol=2) "LHS Variable" labelvec(1)
report(atrow=3,atcol=2) "Individuals" NSecs
if unbal==0
   report(atrow=4,atcol=2) "Time Periods" TPeriods
else
   report(atrow=4,atcol=2) "Time Periods" "Unbalanced"
report(atrow=5,atcol=2) "Dynamic Lags" dlags
if tdum == 0 ; { ; report(atrow=6,atcol=2,span) "common time dummies NOT included" ; }
if tdum == 1 ; { ; report(atrow=6,atcol=2,span) "common time dummies included" ; }
report(atrow=8,atcol=2) "RHS Variable" "H0 Coefficient"
do k=1,m
   report(atrow=8+k,atcol=2) labelvec(k+1) bvec(k)
end do k
compute row=8+m+2
report(atrow=row,atcol=1) "Member" "Variable" "Coefficient" "t-Statistic"
if print==3 {
   do i=firstInd,lastInd
      do k=1,m
         if k==1
            report(row=new,atcol=1) "No."+i labelvec(k+1) betas(i)(k) tstats(i)(k)
         else
            report(row=new,atcol=2)         labelvec(k+1) betas(i)(k) tstats(i)(k)
      end do k
   end do i
   report(row=new,atcol=1,span) "*************"
}
do k=1,m
   if k==1
      report(row=new,atcol=1) "Group" labelvec(k+1) %beta(k) pdltsum(k)
   else
      report(row=new,atcol=2)         labelvec(k+1) %beta(k) pdltsum(k)
end do k
if print>1
   report(action=show)

if %defined(ibetas) {
   dim ibetas(nsecs,m)
   ewise ibetas(i,k)=betas(i)(k)
}
if %defined(itstats) {
   dim itstats(nsecs,m)
   ewise itstats(i,k)=tstats(i)(k)
}
end


