*
* @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
*
*   DET=NONE/[CONSTANT]/TREND
*   TDUM/[NOTDUM]
*      Use TDUM to subtract out common time effects
*
*   SMPL=Dummy series with 0's in entries to skip [include all]
*
*   UNWEIGHTED/[NOUNWEIGHTED]
*      Use UNWEIGHTED for unweighted statistics
*   LAGS=(maximum number of) augmenting lags in the ADF regressions [Schwert's]
*   CRIT=[FIXED]/GTOS/AIC/BIC/HQ
*     Criterion to use to select the number of lags. This is done on an
*     individual by individual basis. CRIT=FIXED uses the input LAGS on
*     each. CRIT=GTOS starts with LAGS and drops lags until it hits one
*     which has the marginal t with significance level set by the SLSTAY
*     option. CRIT=AIC/BIC/HQ select the lag length using AIC, BIC (or SBC)
*     or HQ.
*   SLSTAY=significance level to keep lag in model with METHOD=GTOS [.10]
*
*   [PRINT]/NOPRINT
*   TITLE=title of report ["Panel Cointegration Tests"]
*
* Maintained for compatibility. Use LAGS and DET options instead.
*
*   MLAG  = starting value for maximum number of lags in the ADF based tests [3]
*   TREND/[NOTREND]
*      Use TREND to include heterogeneous member specific trends
*
* 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
*
* 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.
*   04/2012 Extensive revisions by Tom Doan, Estima
*           Switch to use of REPORT
*           Add standard DET option (deprecating the TREND option)
*           Add CRIT and SLSTAY option
*           Add TITLE and PRINT options
*           Add SMPL option
*           Add check for more than 8 endogenous variables (correction
*           tables don't go higher than that).
*   05/2017 Expand description of N(0,1) test statistics
*
************************************************************************
procedure PanCoint start end
type integer start end

option integer  lags
option choice   crit        1  fixed gtos aic bic hq
option real     slstay     .10
option series   smpl
option choice   det         2  none constant trend
option switch   unweighted  0
option switch   tdum        0
option switch   print       1
option string   title
*
option integer  mlag
option switch   trend		 0
option integer  block

*** 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 m maxlag klag
local series        time ehat ehatlag icount u
local series        dehat nhat destar estar uresid
local integer       llags startl endl nsecs
local real          su st lambda mtratio L11var sum
local real          marginal_t logl ic bestic
local vect[integer] datavec rhs tvec
local vect[integer] lagsused
local integer       unbal
local report        preport
local vect[int]     detvars
local integer       tablenum
local rect          sweepxx
local string        ltitle descript
local integer       p slot isec pmax tend tmax starti endi regnobs
local real          pt nt
*
* Lag criterion descriptors
*
fixed vect[string] lagcriterion(5)
Fixed
General-to-specific
AIC
BIC/SBC
Hannan-Quinn
*
* Lookup tables for adjustments without (element 1) and with (2) 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
*
* Check for valid settings
*
if %defined(block) {
   disp "###@PANCOINT(BLOCK=nnn). No longer supported. Use PANEL calendar scheme"
   return
}
if det==1 {
   disp "###@PANCOINT doesn't support DET=NONE"
   return
}
set time = t
*
* Figure out which deterministics are requested.
*
if trend
   compute detvars=||constant,time||
else
if det==3
   compute detvars=||constant,time||
else
   compute detvars=||constant||

if %size(detvars)==1
   compute tablenum=1,descript="Constant"
else
   compute tablenum=2,descript="Constant+Trend"
*
* Pull in list of series. Set m to the number of variables - 1
*
enter(varying) datavec
compute m=%rows(datavec)-1
if m>7 {
   display "###@PANCOINT cannot handle more than 8 endogenous variables"
   return
}
*
* 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. These are Im-Pesaran-Shin (where
* appropriate) for m=0 (unit root tests) and Pedroni's for m>0
*
if m>=1 {
   compute muvec=||tab(tablenum)(m,1),tab(tablenum)(m,3),tab(tablenum)(m,5),tab(tablenum)(m,5),tab(tablenum)(m,7),tab(tablenum)(m,9),tab(tablenum)(m,9)||
   compute vivec=||tab(tablenum)(m,2),tab(tablenum)(m,4),tab(tablenum)(m,6),tab(tablenum)(m,6),tab(tablenum)(m,8),tab(tablenum)(m,10),tab(tablenum)(m,10)||
}
else
if tablenum==1 {
   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||
}
*
* 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) datavec(1) start end u
# rhs
compute startl=%regstart(),endl=%regend(),regnobs=%nobs
*
* 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)
*
dim dvec(m+1) rvec(m+1)
dim statvec(7)
dim tsizes(nsecs) lagsused(nsecs)
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)
*******************************************************************************
*
* Setup data variables
*
if tdum == 1 {
   *
   * Extract fixed time effects if requested. Write transformed series
   * to TDUMSERIES so we don't overwrite original data.
   *
   dim tdumseries(m+1)
   do K=1,m+1
      panel(entry=1.0,time=-1.0,smpl=%valid(u)) datavec(K) startl endl tdumseries(K)
   end do K
   *
   * Reset DATAVEC to point to the new series:
   *
   ewise datavec(i) = tdumseries(i)
}
*
* Figure out the requested lag limit for ADF regressions
*
if %defined(mlag)
   compute pmax=mlag
else
if %defined(lags)
   compute pmax=lags
else
   compute pmax=fix(%round(12*(.01*tmax)^.25,0))
*
* 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.
*
compute isec=0
do J=starti,endi
   if icount(j)==0
      next
   compute isec=isec+1
   compute tsizes(isec)=icount(isec)
   *
   * Copy datavec's to the dvec series and difference
   *
   do k=1,m+1
      set dvec(k) j//1 j//tend = %if(%valid(u),datavec(k){0},%na)
      set rvec(k) j//2 j//tend = dvec(k){0}-dvec(k){1}
   end do k
   ****************************************************************************
   *
   * Do individual cointegrating regressions
   *
   *
   * ehat are residuals from cointegrating regression with deterministics
   *
   if m >= 1 {
      linreg(noprint) dvec(1) j//1 j//tend ehat
      # dvec(2) to dvec(m+1) detvars
   }
   else {
      linreg(noprint) dvec(1) j//1 j//tend ehat
      # detvars
   }
   diff ehat j//2 j//tend dehat
   *
   * nhat are residuals from differences on differences
   *
   if m>=1 {
      linreg(noprint) rvec(1) j//2 j//tend nhat
      # rvec(2) to rvec(m+1)
   }
   else {
      linreg(noprint) dehat j//2 j//tend nhat
      # ehat
   }
   *
   * Set lagsused to 0, which will be the value if we don't find anything
   * significant at 1 or higher.
   *
   if crit==1 {
      compute lagsused(isec)=pmax
      goto chosen
   }
   compute lagsused(isec)=0
   cmom j//(pmax+2) j//tend
   # dehat ehat{1} dehat{1 to pmax}
   *
   * Start with the full regression
   *
   compute sweepxx=%sweeplist(%cmom,%seq(2,%ncmom))
   compute slot=%ncmom
   do p=pmax,(crit==2),-1
      compute %ndf=%nobs-(slot-1)
      if crit==2 {
         *
         * Compute the (absolute) t-statistic from the information in
         * <<sweepxx>>. The pivot is at <<slot>>, which is the number of
         * regressors + 1 for the dependent variable.
         *
         compute marginal_t=abs(sweepxx(slot,1))/$
           sqrt(sweepxx(slot,slot)*sweepxx(1,1)/%ndf)
         *
         * If that's significant at the requested level, break the loop.
         *
         if marginal_t>%invttest(slstay,%ndf) {
            compute lagsused(isec)=p
            break
         }
      }
   		else
      if crit==3.or.crit==4.or.crit==5 {
         compute logl=-.5*%nobs*log(sweepxx(1,1)/%nobs)
         if crit==3
            compute ic=2.0*logl-2*p
         else
         if crit==4
            compute ic=2.0*logl-log(%nobs)*p
         else
            compute ic=2.0*logl-2*log(log(%nobs))*p
         if p==pmax.or.ic>bestic
            compute bestic=ic,lagsused(isec)=p
      }
      *
      * Sweep on the last lag to take it out of the regression.
      *
      compute sweepxx=%sweep(sweepxx,slot)
      compute slot=slot-1
   end do p
   :chosen
   *
   * Compute lag width for Bartlett kernel
   *
   compute klag = fix(%round((4.0*(tsizes(isec)/100.0)^(2.0/9.0)),0))
   *
   * Compute individual member sample stats
   *
   * Sweep out the augmenting lags and save the cross products of lagged
   * ehat with itself and with dehat.
   *
   *
   compute p=lagsused(isec)
   if (p == 0) {
      cmoment(noprint) j//(2+p) j//tend
      # dehat ehat{1}
      compute sweepxx=%cmom
   }
   else {
      cmoment(noprint) j//(2+p) j//tend
      # dehat ehat{1} dehat{1 to p}
      compute sweepxx=%sweeplist(%cmom,%seq(3,%ncmom))
   }
   compute nadf(isec) = sweepxx(2,1)
   compute dadf(isec) = sweepxx(2,2)
   *
   * Sweep out the lagged residual and save the sigma-hat squared.
   *
   compute sweepxx=%sweep(sweepxx,2)
   compute sadf(isec) = sweepxx(1,1)/(%nobs-%ncmom+1)
   *
   * Compute regression w/o augmenting lags and get the long-run
   * variance of its residuals.
   *
   linreg(noprint) dehat j//2 j//tend uresid
   # ehat{1}
   cmoment(noprint,lastreg) j//2 j//tend
   compute d(isec) = %cmom(1,1)
   compute n(isec) = %cmom(2,1)
   mcov(lwindow=newey,lags=klag) j//2 j//tend
   # uresid
   compute su = sqrt(%rss/(tsizes(isec)-1.0))
   compute st = sqrt(%cmom(1,1)/(tsizes(isec)-1.0))
   compute lambda = 0.5*(st^2-su^2)
   compute S(isec) = st^2
   *
   * Compute the long run variance of nhat
   *
   mcov(lwindow=newey,lags=klag,noprint) j//2 j//tend
   # nhat
   compute L11var = (1.0/(tsizes(isec)-1.0))*%cmom(1,1)
   *
   * Construct numerator and denominator terms
   *
   if unweighted == 1 {
      compute Lbar(isec) = L11var
      compute L11var = 1.0
   }

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

   compute GR(isec)    = N2(isec)/D2(isec)
   compute GP(isec)    = N2(isec)/sqrt(D2(isec)*S(isec))
   compute GA(isec)    = Nadf(isec)/sqrt(Dadf(isec)*Sadf(isec))
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 %defined(title)
   compute ltitle=title
else
   compute ltitle="Panel Cointegration Tests"
*
report(use=preport,action=define,title=ltitle)
report(use=preport,atrow=1,atcol=1,span) ltitle

if m>=1 {
   report(use=preport,row=new,atcol=1) "panel v-stat"    statvec(1)
   report(use=preport,row=new,atcol=1) "panel rho-stat"  statvec(2)
   report(use=preport,row=new,atcol=1) "panel pp-stat"   statvec(3)
   report(use=preport,row=new,atcol=1) "panel adf-stat"  statvec(4)
   report(use=preport,row=new)
   report(use=preport,row=new,atcol=1) "group rho-stat"  statvec(5)
   report(use=preport,row=new,atcol=1) "group pp-stat"   statvec(6)
   report(use=preport,row=new,atcol=1) "group adf-stat"  statvec(7)
   report(use=preport,row=new)
}

if m==0 {
   report(use=preport,row=new,atcol=1,span) "raw panel unit root test results"
   report(use=preport,row=new)
   report(use=preport,row=new,atcol=1) "Levin-Lin rho-stat"   statvec(2)
   report(use=preport,row=new,atcol=1) "Levin-Lin t-rho-stat" statvec(3)
   report(use=preport,row=new,atcol=1) "Levin-Lin ADF-stat"   statvec(4)
   report(use=preport,row=new)
   report(use=preport,row=new,atcol=1) "IPS ADF-stat"         statvec(7)
   report(use=preport,row=new,atcol=1,span) "(using large sample adjustment values)"
   report(use=preport,row=new)
}

report(use=preport,row=new,atcol=1) "Individuals" NSecs
if unbal
   report(use=preport,row=new,atcol=1) "Unbalanced"
else
   report(use=preport,row=new,atcol=1) "Time Periods" tmax
report(use=preport,row=new,atcol=1) "Regressors"   M

report(use=preport,row=new)
report(use=preport,row=new,atcol=1,span) "All reported values are distributed N(0,1) under the null"
report(use=preport,row=new,atcol=1,span) "of unit root or no cointegration with rejection in the left tail."
report(use=preport,row=new,atcol=1,span) "Critical values are 10%=-1.28,5%=-1.64,1%=-2.33"
report(use=preport,row=new)
if unweighted == 0
{
   report(use=preport,row=new,atcol=1,span) "Panel stats are weighted by long run variances"
   report(use=preport,row=new,atcol=1,span) "(see OBES reference for details)"
}
if unweighted == 1
{
   report(use=preport,row=new,atcol=1,span) "Panel stats are unweighted by long run variances"
   report(use=preport,row=new,atcol=1,span) "(see ET reference for details)"
}
report(use=preport,action=format,picture="*.##",atcol=2,align=decimal)
if print
   report(use=preport,action=show)
end
