*
* @regPCSE
* computes a Panel-Corrected version of an OLS regression on a panel
* data set. This should be used after the LINREG which does the OLS
* estimates.
*
* Reference: Beck and Katz (1995), "What to Do (and Not to Do) with
* Time-Series-Cross Section Data in Comparative Politics," American
* Political Science Review, vol 89, no 3, pp 634-647.
*
* Options:
*  METHOD=[PCSE]/IHAC/PHAC
*  LAGS=number of lags in IHAC and PHAC. A Bartlett (Newey-West) window
*   is used.
*  TITLE=title for regression output ["OLS with (description) Covariance Matrix"]
*
* Revision Schedule:
*  02/2008 Written by Tom Doan, Estima
*  05/2012 Revised to add the METHOD option for IHAC and PHAC.
*
procedure regPCSE
*
option choice   method  1 pcse ihac phac
option string   title
option integer  lags    0
*
local rect 			 xt xtlag xlx
local symm 			 xvx pcsexx sigma
local integer   i j l
local integer   startl endl np tp rnobs
local vector    u id xv wt
local real      tx
local string    ltitle
local vect[vect[int]] identries
*
compute rnobs =%nobs
compute startl=%regstart(),endl=%regend()
panel(effects=indiv,id=id,identries=identries,smpl=%valid(%resids))
*
compute np=%size(id)
sstats(max) 1 np %size(identries(t))>>tx
compute tp=fix(tx)
if rnobs<>np*tp {
   disp "###@REGPCSE requires balanced sample"
   return
}
if %defined(title)
   compute ltitle=title
dim u(np)
dim xt(np,%nreg) xtlag(np,%nreg)

if (method==1) {
   *
   * Do Beck & Katz standard errors
   *
   compute sigma=%zeros(np,np)
   do t=1,tp
      ewise u(i)=%resids(identries(i)(t))
      compute sigma=sigma+%outerxx(u)
   end do t
   *
   compute sigma=sigma/tp
   *
   compute xvx=%zeros(%nreg,%nreg)
   do t=1,tp
      ewise xt(i,j)=xv=%eqnxvector(0,identries(i)(t)),xv(j)
      compute xvx=xvx+%mqform(sigma,xt)
   end do t
   if .not.%defined(title)
      compute ltitle="OLS with PCSE Covariance Matrix"
}
*
else {
   dim wt(lags)
   ewise wt(i)=1.0-i/(lags+1.0)
   compute xvx=%zeros(%nreg,%nreg)
   if .not.%defined(title) {
      if (method==2)
         compute ltitle="OLS with Individual HAC Covariance Matrix"
      else
         compute ltitle="OLS with Panel HAC Covariance Matrix"
	 }
   do t=1,tp
      ewise xt(i,j)=xv=%eqnxvector(0,identries(i)(t)),xv(j)*%resids(identries(i)(t))
      do l=0,lags
         if (t+l>tp)
            break
         if l==0 {
            if method==2
               compute xvx=xvx+tr(xt)*xt
            else
               compute xvx=xvx+tr(xt)*%fill(np,np,1.0)*xt
         }
         else {
            ewise xtlag(i,j)=xv=%eqnxvector(0,identries(i)(t+l)),xv(j)*%resids(identries(i)(t+l))
            if method==2
               compute xlx=tr(xtlag)*xt
            else
               compute xlx=tr(xtlag)*%fill(np,np,1.0)*xtlag
            compute xvx=xvx+wt(l)*(xlx+tr(xlx))
         }
      end do l
   end do t
}
compute pcsexx=%xx*xvx*%xx
*
linreg(create,lastreg,title=ltitle,$
   covmat=pcsexx,form=chisquare)
end
