*
* @panelfm( options ) start end
* # list of variables (list dependent variable first)
*
* Procedure for multivariate group mean panel FMOLS estimation in
* heterogeneous panels.
*
* 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 [depends upon number of time periods]
*
*   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.
*   TITLE=title of report ["Mean Group Panel FM Estimation"]
*
*   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:
*   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".
*   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 panelfm  start end
type integer start end
*
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 choice  det      2   none   constant trend
option series  smpl     0
option integer lags
option string  title
option rect    *ibetas
option rect    *itstats
*
local vec[series] xvec dxvec dvec rvec
local rect        g
local vec[series] ytilde
local series      u ehat
local vect[int]   datavec rhs
local integer     nsecs
local integer     mlag
local integer     i j k m row
local series      time
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 report      preport
local vect[int]   detvars
local vect        lw alepht
local string      descript ltitle
local integer     startl endl regnobs starti endi
local series      icount
local real        pt nt
local integer     tmax tend isec
local rect        cmomx
local vect        prjcoe
local real        sigma1
local series      ydagger
local string      ltitle
local equation    baseeq
*
* Globals
*
dec rect %%ibetas %%istderrs %%itstats
dec vect[symm] %%ixx
*
* Check for valid settings
*
if %defined(block) {
   disp "###@PANELFM(BLOCK=nnn). No longer supported. Use PANEL calendar scheme"
   return
}
if det==1 {
   disp "###@PANELFM 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
*
* 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 xvec(m) dxvec(m) dvec(m+1) rvec(m+1) ytilde(m+1) bvec(m)
dim tfm(%nreg)
*
* 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(entry=1.0,time=-1.0,smpl=%valid(u)) datavec(K) start end ytilde(k)
   else
      set ytilde(k) start end = %if(%valid(u),datavec(k){0},%na)
end do k
*
* Zero out the matrices for sums across individuals
*
compute pfmbsum1 = %zeros(%nreg,1)
compute pfmbsum2 = %zeros(%nreg,1)
compute pfmbsum3 = %zeros(%nreg,1)
compute pfmtsum  = %zeros(%nreg,1)
compute pfmvsum1 = %zeros(%nreg,%nreg)
compute pfmvsum2 = %zeros(%nreg,%nreg)
compute pfmhsum2 = %zeros(%nreg,%nreg)
compute pfmhsum3 = %zeros(%nreg,%nreg)
*
dim %%ibetas(%nreg,nsecs) %%istderrs(%nreg,nsecs) %%itstats(%nreg,nsecs)
dim %%ixx(nsecs)
*
* Figure out how many lags to use for the Bartlett window.
*
if %defined(lags)
   compute mlag=lags
else
   compute mlag=fix(%round((4.0*(tmax/100.0)^(2.0/9.0)),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.
*
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 the preliminary regression of the LHS variable on the RHS
   * variables and the deterministics.
   *
   linreg(noprint) dvec(1) j//2 j//tend ehat
   # dvec(2) to dvec(m+1) detvars
   *
   * Replace rvec(1) with the first stage residuals
   *
   set rvec(1) j//2 j//tend = ehat
   *
   * Compute the long run variance matrix of the residuals and the RHS
   * differenced endogenous variables.
   *
   mcov(lwindow=bartlett,lags=mlag) j//2 j//tend
   # rvec(1) to rvec(m+1)
   compute %cmom=%cmom/%nobs
   *
   * Sweep out the differenced variables and pull out the coefficients of
   * rvec(1) on those. The residual variance is needed later to compute
   * the final covariance matrix, so hang onto that as well.
   *
   compute cmomx =%sweeplist(%cmom,%seq(2,%ncmom))
   compute prjcoe=%xsubmat(cmomx,2,%ncmom,1,1)
   compute sigma1=cmomx(1,1)
   *
   * Create the bias-adjusted dependent variable
   *
   set ydagger j//2 j//tend = dvec(1)-%dot(prjcoe,%xt(dxvec,t))
   *
   * Compute the asymmetrical (current and leads only) version
   * of the long-run variance.
   *
   dim lw(2*mlag+1)
   ewise lw(i)=%if(i>=mlag+1,float(2*(mlag+1)-i)/(mlag+1),0.0)
   mcov(lwform=lw,rmatrix=g) j//2 j//tend
   # rvec(1) to rvec(m+1)
   *
   * Compute the cross-product correction. (We don't scale the MCOV matrix
   * by %nobs since we need to multiply it up when we use it.)
   *
   compute alepht=%xsubmat(g,2,%ncmom,1,1)-$
      %xsubmat(g,2,%ncmom,2,%ncmom)*prjcoe
   *
   * Generate a cross product matrix with the original variables but with
   * ydagger replacing the y
   *
   cmom(noprint) j//2 j//tend
   # ydagger xvec(1) to xvec(m) detvars
   *
   * Adjust the cross products between the dependent variable and the RHS
   * endogenous variables
   *
   do i=1,m
      compute %cmom(1,i+1)=%cmom(1,i+1)-alepht(i)
   end do i
   *
   * Do the regression by sweeping out all RHS variables
   *
   compute cmomx=%sweeplist(%cmom,%seq(2,%ncmom))
   *
   * 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 bfm    = %xsubmat(cmomx,2,%ncmom,1,1)
   compute vfm    = sigma1*%xsubmat(cmomx,2,%ncmom,2,%ncmom)
   compute hfm    = inv(vfm)
   compute dfm    = inv(%sqrt(%diag(%xdiag(vfm))))
   *
   * Compute t-stat vs the hypothesized vector
   *
   ewise tfm(i)   = (bfm(i)-%if(i<=m,bvec(i),0.0))*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
   *
   * Save the individual coefficients, standard errors and t-statistics
   * in the global variables.
   *
   compute %psubmat(%%ibetas,1,isec,bfm)
   compute %psubmat(%%itstats,1,isec,tfm)
   compute %psubmat(%%istderrs,1,isec,%sqrt(%xdiag(vfm)))
   compute %%ixx(isec)=vfm
   *
   * Update the infobox
   *
   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 global variables for overall estimator
*
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
*
compute %stderrs=%sqrt(%xdiag(%xx))
compute %ngroup =nsecs
*
* 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 FM 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"
if %defined(lags)
   report(use=preport,atrow=5,atcol=2) "Kernel Width" lags
else
   report(use=preport,atrow=5,atcol=2) "Kernel Width" "Computed"

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

