* SHAZAM PROCEDURE FOR ESTIMATION OF THE RANDOM COEFFICIENT MODEL
*           FOR POOLED TIME-SERIES CROSS-SECTION DATA
SET NOECHO
PROC RANDCOEF
* The assumptions are cross-section heteroskedasticity but no
* correlation between cross-sections and no serial correlation.
* Estimation is by GLS. The estimation procedure is a 2-step
* procedure. The first step is OLS estimation of the individual
* cross-sections. These results are used to form the GLS estimator
* in the second step.
* References: Judge, Griffiths, Hill, Lutkepohl, and Lee, 1985,
*                 The Theory and Practice of Econometrics, 2nd Edition,
*                 Chapter 13.5, pages 538+.
*             Dielman, Terry, E., Pooled Cross-Sectional and Time
*                 Series Data Analysis, 1989, Chapter 5.
*=================== Procedure variables =================================
*  The data must be arranged with the time series for the first 
*  cross-section followed by the time series for the second 
*  cross-section etc. That is, in the form required by the SHAZAM
*  POOL command. The panels must NOT contain missing observations but
*  the time series lengths may differ across cross-section units.
*  Input variables must be defined as follows:
*   NCROSS: number of cross section units (N)
*   NTIME:  an array of order (Nx1) with the number of observations 
*           in each cross-section.
*   BEG:    the position of the first observation (eg. BEG: 1)
*   DEPVAR: dependent variable
*   INDEPS: list of independent variables. The list must first give the
*           variables with fixed coefficients (if desired) followed by
*           the variables with random coefficients.
*           A constant is automatically included and is always assumed
*           to be a random coefficient.
*   NFIX:   The number of fixed coefficients. If all coefficients are
*           random then specify  NFIX: 0
*  Output variables must be defined as follows:
*   RESID:  Transformed residuals
*   COEF:   Estimated coefficients. If COEF: BETA is specified then
*           on output BETA will contain the mean parameters and
*           BETA1, BETA2, ... will contain the estimated parameters
*           for cross-section 1, cross-section 2 etc.
*   STDERR: Estimated standard errors. If STDERR: STD is specified then
*           on output STD will contain the estimated std.errors for the
*           mean parameters and STD1, STD2, ... will contain the 
*           std.errors for the parameter estimates for cross-section 1, 
*           cross-section 2 etc.
*========================================================================
SET NODOECHO NOOUTPUT
GEN1 N2_=[BEG]-1
DO #=1,[NCROSS]
GEN1 N2_=N2_+[NTIME]:#
ENDO
* OLS estimation
SAMPLE [BEG] N2_
OLS [DEPVAR] [INDEPS]
GEN1 NOBS_=$N
GEN1 K_=$K
GENR ONE_=1
GENR [RESID]=0
SAMPLE 1 K_
GENR [COEF]=0
* Initialize working arrays
GEN1 KR_=K_-[NFIX]
IF1(KR_.LE.0) KR_=1
GEN1 KRB_=K_-KR_+1
DIM SB_ KR_ KR_ SUM_ KR_ SSB_ KR_ KR_ G_ 1 WEIGHT_ NOBS_ 
* Estimate each cross-section by OLS
GEN1 N2_=[BEG]-1
DO #=1,[NCROSS]
* Check that enough observations
SET NOWARN
IF(([NTIME]:#).LT.K_)
STOP
SET WARN
GEN1 N1_=N2_+1
GEN1 N2_=N1_+[NTIME]:#-1
SAMPLE N1_ N2_
OLS [DEPVAR] [INDEPS] / COEF=[COEF]# COV=ECOV#_
GEN1 SIG2#_=$SIG2
GENR WEIGHT_=1/SIG2#_
COPY [COEF]# BRAN_ / FROW=KRB_;K_ TROW=1;KR_
* The next statement is needed because ECOV is symmetric.
MATRIX ECOV#_=ECOV#_
COPY ECOV#_ CRAN_ / FROW=KRB_;K_ FCOL=KRB_;K_ TROW=1;KR_ TCOL=1;KR_
MATRIX SB_=SB_+BRAN_*BRAN_'
MATRIX SUM_=SUM_+BRAN_
MATRIX SSB_=SSB_+CRAN_
MATRIX IECOV#_=INV(CRAN_)
ENDO
* Estimation assuming cross-section heteroskedasticity 
* and equal coefficients for all cross-sections (Judge, p.543)
SAMPLE [BEG] N2_
OLS [DEPVAR] [INDEPS] / COEF=BETAW_ WEIGHT=WEIGHT_
GEN1 LLF0_=$LLF
* Construct test (Judge p. 542-543, Dielman, p. 99-102)
COPY BETAW_ BMAN_ / FROW=KRB_;K_ TROW=1;KR_
DO #=1,[NCROSS]
COPY [COEF]# BRAN_ / FROW=KRB_;K_ TROW=1;KR_
MATRIX BRAN_=BRAN_-BMAN_
MATRIX G_=G_+BRAN_'IECOV#_*BRAN_
ENDO
SAMPLE 1 1
GEN1 DF_=KR_*([NCROSS]-1)
DISTRIB G_ / TYPE=CHI CDF=PVALUE_ DF=DF_
GEN1 PVALUE_=1-PVALUE_
** Test for H0: the individual coefficients are NOT random.
** Under H0 the statistic G is asy. Chi-square with d.f.=DF
* Under the alternative the coefficients may be regarded as fixed
* and different or as random and different. Thus this test may be 
* regarded as an indirect test for randomness.
FORMAT(F15.4,F12.0,F16.4)
PRINT G_ DF_ PVALUE_ / FORMAT
* Construct the DELTA matrix (the covariance matrix for the
* the coefficients)  (Judge p.542)
MATRIX SB_=(SB_-SUM_*SUM_'/[NCROSS])/([NCROSS]-1) 
MATRIX DRAN_=SB_ - SSB_/[NCROSS]
MATRIX DET_=DET(DRAN_)
* Make an adjustment if the matrix is negative definite (Judge p. 542)
SET NOWARN
IF(DET_.LE.0) 
MATRIX DRAN_=SB_
SET WARN
DIM DELTA_ K_ K_ VARCOV_ K_ K_ 
COPY DRAN_ DELTA_ / TROW=KRB_;K_ TCOL=KRB_;K_ FROW=1;KR_ FCOL=1;KR_
* Construct the GLS estimator (Judge p. 539-540)
DO #=1,[NCROSS]
MATRIX VARCOV_=VARCOV_+INV(DELTA_+ECOV#_)
ENDO
MATRIX VARCOV_=INV(VARCOV_)
DO #=1,[NCROSS]
MATRIX [COEF]=[COEF]+INV(DELTA_+ECOV#_)*[COEF]#
ENDO
MATRIX [COEF]=VARCOV_*[COEF]
*----------------------------------------------------------------------
* Calculate the transformed residuals from GLS estimation
GEN1 NT_=[NTIME]:1
DIM X_ NT_ K_ Y_ NT_ P_ NT_ NT_
GEN1 NTB_=NT_
GEN1 N2_=[BEG]-1
DELETE DET_
GEN1 DET_=0
DO #=1,[NCROSS]
GEN1 NT_=[NTIME]:#
GEN1 N1_=N2_+1
GEN1 N2_=N1_+NT_-1
* If unequal observations then allocate arrays appropriately
SET NOWARN
IF(NT_.NE.NTB_)
DELETE X_ Y_ P_
IF(NT_.NE.NTB_)
DIM X_ NT_ K_ Y_ NT_ P_ NT_ NT_
SET WARN
COPY [INDEPS] ONE_ X_ / FROW=N1_;N2_ FCOL=1;K_ TROW=1;NT_ TCOL=1;K_
COPY [DEPVAR] Y_ / FROW=N1_;N2_ TROW=1;NT_
MATRIX P_=FACT(X_*DELTA_*X_'+SIG2#_*IDEN(NT_))
MATRIX Y_=P_*Y_-P_*X_*[COEF]
COPY Y_ [RESID] / TROW=N1_;N2_ FROW=1;NT_
MATRIX Y_=LOG(DIAG(P_))
STAT Y_ / SUMS=G_ BEG=1 END=NT_
GEN1 DET_=DET_+G_
ENDO
* Get the sum of the squared residuals
SAMPLE [BEG] N2_
GENR ONE_=[RESID]*[RESID]
STAT ONE_ / SUMS=G_
GEN1 LLF_=-NOBS_*(1+LOG(2*$PI))/2 + DET_ - NOBS_*LOG(G_/NOBS_)/2
** Value of the log-likelihood function
PRINT LLF_
SAMPLE 1 1
GEN1 LRTEST_=2*(LLF_-LLF0_)
GEN1 DF_=KR_*(KR_+1)/2
DISTRIB LRTEST_ / TYPE=CHI CDF=PVALUE_ DF=DF_
GEN1 PVALUE_=1-PVALUE_
** Likelihood ratio test for H0: the individual coefficients are NOT random.
** Under H0 the statistic LRTEST is asy. Chi-square with d.f.=DF
FORMAT(F15.4,F12.0,F16.4)
PRINT LRTEST_ DF_ PVALUE_ / FORMAT
*---------------------------------------------------------------------
** GLS estimation of the mean of the model parameters
** (the intercept parameter is last).
FORMAT(F15.0)
PRINT NOBS_ / FORMAT
SAMPLE 1 K_
MATRIX [STDERR]=SQRT(DIAG(VARCOV_))
MATRIX TRATIO_=[COEF]/[STDERR]
FORMAT(3F15.4)
PRINT [COEF] [STDERR] TRATIO_ / FORMAT
*---------------------------------------------------------------------
* Now get the individual estimates using the method in 
* Kadiyala and Oberhelman, Response Predictions in Regressions on
* Panel Data, Communications in Statistics - Theory and Methods,
* Vol. 11, 1982, pages 2699-2714.
* This is equivalent to the estimator described in Judge, p. 541.
* Kadiyala and Oberhelman note that this estimator should be more
* efficient than OLS applied to each cross-section when
* N is "large" and T is "small".
** Estimates of the coefficients for each cross-section
DO #=1,[NCROSS]
MATRIX MAT_=INV(DELTA_+ECOV#_)
MATRIX W1_=DELTA_*MAT_
MATRIX W2_=ECOV#_*MAT_
MATRIX [COEF]#=W1_*[COEF]#+W2_*[COEF]
MATRIX ECOV#_=ECOV#_*W1_'+W2_*VARCOV_*W2_'
MATRIX [STDERR]#=SQRT(DIAG(ECOV#_))
MATRIX TRATIO_=[COEF]#/[STDERR]#
SAMPLE 1 1
GEN1 I_=$DO
GEN1 NOBS_=[NTIME]:#
FORMAT(2F15.0)
PRINT I_ NOBS_ / FORMAT 
SAMPLE 1 K_
FORMAT(3F15.4)
PRINT [COEF]# [STDERR]# TRATIO_ / FORMAT
ENDO
DELETE FORMAT
DELETE / ALL_
SET DOECHO
PROCEND
SET ECHO