*
* @pancoint( options ) start end
* # list of variables (list dependent variable first)
*
* Procedure for Cointegration Tests in Heterogeneous Panels with
* Multiple Regressors ("Pedroni tests").
*
* 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
*
* Supplementary card:
*   If you list just one variable, it will do a raw panel unit root test
*
* Options:
*   mlag  = starting value for maximum number of lags in the ADF based tests [3]
*   block = time dimension [from CALENDAR(PANELOBS)] (see note below for unbalanced panels)
*   trend/[notrend] Use TREND to include heterogeneous member specific trends
*   tdum/[notdum]   Use TDUM to subtract out common time effect
*   unweighted/[nounweighted]  Use UNWEIGHTED for unweighted statistics
*
* 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:
*
* For distribution theory, Monte Carlos and m = 1 critical values, and discussion
* of unweighted versions of pooled stats:
*
* Pedroni, Peter (2004) "Panel Cointegration; Asymptotic and Finite
* Sample Properties of Pooled Time Series Tests with an Application to
* the PPP Hypothesis," Econometric Theory, 20, 597-625.
*
* Working paper version: Pedroni, Peter (1995) Indiana University
* Working Papers in Economics, No. 95-013. Revised 4/97,8/01,8/02.
*
* For general discussion, ADF based tests, and m>1 critical values:
*
* Pedroni, Peter (1999) "Critical Values for Cointegration Tests in
* Heterogeneous Panels with Multiple Regressors," Oxford Bulletin of
* Economics and Statistics, 61, 653-70.
*
* Revisions:
*   Written by Peter Pedroni, Williams College: ppedroni@williams.edu
*   Extended by author from original version to allow for common time effects
*   in unbalanced panels, add an option for weighted vs unweighted stats,
*   and automatic selection of bandwidth for kernels
*
*   06/2004 Converted to procedure by Tom Doan, Estima
*   07/2004 Updated by Tom Maycock, Estima to include the new revisions
*           described above.
*   12/2004 Added code to deal with missing values in the first period(s) of
*           an individual.
*   06/2005 Minor tweak to MOVE command to be more technically correct. Change
*           shouldn't affect results.
*   07/2007 Add TDUM and UNWEIGHTED options.
*   07/2009 Uncommented the "save calendar" instructions, so that retaining
*            the original CALENDAR is the default behavior.
*   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".
*             Plus other changes needed accordingly.
*
************************************************************************
procedure PanCoint startp endp
type integer startp endp

option integer mlag 3
option switch  trend
option integer block

*** 07/2007 Revision: New options:
option switch unweighted 0
option switch tdum 0


*** MAIN SOURCE CODE ***

local vector muvec vivec tvecreal
local vec[series] dvec rvec tdumseries
local vec[real] N D N2 D2 N3 D3 Nadf Dadf Sadf S N3adf D3adf $
  S3adf S3  V3 DD  DD3  NN3 GR GP GA Lbar statvec tsizes
local integer pull i j k maxlag klag
local series time ehat ehatlag icount u
local series dehat nhat destar estar uresid
local integer llags start end nsecs
local real su st lambda mtratio L11var sum
local vect[integer] datavec rhs tvec
local integer m
local integer unbal
local integer tperiods offset tperiodsmax entry i
local vector[real] offsettemp
local integer     firstInd LastInd
*
* Lookup tables with and without trend
*

fixed vect[rect] tab(2)(7,10)

   8.62,  60.75,   -6.02,  31.27,  -1.73,
   0.93,  -9.05,   35.98,  -2.03,   0.66
 11.754, 104.546, -9.495,  57.610, -2.177,
  0.964, -12.938,  51.49,  -2.453,  0.618
 15.197, 151.094, -13.256, 81.772, -2.576,
  0.923, -16.888,  67.123, -2.827,  0.585
 18.910, 190.661, -17.163, 99.331, -2.930,
  0.843, -20.841,  81.835, -3.157,  0.560
 22.715, 231.864, -21.013,119.546, -3.241,
  0.800, -24.775,  98.278, -3.452,  0.553
 26.603, 270.451, -24.944,134.341, -3.531,
  0.750, -28.720, 113.131, -3.726,  0.542
 30.457, 293.431, -28.795,144.615, -3.795,
  0.685, -32.538, 126.059, -3.976,  0.525


  17.86, 101.68,  -10.54,   39.52,  -2.29,
   0.66, -13.65,   50.91,   -2.53,   0.56
 21.162, 160.249, -14.011,  64.219, -2.648,
  0.690, -17.359,  66.387,  -2.872,  0.555
 24.556, 198.167, -17.600,  83.815, -2.967,
  0.686, -21.116,  81.832,  -3.179,  0.548
 28.046, 239.425, -21.287, 103.905, -3.262,
  0.688, -24.930,  97.362,  -3.464,  0.543
 31.738, 276.997, -25.130, 124.613, -3.545,
  0.686, -28.849, 113.145,  -3.737,  0.538
 35.537, 310.982, -28.981, 138.227, -3.806,
  0.654, -32.716, 127.989,  -3.986,  0.530
 39.231, 348.217, -32.756, 154.378, -4.047,
  0.638, -36.494, 140.756,  -4.217,  0.518
*
* Pull in list of series. Set m to the number of variables - 1
*
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)
*
* Generate appropriate lookup values
*
if m>=1 {
   compute pull=trend+1
   compute muvec=||tab(pull)(m,1),tab(pull)(m,3),tab(pull)(m,5),tab(pull)(m,5),tab(pull)(m,7),tab(pull)(m,9),tab(pull)(m,9)||
   compute vivec=||tab(pull)(m,2),tab(pull)(m,4),tab(pull)(m,6),tab(pull)(m,6),tab(pull)(m,8),tab(pull)(m,10),tab(pull)(m,10)||
}
else
if trend==0 {
   compute muvec = ||0.00, -3.00, -1.53, -1.53, 0.00, 0.00, -1.54||
   compute vivec = ||0.00, 10.20, 1.25, 1.25, 0.00, 0.00, 0.71||
}
else {
   compute muvec = ||0.00, -7.50, -1.94, -1.94, 0.00, 0.00, -2.18||
   compute vivec = ||0.00, 5.76, 1.62, 1.62, 0.00, 0.00, 0.56||
}

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

compute tperiodsmax = tperiods

*
* 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

*** Fix to allow starting individual>1
compute firstInd = (start/tperiods)+1
compute lastInd = firstInd+(nsecs-1)


display ' '
display '      Currently computing panel statistics. Please wait. '
display ' '



*** 07/2007 Revision: New treatment for observation counts
***                 and determining if unbalanced:

*
* Run an OLS mainly to figure out what data points are available.
linreg(noprint) datavec(1) start end u
# rhs

* NOTE: If using an old version, you can comment out the lines
* used to save and restore the original CALENDAR setting:
*
compute savecal = %calendar()
*
dim tvec(nsecs)
calendar(panelobs=tperiods)

*** Fix to allow starting individual>1
**** panel(icount=1,compress) u / icount
panel(icount=1,compress) u start end icount

*** Fix to allow starting individual>1
*** ewise tvec(i) = fix(icount(i))
ewise tvec(i) = fix(icount(i+firstind-1))

*
calendar(recall=savecal)
*
*
* Create in tsizes the observation counts for the individuals
*
dim tvecreal(nsecs) tsizes(nsecs)
ewise tvecreal(i) = tvec(i)
compute unbal = fix(%minvalue(tvecreal)<%maxvalue(tvecreal))
if unbal
   ewise tsizes(i)=tvec(i)
else
   ewise tsizes(i)=TPeriods


dim dvec(m+1) rvec(m+1)
dim N(Nsecs) D(Nsecs) N2(Nsecs) D2(Nsecs) N3(Nsecs) D3(Nsecs) $
  Nadf(Nsecs) Dadf(Nsecs) Sadf(Nsecs) S(Nsecs) N3adf(Nsecs) D3adf(Nsecs) $
  S3adf(Nsecs) S3(Nsecs) V3(Nsecs) DD(Nsecs) DD3(Nsecs) NN3(Nsecs) $
  GR(Nsecs) GP(Nsecs) GA(Nsecs) Lbar(Nsecs) statvec(7)


               ***  SETUP DATA VARIABLES ***


if tdum == 1
{
 dim tdumseries(m+1)
 compute savecal = %calendar()
 calender(panelobs=tperiods)
 * Write transformed series to TDUMSERIES
 * so we don't overwrite original data:
 do K=1,m+1
  panel(entry=1,time=-1) datavec(K) / tdumseries(K)
 end do K
 calander(recall=savecal)
 * Reset DATAVEC to point to the new series:
 ewise datavec(i) = tdumseries(i)
}

* Dimension array used to hold temporary data when
* determining offset:
dim offsettemp(m+1)
do J=firstInd,lastInd

   *** Fix to allow starting individual>1
   *** While J runs from firstInd,lastInd,
   *** IndCounter will run from J to Nsecs,
   *** for arrays that need to start from 1.
   *** Replace previous references to (J) with (IndCounter)

   compute IndCounter = J-firstInd+1

   ***  if unbal == 1 ; {; compute Tperiods = fix(Tsizes(J)); }
   if unbal == 1 ; {; compute Tperiods = fix(Tsizes(IndCounter)); }

   * Need to allow for possibility of NA's at the beginning of
   * individual. Compute offset from entry indiv//1 to first entry
   * in that individual for which all series contain data:

   do entry=((j-1)*tperiodsmax+1),j*tperiodsmax
      ewise offsettemp(i) = ([series]datavec(i))(entry)
      * If using Version 6 or later, you can replace the next
      * 5 lines with the single command:
      *if %valid(offsettemp)

      compute sum = 0
      do i=1,%rows(offsettemp)
        compute sum = sum + offsettemp(i)
      end do i
      if %valid(sum)
        {
        compute offset = entry-((j-1)*tperiodsmax)
        break
        }
   end do

   do K=1,M+1
*      move datavec(K) ((j-1)*Tperiods+offset) (j*Tperiods+offset-1) dvec(K) 1
* June 05 change to MOVE command--avoid copying over extra data:
*      move datavec(K) ((j-1)*Tperiodsmax+offset) (j*Tperiodsmax+offset-1) dvec(K) 1
      move datavec(K) ((j-1)*Tperiodsmax+offset) (j*Tperiodsmax) dvec(K) 1

      diff dvec(K) 2 Tperiods rvec(K)
   end do K
   set time = T

         *** DO INDIVIDUAL COINTEGRATING REGRESSIONS ***

  if m >= 1
  {
   linreg(noprint) dvec(1) 1 Tperiods ehat
   # dvec(2) to dvec(m+1) constant
  }

  if trend == 1 .and. m>=1
  {
   linreg(noprint) dvec(1) 1 Tperiods ehat
   # dvec(2) to dvec(m+1) constant time
   }

  if m == 0
  {
   linreg(noprint) dvec(1) 1 Tperiods ehat
   # constant
  }

  if m == 0 .and. trend == 1
  {
   linreg(noprint) dvec(1) 1 Tperiods ehat
   # constant time
  }

  if m>=1
  {
   linreg(noprint) rvec(1) 2 Tperiods nhat
   # rvec(2) to rvec(m+1)
  }

  diff ehat / dehat

  if m==0
  {
   linreg(noprint) dehat 2 Tperiods nhat
   # ehat
  }

          *** COMPUTE ADF LAG TRUNCATIONS ***

   do llags=mlag,1,-1
      linreg(noprint) dehat 2 Tperiods
      # ehat{1} dehat{1 to llags}
      compute mtratio = %tstats(llags+1)
      if abs(mtratio) >= 1.64
         { ; compute maxlag=llags ; break ; }
   end do llags

   if llags == 1 .and. abs(mtratio) < 1.64
      compute maxlag = 0

***   compute klag = fix(%round((4.0*(Tvec(J)/100.0)**(2.0/9.0)),0))
   compute klag = fix(%round((4.0*(Tvec(IndCounter)/100.0)**(2.0/9.0)),0))



            *** COMPUTE INDIVIDUAL MEMBER SAMPLE STATS ***

   if maxlag == 0 {
      linreg(noprint) dehat 2 Tperiods
      # ehat{1}
      compute sadf(IndCounter) = %seesq
      cmoment(noprint) 2+maxlag Tperiods
      # dehat ehat{1}
      compute nadf(IndCounter) = %cmom(2,1)
      compute dadf(IndCounter) = %cmom(2,2)
   }

   if maxlag >= 1 {
      linreg(noprint) dehat 2 Tperiods
      # ehat{1} dehat{1 to maxlag}
      compute sadf(IndCounter) = %seesq
      linreg(noprint) dehat 2 Tperiods destar
      # dehat{1 to maxlag}
      set ehatlag 1 Tperiods = ehat{1}
      linreg(noprint) ehatlag 2 Tperiods estar
      # dehat{1 to maxlag}
      cmoment(noprint) 2+maxlag Tperiods
      # destar estar
      compute nadf(IndCounter) = %cmom(2,1)
      compute dadf(IndCounter) = %cmom(2,2)
   }

   linreg(noprint) dehat 2 Tperiods uresid
   # ehat{1}
   compute d(IndCounter) = 1.0/%xx(1,1)
   cmoment(noprint,lastreg) 2 Tperiods
   compute n(IndCounter) = %cmom(2,1)
   mcov(damp=1.0,lags=klag) 2 Tperiods
   # uresid
   compute su = sqrt(%seesq*%ndf/(Tvec(IndCounter)-1.0))
   compute st = sqrt(%cmom(1,1)/(Tvec(IndCounter)-1.0))
   compute lambda = 0.5*(st**2-su**2)
   compute S(IndCounter) = st**2

   mcov(damp=1.0,lags=klag,noprint) 2 Tperiods
   # nhat
   compute L11var = (1.0/(tvec(IndCounter)-1.0))*%CMOM(1,1)


          *** CONSTRUCT NUMERATOR AND DENOMINATOR TERMS ***

   if unweighted == 1
   {
     compute Lbar(IndCounter) = L11var
     compute L11var = 1.0
    }

   compute N2(IndCounter) = N(IndCounter) - lambda*Tvec(IndCounter)
   compute D2(IndCounter) = D(IndCounter)
   compute N3(IndCounter) = N2(IndCounter)/L11var
   compute D3(IndCounter) = D2(IndCounter)/L11var
   compute N3adf(IndCounter) = Nadf(IndCounter)/L11var
   compute D3adf(IndCounter) = Dadf(IndCounter)/L11var
   compute S3adf(IndCounter) = Sadf(IndCounter)/L11var
   compute S3(IndCounter) = S(IndCounter)/L11var
   compute DD(IndCounter) = D(IndCounter)*S(IndCounter)
   compute DD3(IndCounter) = D3(IndCounter)*S3(IndCounter)

   compute GR(IndCounter) = N2(IndCounter)/D2(IndCounter)
   compute GP(IndCounter) = N2(IndCounter)/sqrt(D2(IndCounter)*S(IndCounter))
   compute GA(IndCounter) = Nadf(IndCounter)/sqrt(Dadf(IndCounter)*Sadf(IndCounter))
end do J


        *** CONSTRUCT PANEL COINTEGRATION STATISTICS ***

ewise V3(I) = D3(I)*(1.0/(Tsizes(I)**2))
ewise NN3(I) = N3(I)*(Tsizes(I)/1.0)
ewise GR(I) = GR(I)*(Tsizes(I)/1.0)


compute statvec(1) = sqrt(Nsecs**3)/%sum(V3)
if unweighted == 1
{
 compute statvec(1) = (sqrt(Nsecs**3)*%sum(Lbar))/(%sum(V3)*Nsecs)
}
compute statvec(2) = (sqrt(Nsecs)*%sum(NN3))/%sum(D3)
compute statvec(3) = %sum(N3)/(sqrt(%sum(D3)*%sum(S3)/Nsecs))
compute statvec(4) = %sum(N3adf)/(sqrt(%sum(D3adf)*(%sum(S3adf)/Nsecs)))
compute statvec(5) = %sum(GR)/sqrt(Nsecs)
compute statvec(6) = %sum(GP)/sqrt(Nsecs)
compute statvec(7) = %sum(GA)/sqrt(Nsecs)

do K=1,7
   compute statvec(K) = (statvec(K) - muvec(K)*sqrt(Nsecs))/sqrt(vivec(K))
end do K


                *** DISPLAY RESULTS ***

if m>=1 {
   display '                           RESULTS:'
   display '         ******************************************** '
   display '               panel v-stat     = ' statvec(1)
   display '               panel rho-stat   = ' statvec(2)
   display '               panel pp-stat    = ' statvec(3)
   display '               panel adf-stat   = ' statvec(4)
   display ' '
   display '               group rho-stat   = ' statvec(5)
   display '               group pp-stat    = ' statvec(6)
   display '               group adf-stat   = ' statvec(7)
   display ' '

}

if m==0 {
   display '                           RESULTS:'
   display '         ******************************************** '
   display '             -raw panel unit root test results-       '
   display ' '
   display '             Levin-Lin rho-stat   = ' statvec(2)
   display '             Levin-Lin t-rho-stat = ' statvec(3)
   display '             Levin-Lin ADF-stat   = ' statvec(4)
   display ' '
   display '                 IPS ADF-stat   = ' statvec(7)
   display '             (using large sample adjustment values)    '
   display ' '
}

if unbal == 0
   display '        Nsecs =' Nsecs ', Tperiods =' Tperiods ', no. regressors =' M
if unbal == 1
   display '        Nsecs =' Nsecs ', Tperiods = (unbalanced)' ', no. regressors =' M


display ' '
display '            All reported values are distributed N(0,1)'
display '            under null of unit root or no cointegration'
display ' '
if unweighted == 0
{
display '            Panel stats are weighted by long run variances'
display '                (see OBES reference for details)'
}
if unweighted == 1
{
display '            Panel stats are unweighted by long run variances'
display '                 (see ET reference for details)'
}

display '         ******************************************** '

end










