// CORNWELL SCHMIDT AND SICKLES (1990) RANDOM EFFECTS
// CODED BY PHILL WHEAT, SENIOR RESEARCH FELLOW, UNIVERISTY OF LEEDS
// P.E.WHEAT@ITS.LEEDS.AC.UK
// 31/03/2015
// CONVERTED TO A STATA PROGRAM, AND OTHER MINOR CHANGES BY ALEX STEAD, m2i RESEARCH FELLOW IN EFFICIENCY ANALYSIS, UNIVERSITY OF LEEDS
// A.D.STEAD@LEEDS.AC.UK
// 07/10/2016

/* THIS CODE CAN BE USED TO ESTIMATE THE TIME VARYING RANDOM EFFECTS MODEL OF 
CORNWELL, SCHMIDT AND SICKLES (1990) Journal of Econometrics 46 (1990) 185-200.
THE OUTPUT IS A VECTOR, BGLS, OF PARAMETER ESTIMATES, A COVARIANCE MATRIX 
WITH RESPECT TO THE ESTIMATES, GLSVC, AND VARIABLES, INEFF AND EFF, BEING THE 
INEFFICIENCY AND EFFICIENCY PREDICTIONS RESPECTIVELY. 

THIS CODE WAS DEVELOPED ORIGINALLY FROM CODE SOURCED FROM THE WEBSITE OF ROBIN SICKLES (WONG SONG, 2002) AND CODED
IN LIMDEP IN 2012. SUBSEQUENT TO THIS CODING, SICKLES RELEASED A REVISED VERSION OF THE CODE (Wong Sung, 2014)
CORRECTING SOME ERRORS IN THE ORIGINAL. THIS CODE INCLUDES THESE CHANGES. AS WELL AS SOME DEGREES OF 
FREEDOM CORRECTION CHANGES, THE MAIN CHANGE IS LINES 107 TO 113 (THIS CODE). SONG HAS CORRECTED THE CODE
TO ALLOW THE WITHIN ESTIMATOR TO BE USED. IN THE LIMDEP CODE AT THE TIME WE ADOPTED A POOLED ESTIAMTOR AS
THE WITHIN GAVE VERY STRANGE RESULTS (BECAUSE THE TIME VARYING RANDOM COMPOENNTS WERE NOT EXPLICITLY 
TAKEN INTO ACCOUNT). FOR COMPARISION I HAVE RETAINED THE POOLED CODE BUT COMMENTED IT OUT OF THE SCRIPT.
IF THIS IS REINSTATED (AND THE WITHIN ELEMENTS COMMENTED OUT) THEN THE CODE REPLICATES THE LIMDEP
CODE (WITH THE EXCEPTION OF A VERY SMALL EFFECT OF SOME DEGREES OF FREEDOM CHANGE ON LINE 97) */

quietly{
set matsize 800

program cssre, eclass
version 11
if !replay() {
syntax varlist(numeric ts) [if] [in] [, Level(cilevel)]
marksample touse

gettoken depvar indepvars : varlist
_fv_check_depvar `depvar'

fvexpand `indepvars'
local cnames `r(varlist)'

tempname qtime nreg colxm colwmat rowind zcoln qmat tempold c d f g j k rowind1 temp1 kend temp2 pmat  rowp mmat intwith int2with bwith  int3with  reswith ssgwint ssgtop ssgw intbet int2bet bbet resp pw respoo deltamat f1 c1 c2 c3 c4 c5 omega d1 d2 d3 bgls glsvc a1 a2 a3 a4 a5 a6 a7 a8 a9 observ

capture drop alpha eff ineff
gen cons=1
xtset
gen panelvar=`r(panelvar)'
gen timevar=`r(timevar)'
replace timevar=timevar-r(tmin)+1
gen timevar_sq=timevar^2
local qtime "cons timevar timevar_sq"
local cons "cons"

mkmat `depvar', matrix(ymat) //y vector
mkmat `indepvars', matrix(xmat)	//x matrix
mkmat `qtime', matrix(wmat)	//q matrix
mkmat `cons' `indepvars', matrix(commat)	//x and q matrix


quietly tab panelvar, matcell(indicat)	//how many obs for each firm?
scalar nreg=rowsof(xmat)	//number of observations
scalar colxm=colsof(xmat)	//number of xvars
scalar colwmat=colsof(wmat)	//number of time vars
scalar rowind=rowsof(indicat)	//number of firms
scalar zcoln=rowind*colwmat		//number of time parameters

matrix qmat=J(nreg,zcoln,0)		//n by no. of time parameters matrix of zeros

scalar tempold=0

local i=1
scalar k=1
scalar j=1
scalar rowind1=rowind+1

while (`i'<rowind1) {			//for firm i ...
scalar j=(colwmat*(`i'-1))+1	//define j as 3*i-2
scalar temp1=indicat[`i',1]		//take the number of observations for that firm
scalar kend=k+temp1-1			//
matrix temp2=wmat[k..kend,1..colwmat]	
matrix qmat[k,j]=temp2			//
scalar k=k+temp1
local i=`i'+1
}


// P transform of Q matrix (note Q matrix is not the Q of standard panel data
// instead Q is the diagonally stacked W matrix, M is the Q transform (group means - trends))


matrix pmat=qmat*invsym(qmat'*qmat)*qmat'
scalar rowp=rowsof(pmat) 
matrix mmat=I(rowp)-pmat

// compute within estimator

matrix intwith= xmat'*mmat*xmat
matrix intwith=(intwith+intwith')/2
matrix int2with=invsym(intwith)
matrix bwith=int2with*xmat'*mmat*ymat
matrix int3with=mmat'*xmat
matrix reswith=mmat'*ymat-(int3with*bwith)
matrix ssgwint=reswith' * reswith
scalar ssgtop=ssgwint[1,1]
scalar ssgw=ssgtop/(rowp-(rowind*colwmat)-colxm) // new sickles code indicates instead of (rowp -rowind-colxm) note previous annotation was (rowp -colwat-colxm)
// scalar ssgw=ssgtop/(rowp-rowind-colxm)

// compute between estimator

matrix intbet=xmat'*pmat*xmat
matrix int2bet=invsym(intbet)
matrix bbet=int2bet*xmat'*pmat*ymat

// compute pooled residuals (using within para vector - *this is a change to the LIMDEP code reflecting an update from Sickles)

matrix resp=ymat-xmat*bwith
matrix pw=wmat*invsym(wmat'*wmat)*wmat'
// matrix respoo=resp-pw*resp

// test for previous model - ONLY UNCOMMENT OUT IF WISH TO REPLICATE LIMDEP RESULTS

matrix respoo=ymat-commat*invsym(commat'*commat)*commat'*ymat

// Now create delatamat


matrix deltamat=J(colwmat,colwmat,0)
local i=1
scalar k=1
while (`i'<rowind1) {
scalar temp1=indicat[`i',1]
scalar kend=k+temp1-1
matrix c=wmat[k..kend, 1..colwmat]
matrix d=respoo[k..kend, 1..1]
matrix f=invsym(c'*c)*c'*d*d'*c*invsym(c'*c)-ssgw*invsym(c'*c)
matrix g=f*temp1
matrix deltamat=deltamat+g
scalar k= k+temp1
local i=`i'+1
}

scalar f1=1/(rowp)
matrix deltamat=f1*deltamat


// create omega

matrix c1=I(rowind)
matrix c2=c1#deltamat
matrix c3=I(rowp)
matrix c4=c3*ssgw
matrix c5=qmat*c2*qmat'
matrix omega=c5+c4

// estimate beta 


matrix d1=invsym(omega)
matrix d2=commat'*d1*commat
matrix d3=invsym(d2)
matrix bgls= d3*commat'*d1*ymat
matrix bgls=bgls'
matrix glsvc=d3

// Sickles' code now computes the residuals based on subtracting both the 
// xs and the ws. But it makes no difference to the efficiencies becuase they are relative.
// to verify run //out code below

matrix a1=J(rowp,1,0)
scalar a2=colwmat+1
scalar a3=colwmat+colxm
//matrix bglsx=bgls[a2..a3,1..1]
//matrix a4=ymat-xmat*bglsx
matrix bgls=bgls'
matrix a4=ymat-commat*bgls
matrix bgls=bgls'

local i=1
scalar k=1
while (`i'<rowind1) {
scalar a5=indicat[`i',1]
scalar kend=k+a5-1
matrix a6=a4[k..kend, 1..1]
matrix a7=wmat[k..kend, 1..colwmat]
matrix a8=invsym(a7'*a7)*a7'*a6
matrix a9=a7*a8
matrix a1[k,1]=a9
scalar k= k+a5
local i=`i'+1
}

// Now compute the efficiency scores

// New matrix for reordering

svmat a1, name(alphait)
egen observ = fill(1 2)
sort timevar panelvar
by timevar: egen minyear=min(alphait)
gen ineff=alphait-minyear
sort panelvar timevar
gen eff=exp(-ineff)      

// Display output table
ereturn post bgls glsvc
display ""
display "Time-varying random-effects model (CSS GLS)"
display ""
ereturn display
drop cons panelvar timevar timevar_sq observ minyear
rename alphait alpha
}
end
}