*
* @panelfm( options ) start end
* # list of variables (list dependent variable first)
*
* Procedure for multivariate group mean panel FMOLS 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.
*
*   lags=number of lags to use in the Bartlett kernel [depends upon number of time periods]
*
* 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 Theory and Monte Carlos:
*
*   Pedroni, Peter (2000) "Fully Modified OLS for Heterogeneous Cointegrated
*     Panels," Advances in Econometrics, Vol. 15, 93-130, NONSTATIONARY PANELS, PANEL
*     COINTEGRATION AND DYNAMIC PANELS, JAI Press.
*
*   Working paper version: Pedroni, Peter (1996), Indiana University
*     Working Papers in Economics, No. 96-020.
*
*   For general discussion and application:
*
*   Pedroni, Peter (2001) "Purchasing Power Parity Tests in Cointegrated Panels,"
*     Review of Economics and Statistics, 83, 727-731.
*
* Revision Schedule:
*   Originally written by Peter Pedroni, Williams College: ppedroni@williams.edu
*   Includes revisions made by Prof. Pedroni to his original program that allow for
*   common time effects in unbalanced panels, automatic selection of bandwidth for
*   kernels, and automatic adjustment of display for number of regressors
*
*   06/2004 Converted to procedure by Tom Doan, Estima
*   09/2005 Changed: compute term1 = tr(G21+V21)
*                To: compute term1 = tr(G12+V21)
*           to match Pedroni's original program
*   04/2008 Code modernized to use REPORT and asymmetric MCOV.
*           LAGS, 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 panelfm  startp endp
type integer startp endp
*
option integer block
option switch  tdum     0
option vector  bveco
option choice  average  1   simple sqrt   precision
option choice  print    3   none   short  full
option integer lags     0
option rect    *ibetas
option rect    *itstats
*
local integer tperiods jperiods start end
*
local real        maxt mint
local vect        W
local vec[real]   itemp1 O21 G12 V21
local rect[real]  gammastar gXX O22 G22 V22
local vec[series] xvec dxvec dvec rvec
local rect        g
local vec[series] ytilde
local series      u ehat
local vect[int]   datavec rhs savecal
local integer     nsecs
local integer     mlag
local integer     i j k m row
local symm        v omega
local rect        term1 term2 itemp2
local series      ystar tvec
local real        o11
local vec         bvec
local vec[label]  labelvec
local integer     unbal
local vect        pfmbsum1 pfmbsum2 pfmbsum3 pfmtsum bfm tfm
local symm        pfmvsum1 pfmvsum2 pfmvsum3 pfmhsum2 pfmhsum3 vfm hfm dfm
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 xvec(m) dxvec(m) dvec(m+1) rvec(m+1) ytilde(m+1) tfm(m) bvec(m)



                        *** MAIN SOURCE CODE ***

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

*
* Zero out the matrices for sums across individuals
*

compute pfmbsum1 = %zeros(m,1)
compute pfmbsum2 = %zeros(m,1)
compute pfmbsum3 = %zeros(m,1)
compute pfmtsum  = %zeros(m,1)
compute pfmvsum1 = %zeros(m,m)
compute pfmvsum2 = %zeros(m,m)
compute pfmhsum2 = %zeros(m,m)
compute pfmhsum3 = %zeros(m,m)

infobox(action=define,progress,lower=1,upper=NSecs) "Computing Panel Statistics"
do J=firstInd,lastInd
   *
   * Figure out how many lags to use
   *
   if lags>0
      compute mlag=lags
   else
      compute mlag = fix(%round((4.0*(Tvec(J)/100.0)**(2.0/9.0)),0))
   *
   * 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 the preliminary regression of the LHS variable on the RHS variables and
   * the CONSTANT.
   *
   linreg(noprint) dvec(1) 2 jperiods ehat
   # dvec(2) to dvec(m+1) constant
   *
   * Replace rvec(1) with the first stage residuals
   *
   set rvec(1) =  ehat
   *
   * Compute the variance and long-run variance matrices
   *
   vcv(matrix=v,noprint) 2 jperiods
   # rvec(1) to rvec(m+1)
   mcov(lwindow=bartlett,lags=mlag,noprint) 2 jperiods
   # rvec(1) to rvec(m+1)
   compute Omega =%cmom/%nobs
   *
   * Pull out the information that we need from those
   *
   compute O11   =Omega(1,1)
   compute O21   =%xsubmat(Omega,2,m+1,1,1)
   compute V21   =%xsubmat(V,2,m+1,1,1)
   compute O22   =%xsubmat(Omega,2,m+1,2,m+1)
   compute V22   =%xsubmat(V,2,m+1,2,m+1)
   *
   * Do the asymmetrical window using only the current and leads.
   * (The LWFORM vector on MCOV has leads first, then 0, then lags.)
   *
   dimension w(2*mlag+1)
   ewise w(i)=%if(i<=mlag+1,0.0,1.0-(i-mlag-1)/(mlag+1.0))
   mcov(lwform=w,lags=mlag,noprint,rmatrix=g) 2 jperiods
   # rvec(1) to rvec(m+1)
   compute g=g/%nobs
   compute g22=%xsubmat(g,2,m+1,2,m+1)
   compute g12=%xsubmat(g,2,m+1,1,1)
   *
   * Compute the modified estimator
   *
   set ystar 2 jperiods = dvec(1)-%dot(tr(o21)*inv(o22),%xt(dxvec,t))
   compute term1 = tr(G12+V21)
   compute term2 = tr(o21)*inv(o22)*(tr(G22)+V22)
   compute gammastar = term1 - term2
   linreg(noprint) ystar 2 jperiods
   # xvec(1) to xvec(m) constant
   *
   * 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 gxx    = %xsubmat(%xx,1,m,1,m)
   compute itemp1 = %xsubvec(%beta,1,m)
   compute itemp2 = (jperiods-1.0)*gammastar*gXX
   compute bfm    = itemp1 - tr(itemp2)
   compute vfm    = o11*gxx
   compute hfm    = inv(vfm)
   compute dfm    = inv(%sqrt(%diag(%xdiag(vfm))))
   *
   * Compute t-stat vs the hypothesized vector
   *
   ewise tfm(i)   = (bfm(i)-bvec(i))*dfm(i,i)
   *
   * Keep track of several combined measures:
   *  pfmbsum1 is the arithmetic average
   *  pfmbsum2 is weighted by the square root of the diagonals of the precision
   *  pfmbsum3 is weighted by the precision
   *  pfmvsum1 is the sum of the variance matrices
   *  pfmvsum2 is the sum of the #2 rescaled variance matrices
   *  pfmhsum2 is the sum of the #2 weights
   *  pfmhsum3 is the sum of the #3 weights
   *  pfmtsum  is the sum of the t-statistics
   *
   compute pfmbsum1 = pfmbsum1  + bfm
   compute pfmbsum2 = pfmbsum2  + dfm * bfm
   compute pfmbsum3 = pfmbsum3  + hfm * bfm
   compute pfmvsum1 = pfmvsum1  + vfm
   compute pfmvsum2 = pfmvsum2  + dfm * vfm * dfm
   compute pfmhsum2 = pfmhsum2  + dfm
   compute pfmhsum3 = pfmhsum3  + hfm
   compute pfmtsum  = pfmtsum   + tfm
   gset betas  j j = bfm
   gset tstats j j = tfm
   infobox(action=modify,current=j)
end do j
infobox(action=remove)
*
* Standardize the estimates as appropriate
*
compute pfmtsum  = (1.0/sqrt(NSecs))*pfmtsum
compute pfmbsum1 = (1.0/NSecs)*pfmbsum1
compute pfmbsum2 = inv(pfmhsum2)*pfmbsum2
compute pfmbsum3 = inv(pfmhsum3)*pfmbsum3
compute pfmvsum1 = (1.0/NSecs)**2*pfmvsum1
compute pfmvsum2 = inv(pfmhsum2)*pfmvsum2*inv(pfmhsum2)
compute pfmvsum3 = inv(pfmhsum3)
*
* Define variables
*
compute %tstats=pfmtsum
if average==1
   compute %beta=pfmbsum1,%xx=pfmvsum1
else
if average==2
   compute %beta=pfmbsum2,%xx=pfmvsum2
else
if average==3
   compute %beta=pfmbsum3,%xx=pfmvsum3

***********************************************************************

    *** DISPLAY INDIVIDUAL AND PANEL RESULTS ***

report(action=define)
report(atrow=1,atcol=1,span) "Mean Group Panel FM 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"
if lags==0
   report(atrow=5,atcol=2) "Kernel Width" "Computed"
else
   report(atrow=5,atcol=2) "Kernel Width" lags

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) pfmtsum(k)
   else
      report(row=new,atcol=2)         labelvec(k+1) %beta(k) pfmtsum(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







