options double crt;
name lsdvc 'LSDVc:  bias-corrected fixed effects for dynamic panel';

? Estimates a dynamic panel model with fixed effects
?  y_it = gamma*y_i,t-1 + beta*x_it + alpha_i + e_it
? using the bias correction due to
? Jan Kiviet, Journal of Econometrics 68, 1995, pp.53-78.
?   Contains an example with the "penngrow" benchmark panel data.

? by Clint Cummins 5/2000

? Read balanced data here and set proper sample for regression
const n,94 t,6;
set nt = n*t;
smpl 1,nt;
read(file='penngrow.txt') country year y x;
freq(panel,n=n,t=t,id=country,time=year) ;
select year>1960;
ylag = y(-1);  ? transformation to get rid of explicit lag

? If you want to use the value of blsdvc outside of the Proc,
? make it a global variable:
mmake blsdvc 0 0;  ? the actual dimension given here doesn't matter

? For the penngrow model/dataset above, the results are:
?            Within/LSDV     AH-L       LSDVC
?   y(-1)      .720369     .686996     .800760
?   (s.e.)    (.023677)   (.064469)   (.026718)
?   x          .165589     .148424     .156264
?             (.019332)   (.021146)   (.019472)
? These coefficients agree to at least 6 digits with Gauss
? results based on JOSIM1.GAU.  (The SEs were not checked).

LSDVc y ylag x;   ? Call the proc with the variables.
                  ? Leave out the intercept C.
                  ? Don't use any explicit lags - this makes it
                  ? easier for the Proc to use the DOT command.
print blsdvc;  ? print variable outside the Proc

?=====================
Proc LSDVc yxs;

? Computes the bias corrected fixed effects estimator.
? (also computes the AH-L estimator)
? Input arguments:
?   Dependent variable,
?   Lagged dependent variable (without explicit lag),
?   Other RHS exogenous variables.
?   Do not include the constant term C
?   (it is implicit in these fixed effects models).
? Outputs:
?   blsdvc , if it is created as a global variable before
?            calling this proc.
? Notes:
? 1. Adapted from JOSIM1.GAU, by Ruth Judson, Feb 1996
?    (downloaded from Ann L. Owen's web page, 8/1999)
?      http://academics.hamilton.edu/economics/aowen
?    The citation for the discussion paper is:
?    Judson, Ruth, and Ann L. Owen,
?    "Estimating Dynamic Panel Data Models: A Practical Guide
?    for Macroeconomists," FRBG FEDS 1997-3, 1/1996.
?      http://www.bog.frb.fed.us/pubs/feds/1997/199703/199703pap.pdf
?    The published version is in Economics Letters, October 1999.
? 2. The TSP code here has been numerically tested against Gauss
?    results from the doreg, doah, and docorr sections of JOSIM1.GAU,
?    on the penngrow.txt data, and the answers agree to at least
?    6 digits.
? 3. The initial unbiased estimator used in step 2 below does not
?    have to be AH-L.  It could be something like a 1-step GMM, which
?    would use additional orthogonality conditions.  I haven't
?    attempted to use that at present, because it is more
?    complicated.  For example, there is a choice of the number of
?    lags of the exogenous variables.
? 4. The SEs for the LSDVc are given in Bun and Kiviet, 6/1999:
?      http://www.fee.uva.nl/ke/jfk/bk98a.PDF
?    These are computed here in section 2.  They have not been
?    tested numerically against any other program, but they
?    appear to be in the right range (slightly higher than the
?    LSDV SEs).  Bun and Kiviet's tests indicate the SEs do not
?    perform very well in the rather small samples they considered.
?    It's an area of current research which may show improved results
?    in the future.
? 5. I have only implemented LSDVc for balanced panel data.
?    Kiviet's results could be extended to unbalanced panels,
?    but I haven't attempted it at present.  It would probably not
?    be too tough to compute the LSDVc coefficients, but the SEs
?    would probably be difficult.
? by Clint Cummins 5/2000
?

local smpsav nyx k coefw sesw vcovw s2w ntm1 n tm1 t;

copy @smpl smpsav;
length yxs nyx;
set k = nyx-1;      ? number of coefs (y(-1) , exog)

? 1. Plain (biased) within estimator (LSDV)
panel(notot,nobet,novarc) yxs;
copy @coefw coefw;
copy @sesw sesw;
copy @vcovw vcovw;
copy @s2w s2w;
set ntm1 = @nob;
mat n = nrow(@fixed);
? t = number of time periods the dependent variable
?     is observed (at least t >= 3, since there is a
?     single lag on the RHS, and one more lag is needed for AH-L).
? tm1 = t-1 is the number of observations for each individual
?           in the regression.
set tm1 = ntm1/n;
set t = tm1+1;
print n t;

? 2. Initial unbiased estimator, for gamma, full sample residuals
? using unbiased coefs, and sigeah (est. variance of residuals).
? We use Anderson-Hsiao (AH-L) here, due to its simplicity.
? It would also be possible to use various GMM estimators.

? Create first differenced series
local YL good_AHL d_yxs exog d_xs coefah sesah vcovah s2ah gamma;

list YL yxs(2);
good_AHL = .not.miss(YL(-1));
select(silent) good_AHL;
list(prefix=d_) d_yxs yxs;
local d_yxs;
dot yxs;
   d_. = . - .(-1);  ? first differences of all variables
enddot;
list exog yxs;
list(drop) exog yxs(1) YL;
list(prefix=d_) d_xs exog;
?inst d_yxs INVR YL(-1) d_xs c;
inst d_yxs INVR YL(-1) d_xs;

copy @coef coefah;
copy @ses sesah;
copy @vcov vcovah;
copy @s2 s2ah;
set gamma = coefah(1);
smpl smpsav;

? Create group-demeaned series, for use in making eah (and awbari)
? This is easily done as residuals from a within regression on C.
local dm_yxs i eah sigeah;

list(prefix=dm_) dm_yxs yxs;
local dm_yxs;
set i=0;
dot yxs;
  panel(silent,notot,nobet,novarc) . c;
  dm_. = @resw;
  ?  eah = dm_y - coefah*dm_x
  if (i=0); then;
    eah = @resw;
  else;
    eah = eah - coefah(i)*@resw;
  set i = i+1;
enddot;
? The group-demeaned series can be checked by reproducing
? the original LSDV within results:
? panel(notot,nobet,novarc) dm_yxs;
mat sigeah = (eah'eah)/(ntm1 - k - n);

? 3. bias correction terms
?---------------------------------------------------------------
? The correction consists of 3 terms, and is on p. 64 (Thm 1) in
? Kiviet, Journal of Econometrics 68, 1995
? Use initial consistent estimates from AH-L
?-----------------------------------------------------------------
? NB: Many of these vectors and matrices can be defined/calculated
? less often.  Problems: WBar def, values to use for bias calculation
?-----------------------------------------------------------------
local cvec cmat ici itm1 atmat trcac trcacac atc;

mform(nrow=tm1,ncol=1) cvec=0;
do i=2,tm1;
  set cvec(i) = gamma**(i-2);
enddo;
mform(band,nrow=tm1) cmat = cvec;
mat cmat = tri(cmat)';  ? zero lower triangle, then tranpose

mat ici = sum(cmat);
set itm1 = 1/tm1;
mform(nrow=tm1,type=sym) atmat = itm1;
mat atmat = ident(tm1) - atmat;
mat trcac = tr(cmat'atmat*cmat);
mat trcacac = tr(cmat'atmat*cmat*atmat*cmat);
mat atc = atmat*cmat;

local awtilde awbar1 waw wacaw dm_YL dm_xs id_var awbari;

awtilde = 0;
awbar1 = 0;
mform(nrow=k,type=sym) waw = 0;
copy waw wacaw;
list(prefix=dm_) dm_YL YL;
list(prefix=dm_) dm_xs exog;
trend id_var;  ? 1 2 ... n*tm1
id_var = 1 + int((id_var-1)/tm1) ;  ?  1 1 1 ... 2 2 2 .... N N N ...
do i=1,n;
   select(silent) id_var=i;
   mat awtilde = ser(atc*eah);
   awbar1 = dm_YL - awtilde;
   mmake awbari awbar1 dm_xs;
   mat waw = waw + awbari'awbari;
   mat wacaw = wacaw + awbari'cmat*awbari;
enddo;

local dbar dbarinv qvec corr1-corr3 wacawd kcorr;

mat dbar=waw;
set dbar(1,1) = dbar(1,1) + sigeah*n*trcac;
mat dbarinv = yinv(dbar);
mform(nrow=k,ncol=1) qvec=0;  set qvec(1) = 1;
? The corr1-corr3 terms are slightly different from the ones
? in JOSIM1.GAU -- they include sigeah*dbarinv .
mat corr1 = sigeah*dbarinv* (n/tm1)*ici*(2*qvec - waw*dbarinv*qvec);
?
mat wacawd = wacaw*dbarinv;
?
mat corr2 = sigeah*dbarinv* ( tr(wacawd)*qvec + wacawd*qvec );
mat corr3 = sigeah*dbarinv* sigeah*n*(qvec'dbarinv*qvec) *
            ( (-n/tm1)*ici*trcac + 2*trcacac )*qvec;
print sigeah;
print corr1-corr3;
?
mat kcorr = corr1+corr2+corr3;
print kcorr;

mat blsdvc = coefw + kcorr;
print blsdvc coefw coefah;

? 4. VCOV for LSDVc , following Bun and Kiviet, 6/99
?    (This was not in Ruth Judson's 1996 code)
?   Note that T in Bun and Kiviet equals tm1 in Judson and Owen.
local QiWAW term2-term4 tm2 qvqvp Q1
;
mat QiWAW = n*vcovw/s2w;

mat term2 = sigeah*trcac*(QiWAW*qvec)*(QiWAW*qvec)';  ? PI_T'PI_T = cac

mat qvqvp = qvec*qvec';
set tm2 = tm1-1;
mat Q1 = sigeah*[  (1+gamma**tm2)/(1-gamma)**2 -
                 (2/tm1)*(1-gamma**tm1)/(1-gamma)**3]
               * QiWAW * qvqvp ;
mat term3 = Q1*(n*vcovah/s2ah)*Q1';

local I_tm2 tm20 K_T D_T trpdkc
;
mat I_tm2 = ident(tm2);
mform(nrow=tm2,ncol=1) tm20 = 0;
mmake K_T I_tm2 tm20;  ? lag operator
mmake D_T tm20 I_tm2;
mat D_T = D_T - K_T;   ? first difference operator
mat trpdkc = tr(atc*(D_T'K_T)*cmat);  ? PI_T = atc , L_T*GAMMA_T = cmat

local QWADZ QiZWstar d_YL Zi Wstari dm_Wi
;
mform(nrow=k,ncol=k,type=gen) QZWstar = 0;
copy QZWstar QWADZ;
list(prefix=d_) d_YL YL;
do i=1,n;
   select(silent) (id_var=i) & good_AHL;
   mmake Zi YL(-1) d_xs;
   mmake Wstari d_YL d_xs;
   mat QZWstar = QZWstar + Zi'Wstari;
   select(silent) (id_var=i);
   ?mmake Wi YL exog;
   ?mat QWADZ = QWADZ + Wi'atmat*D_T'Zi;
   mmake dm_Wi dm_YL dm_xs;
   mat QWADZ = QWADZ + dm_Wi'D_T'Zi;
enddo;
mat QZWstar = QZWstar/n;
mat QWADZ = QWADZ/n;

mat term4 = QiWAW*(QWADZ + sigeah*trpdkc*qvqvp)*((QZWstar")')*(Q1');

print QiWAW term2-term4;

mat vlsdvc = (sigeah/N)*( QiWAW + term2 + term3 + term4 );

title 'LSDVc estimates';
tstats(names=(YL,exog)) blsdvc vlsdvc;

smpl smpsav;
endproc;