*
* @paneldols( options ) start end
* # list of variables (list dependent variable first)
*
* Procedure for multivariate group mean panel DOLS estimation in
* heterogeneous panels. Note that this can very quickly exhaust the
* degrees of freedom in a data set. There are 2*DLAGS+1 (by default 5)
* added regressors in the differences for each right side endogenous
* variable plus you lose DLAGS*2+1 data points allowing for lags and
* leads and differences. So 20 observations per individual, two right
* side endogenous variables, DLAGS=2 leaves you with 15 usable
* observations, and 13 regressors (constant + 2 current RHS + 2 x 5
* additional lags and leads on the differences).
*
* Parameters:
*   start,end  = Range of entries to use. Defaults to all available data.
*
* Options:
*   DET=NONE/[CONSTANT]/TREND
*   TDUM/[NOTDUM]   Use TDUM to subtract out common time effects
*   LAGS=number of lags to use in the Bartlett kernel (for computing the variance) [5]
*   DLAGS=number of lags and leads on the differences [2]
*   SMPL=Dummy series with 0's in entries to skip [include all]
*
*   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.
*   BVEC=VECTOR with hypothesized slope coefficients [all zeros]. The
*     t-statistics reported are for tests against this vector.
*   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.
*
* Maintained for compatibility. Use the %%IBETAS, %%ITSTATS and
* %%ISTDERRS instead.
*
*   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.
*
* Don't use. Use a PANEL-dated CALENDAR instead. For unbalanced panels,
* you should enter into your data set either blanks or the value NA for
* missing values to pad out the individual records.
*
*   BLOCK  = time dimension
*
* Variables Defined:
*
*   %BETA        group estimates of coefficients
*   %STDERRS     group estimates of standard errors
*   %TSTATS      group estimates of t-statistics (for testing BVEC if you include it).
*   %XX          estimates of covariance matrix of group estimates
*   %RESIDS      cointegrating residuals (computed at the group mean coefficients)
*
* Except for %%IXX, these are all # of regressors x # of individuals
* matrices. The regressors are the endogenous variables first, followed
* by intercept and (if needed) trend. The individuals direction includes
* only the individuals actually in the sample.
*
*   %%IBETAS     individual coefficients
*   %%ISTDERRS   individual standard errors
*   %%ITSTATS    individual t-statistics (for testing BVEC if you include it).
*   %%IXX        individual covariance matrices. This is a VECTOR[SYMM] with
*                %%IXX(i) as the covariance matrix for individual i.
*
* 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.
*   04/2012 Extensive changes by Tom Doan, Estima
*           Added SMPL option
*           Added standard DET option, including allowance for TREND.
*           Change to save full coefficient vector and covariance matrices,
*           including the determinstics.
*
procedure panelDOLS  start  end
type integer start end
*
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 choice  det      2   none   constant trend
option series  smpl
option string  title
*
option integer block
option rect    *ibetas
option rect    *itstats
*
local integer     Nsecs
local integer     Tperiods
local integer     m
local vect[int]   datavec rhs
local vec[label]  labelvec
local vec         bvec
local integer     unbal

local vec[series] xvec dxvec dvec rvec ytilde
local series      u dolsres
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 integer     starti endi startl endl tmax tend isec
local series      time icount
local vect[int]   detvars
local string      descript ltitle
local integer     regnobs baseregs
local real        pt nt
local report      preport
local string      ltitle
local equation    baseeq
*
* Globals
*
dec rect %%ibetas %%istderrs %%itstats
dec vect[symm] %%ixx
*
* Check for valid settings
*
if %defined(block) {
   disp "###@PANELDOLS(BLOCK=nnn). No longer supported. Use PANEL calendar scheme"
   return
}
if det==1 {
   disp "###@PANELDOLS doesn't support DET=NONE"
   return
}
set time = %period(t)
*
* Figure out which deterministics are requested.
*
if det==3
   compute detvars=||constant,time||
else
   compute detvars=||constant||
*
if %size(detvars)==1
   compute descript="Constant"
else
   compute descript="Constant+Trend"
*
* Pull in list of series. Set m to the number of variables - 1
*
enter(varying) datavec
compute m=%rows(datavec)-1
dim labelvec(m+1)
ewise labelvec(i)=%l(datavec(i))
*
* 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)
*
* Get the data range by running a (silent) OLS. Use SMPL=%VALID(U) after
* for the SMPL after this. (That will take care both of series knocked
* out by the SMPL and those knocked out by missing data.)
*
linreg(noprint,smpl=smpl,define=baseeq) datavec(1) start end u
# rhs detvars
compute startl=%regstart(),endl=%regend(),regnobs=%nobs
compute baseregs=%nreg
*
* Figure out the number of covered individuals and the number of data
* points for each.
*
compute starti=%indiv(startl),endi=%indiv(endl)
panel(icount=1.0,smpl=%valid(u),compress) u startl endl icount
set icount  starti endi = %if(%valid(icount),icount,0.0)
sstats(max) starti endi icount>>pt
sstats      starti endi (icount>0)>>nt
compute tmax=fix(pt),tend=%panelobs()
compute nsecs=fix(nt)
*
compute unbal=(regnobs<>tmax*nsecs)
*
* 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(%nreg)
*
* 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=%rlconcat(reglist,detvars)
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(baseregs,1)
compute pdlbsum2 = %zeros(baseregs,1)
compute pdlbsum3 = %zeros(baseregs,1)
compute pdltsum  = %zeros(baseregs,1)
compute pdlvsum1 = %zeros(baseregs,baseregs)
compute pdlvsum2 = %zeros(baseregs,baseregs)
compute pdlhsum2 = %zeros(baseregs,baseregs)
compute pdlhsum3 = %zeros(baseregs,baseregs)
*
dim %%ibetas(baseregs,nsecs) %%istderrs(baseregs,nsecs) %%itstats(baseregs,nsecs)
dim %%ixx(nsecs)
*
* Run loop over individuals, but skip any which have no data. All the
* statistics are dimensioned <<NSECS>>, so we fill them based upon the
* running value of <<ISEC>> rather than the absolute individual number.
*
infobox(action=define,progress,lower=starti,upper=endi) "Computing Panel Statistics"
compute isec=0
do J=starti,endi
   if icount(j)==0
      next
   compute isec=isec+1
   *
   * Copy the data from ytilde into dvec for this individual. Difference
   * and de-mean the differences.
   *
   do k=1,m+1
      set dvec(k)  j//1 j//tend = %if(%valid(u),ytilde(k){0},%na)
      diff dvec(k) j//2 j//tend rvec(k)
      diff(center) rvec(k) j//2 j//tend
   end do k
   *
   * Copy the subvector of data and differenced data for the RHS variables
   *
   do K=1,M
      set xvec(K)  j//1 j//tend = dvec(K+1){0}
      set dxvec(K) j//2 j//tend = rvec(K+1){0}
   end do K
   *
   * Do DOLS augmented regression
   *
   linreg(noprint) dvec(1) j//(dlags+2) j//(tend-dlags) dolsres
   # reglist
   *
   * Compute the long-run variance of the residuals.
   *
   mcov(lwindow=bartlett,lags=lags,noprint) %regstart() %regend()
   # dolsres
   compute idlrvar = %cmom(1,1)/%nobs
   *
   * 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,baseregs)
   compute vdols  = idlrvar*%xsubmat(%xx,1,baseregs,1,baseregs)
   compute hdols  = inv(vdols)
   compute ddols  = inv(%sqrt(%diag(%xdiag(vdols))))
   *
   * Compute t-stat vs the hypothesized vector
   *
   ewise tdols(i)  = (bdols(i)-%if(i<=m,bvec(i),0.0))*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
   *
   * Save the individual coefficients, standard errors and t-statistics
   * in the global variables.
   *
   compute %psubmat(%%ibetas,1,isec,bdols)
   compute %psubmat(%%itstats,1,isec,tdols)
   compute %psubmat(%%istderrs,1,isec,%sqrt(%xdiag(vdols)))
   compute %%ixx(isec)=vdols
   *
   * Update the infobox
   *
   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 global variables for overall estimator
*
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
compute %stderrs=%sqrt(%xdiag(%xx))
compute %ngroup =nsecs
*
compute %stderrs=%sqrt(%xdiag(%xx))
*
* Define the full-sample residuals
*
set %resids startl endl = %if(%valid(u),%eqnrvalue(baseeq,t,%beta),%na)

    *** DISPLAY INDIVIDUAL AND PANEL RESULTS ***

if %defined(title)
   compute ltitle=title
else
   compute ltitle="Mean Group Panel DOLS Estimation"

report(use=preport,action=define,title=ltitle)
report(use=preport,atrow=1,atcol=1,span) ltitle
report(use=preport,atrow=2,atcol=2) "LHS Variable" labelvec(1)
report(use=preport,atrow=3,atcol=2) "Individuals" NSecs
if unbal==0
   report(use=preport,atrow=4,atcol=2) "Time Periods" tmax
else
   report(use=preport,atrow=4,atcol=2) "Time Periods" "Unbalanced"
report(use=preport,atrow=5,atcol=2) "Dynamic Lags" dlags
if tdum == 0 ; { ; report(use=preport,atrow=6,atcol=2,span) "common time dummies NOT included" ; }
if tdum == 1 ; { ; report(use=preport,atrow=6,atcol=2,span) "common time dummies included" ; }
report(use=preport,atrow=8,atcol=2) "RHS Variable" "H0 Coefficient"
do k=1,m
   report(use=preport,atrow=8+k,atcol=2) labelvec(k+1) bvec(k)
end do k
compute row=8+m+2
report(use=preport,atrow=row,atcol=1) "Member" "Variable" "Coefficient" "t-Statistic"
if print==3 {
   compute isec=0
   do i=starti,endi
      if icount(i)==0
         next
      compute isec=isec+1
      do k=1,m
         if k==1
            report(use=preport,row=new,atcol=1) "No."+i labelvec(k+1) %%ibetas(k,isec) %%itstats(k,isec)
         else
            report(use=preport,row=new,atcol=2)         labelvec(k+1) %%ibetas(k,isec) %%itstats(k,isec)
      end do k
   end do i
   report(use=preport,row=new,atcol=1,span) "*************"
}
do k=1,m
   if k==1
      report(use=preport,row=new,atcol=1) "Group" labelvec(k+1) %beta(k) pdltsum(k)
   else
      report(use=preport,row=new,atcol=2)         labelvec(k+1) %beta(k) pdltsum(k)
end do k
if print>1
   report(use=preport,action=show)

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


