/*
**         PROGRAM FOR WITHIN and GLS ESTIMATION OF PANEL DATA
**
**                          WRITTEN BY
**                        SEUNG CHAN AHN
**                   DEPARTMENT OF ECONOMICS
**                     W.P CAREY SCHOOL OF BUSINESS
**                       TEMPE, AZ 85287
**
*/

/*
**  Computing within and GLS Estimators
**  With Hausman Tests
*/

new  ;

@  You must locate MGIV.COL in the directory you execute this program @
 
	#include mgiv.col ;

@	Open an output file  @

	output file = pan_gls.out reset ;

@	Formatting output file  @

	format /rd 12,4 ;

@   Provide # of observations and # of variables @

    nobs = 336 ;
    nvar = 13   ;

@   Read Data   @

    load dat[nobs,nvar] = auto_1.txt ;
	@ 48 states (N = 48), 1982 - 1988 (T = 7)@


@  Define Variables @

    id       = dat[.,1]  ; @ ID for States @
    year     = dat[.,2]  ; @ year @
    spircons = dat[.,3]  ; @ Spirits consumption @
    unrate   = dat[.,4]  ; @ Unemployment rate @
    perinc   = dat[.,5]  ; @ Personal Income @
    emppop   = dat[.,6]  ; @ Employment/Population Ration @
    beertax  = dat[.,7]  ; @ Tax on Case of Beer @
    mlda     = dat[.,8]  ; @ Minimum Legal Drinking Age @
	vmiles   = dat[.,9]  ; @ Ave. mile per driver @
    jaild    = dat[.,10] ; @ Mandatory Jail Sentense = 1 @
    comserd  = dat[.,11] ; @ Mandotory Jail Sentence @
	allmort  = dat[.,12] ; @ # of Vehicle Fatalities @
	mrall    = dat[.,13] ; @ Vehicle Fatality Rate (VFR) @

@ Creating Time dummy variables @

	v    = {1982.5, 1983.5, 1984.5, 1985.5, 1986.5, 1987.5 };
    dyr  = dummy(year,v);

@   Define N and T @

    t  = 7   ;
    n  = rows(dat)/t ;

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

    yy    = mrall*10000  ; @ dependent var. @
    xx    = beertax~mlda~jaild~comserd~unrate~ln(perinc)~dyr[.,2:7]; @ time-varying indep. @
	zz    = ones(rows(yy),1); @ time-invariant indep. @ 

	vny   = {"VFR"};
	vnx   = {"beertax","mlda","jailed","comserd","unrate","lpinc",
			 "yr83", "yr84", "yr85", "yr86", "yr87", "yr88"};
    vnz   = {"cons"};

@ Exclude year dummy vars. from xx to make pvxx full column @
@ Use xxt for ALT1 test @

    xxt   = beertax~mlda~jaild~comserd~unrate~ln(perinc) ;

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

    clear dat ;

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

@   Define k and g @

    k  = cols(xx)  ;
    kt = cols(xxt) ;
    g  = cols(zz)  ;

@ Creating AM and Mean Variables @

    pvxxt = pvmat1(xxt,n,t)   ;
    pvxx  = pvmat1(xx,n,t)    ;
    pvyy  = pvmat1(yy,n,t)    ;

@   creating Deviation-From-Mean Variables  @

    qvxx = qvmat(xx,n,t) ;
    qvyy = qvmat(yy,n,t) ;

@	OLS estimation  @

	ow   = xx~zz;
	od   = invpd(ow'ow)*(ow'yy);
	oe   = yy - ow*od;
	os2  = (oe'oe)/(rows(ow)-cols(ow));
	ov   = os2*invpd(ow'ow);
	ose  = sqrt(diag(ov));
	orsq = 1-(oe'oe)/(yy'yy-rows(yy)*meanc(yy)^2);

"OLS Estimation Result" ;
"------------------------" ;
" dependent variable: " $vny ;
"";
"variable     coeff.  std. err.   t-st " ;
yyprin = printfm((vnx|vnz)~od~ose~(od./ose),mask,fmt); 
"" ;
"R-Square =" orsq;
"";

@   Within Estimation  @

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

    we   = qvyy - qvxx*wb       ;
    ssq  = (we'we)/(n*t-n-k)   ;
    wc   = ssq*inv(qvxx'qvxx)  ;
    ws   = sqrt(diag(wc)) ;
	wrsq = 1 - (we'we)/(yy'yy-rows(yy)*meanc(yy)^2);

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

@   GLS estimation  @

    bx    = xx~zz     ;
    bd    = invpd(bx'bx)*(bx'yy)     ;
    bee   = pvyy - pvmat1(bx,n,t)*bd ;
    ssqbb = (bee'bee)/(n-k-g) ;
    theta = sqrt(ssq/ssqbb)   ;
    ssqaa = (ssqbb-ssq)/t ;

    yystar = yy - (1-theta)*pvyy ;
    xxstar = xx - (1-theta)*pvxx ;
    zzstar = theta*zz ;

    regstar = xxstar~zzstar ;

    gd  = invpd(regstar'regstar)*(regstar'yystar) ;
    gc  = ssq*invpd(regstar'regstar)              ;
    gs  = sqrt(diag(gc)) ;
	gee  = qvyy -qvxx*gd[1:cols(xx)] ;
	grsq = 1 - gee'gee/(yy'yy-rows(yy)*meanc(yy)^2);

"GLS Estimation Result" ;
"------------------------" ;
" dependent variable: " $vny ;
"";
"variable     coeff.  std. err.   t-st " ;
yyprin = printfm((vnx|vnz)~gd~gs~(gd./gs),mask,fmt); 
"" ;
"R-Square =" grsq;
"";

    " THETA = " theta   ;
    " SIGE2 = " ssq     ;
    " SIGA2 = " ssqaa   ;
    " S.E.R.= " sqrt(ssq) ;
    "" ;

@ CREATING GLS RESIDUALS @

    gee  = yystar - regstar*gd  ;

@ H TEST FOR W VS. GLS @

    gb   = gd[1:k]        ;
    gbc  = gc[1:k,1:k]    ;
    ht   = (wb-gb)'pinv(wc-gbc)*(wb-gb)  ;
    df   = rank(wc-gbc)   ;

    "Hausman Test, p-val, df =" ht cdfchic(ht,df) df ;

@ J TEST FOR W VS. GLS USING PX ONLY @

    axx = qvxx~pvxxt~zz ;
    abb = invpd(axx'axx)*(axx'gee) ;
    ru2 = (axx*abb)'(axx*abb)/(gee'gee) ;
    alt1 = t*n*ru2     ;
    df   = cols(pvxxt) ;

    "ALT1    Test, p-val, df =" alt1 cdfchic(alt1,df) df ;

OUTPUT OFF