*! version 0.2b, Antoni Sureda-Gomila, Monday, April 10, 2006
*       Out of alpha! Now it is still beta
* version 0.1, Antoni Sureda-Gomila, Saturday, April 8, 2006
* Copyright 2006: Antoni Sureda-Gomila
*-----------------------------------------------------------------------
program fmregress, eclass
    version 8
    syntax varlist(ts min=2)            ///
        [                               ///
            ,   NOCONStant              ///
                LAG(integer 0)          ///
                DROPFirst(integer 0)    ///
        ]
    *macro list
    *disp("I'm in")
    quietly capture xt_tis `t'
    if _rc {
        disp as err "Use tsset before fmreg"
        exit 111
    }
    *disp("Tsset checked")
    
    quietly summarize `s(timevar)'
    *disp("I've analyzed `s(timevar)'")
    local tmin = r(min) + `dropfirst'
    local tmax = r(max)
    local T = `tmax' - `tmin' + 1
    local tminplus1 = `tmin' + 1
    *disp("Time limits obtained tmin: `tmin' tmax: `tmax' T: `T'")
    * Regressions
    quietly regress `varlist' if `s(timevar)' == `tmin' ///
        , `noconstant'
    *disp("First regression done")
    matrix betas   = e(b)
    *disp("We have the betas of the first regression")
    matrix NumObs   = `e(N)'
    *disp("We have n of the firt regression")
    quietly forval i = `tminplus1'/`tmax' {
        regress `varlist' if `s(timevar)' == `i', `noconstant'
        matrix betas     = betas\e(b)
        matrix NumObsAux = `e(N)'
        matrix NumObs    = NumObs\NumObsAux
    }
    *disp("Regressions done")
    matrix U = J(rowsof(betas),1,1)
    *disp("Matrix U constructed")
    matrix beta = U'*betas/rowsof(betas)
    *disp("Coefficients estimated")
    matrix TotObsAux = U'*NumObs
    *disp("Total number of observations computed")
    matrix betasMinusBeta = betas - ///
        J(rowsof(betas),colsof(betas),1)*diag(beta)
    matrix V = betasMinusBeta'*betasMinusBeta / ///
            (rowsof(betasMinusBeta)*(rowsof(betasMinusBeta) - 1))
    *scalar d = det(V)
    *display d
    *disp("VCV matrix of coefficients estimated")
    forval j = 1/`lag' {
        *disp("I'm in the loop for Newey-West correction")
        scalar weight   = 1 - `j'/(1 + `lag')
        scalar lIni = 1 + `j'
        scalar lFin = rowsof(betasMinusBeta) - `j'
        matrix Vaux = ///
            betasMinusBeta[lIni...,1...]'*betasMinusBeta[1..lFin,1...] / ///
            ((rowsof(betasMinusBeta) - `j')*(rowsof(betasMinusBeta) - `j' - 1))
        matrix V = V + weight*(Vaux + Vaux')
    }
    *disp("I'm going to create ereturn")
    *scalar d = det(V)
    *display d
    ereturn post beta V
    *disp("Ereturn created")
    ereturn scalar N = TotObsAux[1,1]
    if "`noconstant'" != "" {
        ereturn scalar df_m = colsof(betas)
    }
    else {
        ereturn scalar df_m = colsof(betas) - 1
    }
    ereturn scalar df_r = e(N) - colsof(betas)
    ereturn scalar lag  = `lag'
    ereturn local cmd "fmregress"
    ereturn local depvar "`1'"
    *disp("ereturn filled")
    
    disp _n in gr ///
        `"Fama MacBeth regression with Newey-West standard errors"'
    disp in gr `"Number of obs  ="' in yel %10.0f e(N)
    disp in gr `"Maximum lag    ="' in ye %10.0f e(lag)   
    ereturn display
    
    matrix drop _all
    scalar drop _all
    *disp("I'm out")
end
*-----------------------------------------------------------------------