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

/*
**  COMPUTING MGIV and GMM FOR HT, AM AND BMS ESTIMATOR
**  COMPUTING GMM USING HT and AM INSTRUMENTS
**  HAUSMAN TESTS
**  HANSEN TESTS
**	EHS TESTS
**
*/

	new ;

@	Locate MGIV.COL in the directory you execute this program @

	#include mgiv.col ;

@	Open output file @

	output file = pan_mgiv.out reset ;

@	Format output file @

	format /rd 12,4 ;

@   Provide # of ovservations and # of variables @

    nobs = 4165 ;
	tim  = 7    ;
    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] ;
    UNKNOWN  = dat[.,15] ;
    UNKNOWN2 = dat[.,16] ;
    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  = rows(dat)/tim ; 
    t  = tim      ;

@   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 Instrumental Variables  @

    htiv1 = qvmat(xx,n,t)            ;
    htiv2 = pvmat(xx1,n,t)~zz1       ;
    amiv2 = ammat(xx1,n,t)~zz1       ;

@   Kiefer's Estimation @

    {kb,kcov,kcovh} = kiefer(xx,yy,n,t) ;
    ks  = sqrt(diag(kcov))  ;
    ksh = sqrt(diag(kcovh)) ;


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

	"Kiefer's Within Estimation Results" ;
	"------------------------" ;
	" dependent variable: " $vny ;
	"";
	"variable     coeff.  std. err.   t-st " ;
	yyprin = printfm(vnx~kb~ks~(kb./ks),mask,fmt); 
	"" ;

	"Kiefer's Within Estimation Results (HETERO ADJUSTED)" ;
	"------------------------" ;
	" dependent variable: " $vny ;
	"";
	"variable     coeff.  std. err.   t-st " ;
	yyprin = printfm(vnx~kb~ksh~(kb./ksh),mask,fmt); 
	"" ;

/*
**  HT estimation and Specification Test
*/

@   HT MGIV Estimation  @

        {htd,htcov,htcovh} = mgiv2(xx~zz,htiv1,htiv2,yy,n,t) ;
        hts  = sqrt(diag(htcov))  ;
        htsh = sqrt(diag(htcovh)) ;

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

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

@   Hausman Test for HT Model  @

    {htb,htbcov,htbcovh} = mgiv1(xx,zz,htiv1,htiv2,yy,n,t) ;
    q  = htb - htd[1:k]          ;
    v  = htbcov - htcov[1:k,1:k] ;
    ht = q'pinv(v)*q    ;
    df = minc(rank(v)|(k1-g+g1)) ;

    "Hausman Test for HT MGIV VS. Within MGIV" ;
	"     stat, p-val, df =" ht cdfchic(ht,df) df ;
    "" ;

@   HT GMM Estimation @

        {htd,htcov,han} = gmm(xx~zz,htiv1~htiv2,yy,n,t) ;
        hts  = sqrt(diag(htcov))  ;
        hthant = han ;
        htdf   = cols(htiv1~htiv2) - k - g ;

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

@   Hansen Test for HT Model  @

	"Hansen Test for HT" ;
    "     stat, p-val, df =" hthant cdfchic(hthant,htdf) htdf ;
    "" ;

/*
**  AM estimation and Specification Test
*/

@   AM MGIV Estimation  @

        {amd,amcov,amcovh} = mgiv2(xx~zz,htiv1,amiv2,yy,n,t) ;
        ams  = sqrt(diag(amcov))  ;
        amsh = sqrt(diag(amcovh)) ;

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

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

@   Hausman Test for AM Model  @

    {amb,ambcov,ambcovh} = mgiv1(xx,zz,htiv1,amiv2,yy,n,t) ;
    q  = amb - amd[1:k]          ;
    v  = ambcov - amcov[1:k,1:k] ;
    ht = q'pinv(v)*q    ;
    df = minc( rank(v)|(cols(htiv1~amiv2)-k-g) ) ;

    "Hausman Test for AM MGIV VS. Within MGIV" ;
	"     stat, p-val, df =" ht cdfchic(ht,df) df ;
    "" ;

@   AM GMM Estimation @

        {amd,amcov,han} = gmm(xx~zz,htiv1~amiv2,yy,n,t) ;
        ams  = sqrt(diag(amcov))  ;
        amhant = han ;
        amdf   = cols(htiv1~amiv2) - k - g  ;

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

@   Hansen Test for AM Model  @

	"Hansen Test for AM";
    "     stat, p-val, df =" amhant cdfchic(amhant,amdf) amdf ;
    "" ;

@   EHS Test for AM VS. HT  @

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

	"EHS Test" ;
    "     stat, p-val, df =" ehst cdfchic(ehst,df) df ;
    "" ;

output off