/*
**        PROGRAM FOR HT ESTIMATION OF PANEL DATA
**
**                     WRITTEN BY
**                   SEUNG CHAN AHN
**              DEPARTMENT OF ECONOMICS
**                COLLEGE OF BUSINESS
**                  TEMPE, AZ 85287
**
*/

/*
**  Computing HT, AM and BMS Estimators
**  With Hausman Tests
*/

	new ;

@ Put MGIV.COL in the directory you execute this program! @

	#include mgiv.col ;

@ Open Output file @

	output file = pan_ht.out reset ;

@ Format output file @

	format /rd 12,4 ;


@   Provide # of ovservations and # of variables @

    nobs = 4165 ;
    nvar = 23   ;

@   Read Data   @

    load dat[nobs,nvar] = cr.db ;
@    dat = delif(dat,dat[.,10] .== 1) ; @

@  Define Variables @

    id68  = dat[.,1]  ;
    expp  = dat[.,2]  ;
    expp2 = dat[.,3]  ;
    wks   = dat[.,4]  ;
    occ   = dat[.,5]  ;
    ind   = dat[.,6]  ;
    south = dat[.,7]  ;
    smsa  = dat[.,8]  ;
    ms    = dat[.,9]  ;
    fem   = dat[.,10] ;
    unionm= dat[.,11] ;
    edu   = dat[.,12] ;
    blk   = dat[.,13] ;
    wage  = dat[.,14] ;
    y76      = dat[.,17] ;
    y77      = dat[.,18] ;
    y78      = dat[.,19] ;
    y79      = dat[.,20] ;
    y80      = dat[.,21] ;
    y81      = dat[.,22] ;
    y82      = dat[.,23] ;
    lwage    = ln(wage)  ;
    lwks     = ln(wks)   ;
    dyr      = y77~y78~y79~y80~y81~y82 ;

@   Define N and T @

    n  = 4165/7 ;
@    n  = 528    ; @
    t  = 7      ;

@   Define Dep. Var., Time-varying Reg. and Time-invariant Reg. @

    @ Treating expp as time-invariant @
    rexp = {0,1,2,3,4,5,6}    ;
    rexp = ones(n,1) .*. rexp ;
    rexp = expp - rexp        ;

    yy  = lwage  ;
    xx  = wks~south~smsa~ms~(expp^2)~occ~ind~unionm~dyr ;
    zz  = ones(n*t,1)~fem~blk~edu~rexp ;
    xx1 = wks~south~smsa~ms   ;
    zz1 = ones(n*t,1)~fem~blk  ;

	vny   = {"logwage"};
	vnx   = {"wks", "south", "smsa", "ms", "expp2", "occ", "ind", "unionm",
			 "yr77", "yr78", "yr79", "yr80", "yr81", "yr82"};
    vnz   = {"cons","fem","blk","edu","exp"};

/*
**  From Here, Do Not Change
*/

    clear dat   ;

@   Define k and g @

    k1 = cols(xx1) ;
    g1 = cols(zz1) ;
    k  = cols(xx)  ;
    g  = cols(zz)  ;

@ Creating AM, BMS and Mean Variables @

    qvxx  = qvmat(xx,n,t)  ;
    pvxx  = pvmat(xx,n,t)  ;
    qvyy  = qvmat(yy,n,t)  ;
    pvyy  = pvmat(yy,n,t)  ;
    pvxx1 = pvmat(xx1,n,t) ;
    amxx1 = ammat(xx1,n,t) ;
    bmsxx1= bmsmat(xx1,n,t) ;

@   Within Estimation  @

    wb = invpd(qvxx'qvxx)*(qvxx'qvyy) ;

    we   = qvyy - qvxx*wb       ;
    ssq  = (we'we)/(n*t-n-k)    ;
    wc   = ssq*invpd(qvxx'qvxx) ;
    ws   = sqrt(diag(wc))       ;

	let mask[1,4] = 0 1 1 1;
	let fmt[4,3] =
    "-*.*s" 8 8
    "*.*lf" 10 4
    "*.*lf" 10 4
    "*.*lf" 10 4;


	"Within Estimation Result" ;
	"------------------------" ;
	" dependent variable: " $vny ;
	"";
	"variable     coeff.  std. err.   t-st " ;
	yyprin = printfm(vnx~wb~ws~(wb./ws),mask,fmt); 
	"" ;

@   HT Estimation  @

        iv   = pvxx1~zz1 ;
        wee  = pvyy - pvxx*wb ;

        wg   = invpd( (zz'iv)*invpd( iv'iv )*(iv'zz) )
               *(zz'iv)*invpd( iv'iv )*(iv'wee)        ;

        auxe  = wee - zz*wg         ;
        ssqbb = (auxe'auxe)/(n-g)   ;
        theta = sqrt( ssq/ssqbb )   ;
		ssqaa = (ssqbb-ssq)/t       ;

        yys = qvyy + theta*pvyy ;
        xxs = qvxx + theta*pvxx ;
        zzs = theta*zz          ;

        htiv = qvxx~pvxx1~zz1 ;
        regs = xxs~zzs      ;

        htd   = invpd( (regs'htiv)*invpd( htiv'htiv )*(htiv'regs) )
                *(regs'htiv)*invpd( htiv'htiv )*(htiv'yys)       ;
        htc   = ssq* invpd( (regs'htiv)*invpd( htiv'htiv )*(htiv'regs) ) ;
        hts   = sqrt(diag(htc)) ;


	"Hausman_Taylor Estimation Result" ;
	"------------------------" ;
	" dependent variable: " $vny ;
	"";
	"variable     coeff.  std. err.   t-st " ;
	yyprin = printfm((vnx|vnz)~htd~hts~(htd./hts),mask,fmt); 
	"" ;

	"     THETA = " theta ;
	"     SIGE2 = " ssq   ;
	"	  SIGAA = " ssqaa ;

    @ H TEST FOR HT VS. W @

    gb   = htd[1:k]        ;
    gbc  = htc[1:k,1:k]    ;
    ht   = (wb-gb)'pinv(wc-gbc)*(wb-gb)  ;
    df   = k1 - g + g1     ;

    "     HAUSMAN TEST FOR HT VS. WITHIN  " ;
    "     STATISTIC, P-VAL, DF =" ht cdfchic(ht,df) df ;
    "" ;

    @ HANSEN TEST FOR HT @

    he   = yys - regs*htd    ;
    hthant = n*t*(he'htiv)*invpd(htiv'htiv)*(htiv'he)/(he'he) ;
    htdf   = cols(htiv) - k - g  ;

    "     HANSEN TEST FOR HT " ;
    "     STATISTIC, P-VAL, DF =" hthant cdfchic(hthant,htdf) htdf ;
    "" ;

@   AM Estimation  @

        amiv  = qvxx~zz1~amxx1 ;
        amd   = invpd( (regs'amiv)*invpd( amiv'amiv )*(amiv'regs) )
                *(regs'amiv)*invpd( amiv'amiv )*(amiv'yys)              ;
        amc   = ssq*invpd( (regs'amiv)*invpd( amiv'amiv )*(amiv'regs) ) ;
        ams   = sqrt( diag(amc) )                               ;


	"AM Estimation Result" ;
	"------------------------" ;
	" dependent variable: " $vny ;
	"";
	"variable     coeff.  std. err.   t-st " ;
	yyprin = printfm((vnx|vnz)~amd~ams~(amd./ams),mask,fmt); 
	"" ;

    @ H TEST FOR AM VS. W @

    gb   = amd[1:k]        ;
    gbc  = amc[1:k,1:k]    ;
    ht   = (wb-gb)'pinv(wc-gbc)*(wb-gb)  ;
    nnn  = rank(wc-gbc)|(cols(amiv)-k-g) ; 
    df   = minc(nnn)                     ;

    "     HAUSMAN TEST FOR AM VS. WITHIN " ;
    "     Statistic, p-val, df =" ht cdfchic(ht,df) df ;
    "" ;

    @ H TEST FOR AM VS. HT @

    amdf = cols(amiv) - k - g ;
    ht   = (htd-amd)'pinv(htc-amc)*(htd-amd)    ;
    nnn  = rank(htc-amc)|(amdf-htdf)            ; 
    df   = minc(nnn)                            ;

    "     HAUSMAN TEST FOR AM VS. HT " ;
    "     Statistic, p-val, df =" ht cdfchic(ht,df) df ;
    "" ;

    @ HANSEN TEST FOR AM @

    he   = yys - regs*amd    ;
    amhant = n*t*(he'amiv)*invpd(amiv'amiv)*(amiv'he)/(he'he) ;

    "     HANSEN TEST FOR AM " ;
    "     STATISTIC, P-VAL, DF =" amhant cdfchic(amhant,amdf) amdf ;
    "" ;

    @ EHS TEST FOR AM VS. HT @

    amehst = amhant - hthant ;
    df     = amdf - htdf     ;

    "     EHS TEST FOR AM VS. HT " ;
    "     STATISTIC, P-VAL, DF =" amehst cdfchic(amehst,df) df ;
    "" ;

OUTPUT OFF