*
* @Swamy( options )  depvar start end
* # list of regressors
*
* Computes a GLS matrix weighted estimator for a panel data set. The
* procedure MeanGroup.src does a simple weighted average of the
* coefficient vectors.
*
* Parameters:
*   depvar      dependent variable
*   start end   range to estimate
*
* Note: This procedure uses <<start>> and <<end>> only to determine the
* range of individuals to use in the estimation. It will use all valid
* observations for each individual within that range. Use the SMPL
* option if you want to skip observations for reasons other than lack of
* data.
*
* Options:
*   DELTA =covariance matrix of beta(i)-beta [estimated from data]
*
*   SMPL  =SMPL series [none]
*
*   TITLE =title for output ["Swamy Random Coefficients"]
*   [PRINT]/NOPRINT
*   DEFINE=equation to define [none]
*
* Reference:
*   Swamy(1970), "Efficient Inference in a Random Coefficient Regression
*    Model", Econometrica, vol 38, no 2, pp 311-323.
*
* Variables Defined:
*   %BETA, %XX are for the central estimated model
*   %NGROUP  is the number of individuals in the sample
*   %%IBETAS is an %NREG x %NGROUP matrix with the individual coefficients
*            in the columns.
*   %%IXX    is an %NGROUP vector of covariance matrices for the individual
*            coefficients.
*
* Revision schedule:
*  10/2002 Converted by procedure by Estima
*  04/2011 Added TITLE option
*  05/2012 Added PRINT option. Added DELTA option. Modernized handling of SMPL.
*          Added calculating and saving of individual coefficient vectors
*          and covariance matrices.
*
procedure Swamy depvar start end
type series depvar
type integer start end
*
option equation *define
option series   smpl
option switch   print		1
option string   title
*
local vect[integer] 	reglist
local integer 			startl endl starti endi regnobs
local integer 			i j k
local integer 			nindiv ncoeff
local real           pt nt
local vect[symm] 		xxmats
local vect[vect] 		betas
local vector 			sigmas
local symm 				xxgls delta idelta
local vect 				betagls hb
local equation 		paneleqn
local string         ltitle
local series         icount u
local integer        tmax tend ifill
*
* Globals
*
dec rect %%ibetas
dec vect[symm] %%ixx
*
enter(varying) reglist
*
inquire(regressorlist) startl<<start endl<<end
# depvar reglist
equation paneleqn depvar
# reglist
compute ncoeff=%eqnsize(paneleqn)
*
* Do a silent LINREG to get the SMPL for the data
*
linreg(noprint,smpl=smpl,equation=paneleqn) * start end u
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 nindiv=fix(nt)
*
* Dimension everything.
*
dim xxmats(nindiv) betas(nindiv) sigmas(nindiv)
dim xxgls(ncoeff,ncoeff) betagls(ncoeff) delta(ncoeff,ncoeff)
dim %%ixx(nindiv) %%ibetas(ncoeff,nindiv)
*
* First pass, run OLS regressions over each individual. Save %xx, %beta
* and %seesq.
*
compute delta=%const(0.0)
compute betagls=%const(0.0)
compute ifill  =0
do i=starti,endi
   if icount(i)==0
      next
   linreg(smpl=smpl,noprint,equation=paneleqn) * i//1 i//tend
   compute ifill=ifill+1
   compute xxmats(ifill)=%xx
   compute betas(ifill)=%beta
   compute sigmas(ifill)=%seesq
   compute delta=delta+%outerxx(%beta)
   compute betagls=betagls+%beta
end do i
*
* Estimate the covariance matrix of beta's
*
ewise delta(j,k)=1.0/(nindiv-1)*(delta(j,k)-(1.0/nindiv)*betagls(j)*betagls(k))
*
* Second pass. Compute the GLS covariance matrix. While we're at it,
* replace xxmats with the precision. Since we need it for computing
* individual coefficients, also save the precision of the estimates
* without the delta.
*
compute xxgls=%const(0.0)
do i=1,nindiv
   compute %%ixx(i)=inv(sigmas(i)*xxmats(i))
   compute xxmats(i)=inv(delta+sigmas(i)*xxmats(i))
   compute xxgls=xxgls+xxmats(i)
end do i
compute xxgls=inv(xxgls)
*
* Final pass. Compute the grand coefficient matrix.
*
compute betagls=%const(0.0)
do i=1,nindiv
   compute betagls=betagls+xxgls*xxmats(i)*betas(i)
end do i
if %defined(title)
   compute ltitle=title
else
   compute ltitle="Swamy Random Coefficients"
*
linreg(title=ltitle,print=print,create,equation=paneleqn,$
   define=define,coeffs=betagls,covmat=xxgls,form=chisqr) *
*
* Define the individual coefficients and covariance matrices.
*
compute idelta=inv(delta)
do i=1,nindiv
   compute hb=%%ixx(i)*betas(i)+idelta*betagls
   compute %%ixx(i)=inv(%%ixx(i)+idelta)
   compute %psubmat(%%ibetas,1,i,%%ixx(i)*hb)
end do i
*
compute %ngroup=nindiv
*
end

