*! version 0.2b, Antoni Sureda-Gomila, Monday, April 10, 2006
*       This is a beta
* Copyright 2006: Antoni Sureda-Gomila
*-----------------------------------------------------------------------
program fmivreg, eclass
    version 8
    syntax anything                     ///
        [                               ///
            ,   NOCONStant              ///
                LAG(integer 0)          ///
                DROPFirst(integer 0)    ///
        ]
    *macro list
    *disp("I'm in")
    *disp("`anything'")
    quietly capture xt_tis `t'
    if _rc {
        disp as err "Use tsset before fmreg"
        exit 111
    }
    *disp("Tsset checked")
    local time = "`s(timevar)'"
    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
    *disp("    ivreg `anything' if `time' == `tmin', `noconstant'")
    quietly ivreg `anything' if `time' == `tmin' ///
        , `noconstant'
    local instruments  = "`e(insts)'"
    local instrumented = "`e(instd)'"
    *disp("Instruments `e(insts)'")
    *disp("Instrumented `e(instd)'")
    *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' {
        disp("I'm in the regressions loop")
        disp("ivreg `anything' if `time' == `i', `noconstant'")
        ivreg `anything' if `time' == `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))
    *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')
    }
    
    ereturn post beta V
    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 "fmivreg"
    ereturn local depvar "`1'"
    ereturn local insts "`instruments'"
    ereturn local instd "`instrumented'"
    *disp("ereturn filled")
    
    disp _n in gr ///
        `"IV 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)
    disp in gr `"First time dropped ="' in ye %10.0f `dropfirst'
    ereturn display
    disp in gr "Instrumented:  " _c
    disp "`e(instd)'"
    disp in gr "Instruments:   " _c
    disp "`e(insts)'"
    di in smcl in gr "{hline 78}"
    
    matrix drop _all
    scalar drop _all
    *disp("I'm out")
end
*-----------------------------------------------------------------------