*! 1.0.0  11Aug97  Jeroen Weesie/ICS

program define hausman
        version 5.0

      * parse command line

        parse "`*'", p(",")
        if "`2'" ~= "," & "`2'" ~= "" { exit 198 }
        local mats "`1'"
        mac shift 2
        local opts "`*'"

        parse "`mats'", p(" ")
        local b1 "`1'"
        local V1 "`2'"
        local b2 "`3'"
        local V2 "`4'"

        if "`4'" == "" | "`5'" ~= "" {
                di in re "hausman expects 4 matrix arguments"
                exit 198
        }

        local options "Table Estnames(str) Zero noDisplay"
        parse ", `opts'"

      * matrix conformity

        if rowsof(`b1') ~= 1 | rowsof(`b2') ~= 1 {
                di in re "`b1' and `b2' should be row vectors"
                exit 503
        }
        if colsof(`b1') ~= colsof(`V1') | rowsof(`V1') ~= colsof(`V1') {
                di in re "b1 and V1 are not conformable"
                exit 503
        }
        if colsof(`b2') ~= colsof(`V2') | rowsof(`V2') ~= colsof(`V2') {
                di in re "b2 and V2 are not conformable"
                exit 503
        }

      * determine overlap in (eq:names) -of b1- and -b2-

        local names1 : colnames(`V1')
        local eq1    : coleq(`V1')

        local i1 1
        while `i1' <= colsof(`V1') {
                local name : word `i1' of `names1'
                if "`eq1'" ~= "" {
                        local eq : word `i1' of `eq1'
                        local name "`eq':`name'"
                }
                local i2 = colnumb(`V2',"`name'")
                * colnumb returns "." if -name- does not occur in V2

                if "`i2'" ~= "." & "`zero'" ~= "" {
                        if `V1'[1,`i1']==0 | `V2'[1,`i2']==0 {
                                local i2 "."
                        }
                }

                if "`i2'" ~= "." {              /* -name- occurs in both vectors */
                        local I1 "`I1'`i1' "
                        local I2 "`I2'`i2' "
                        local names "`names'`name' "
                }
                local i1 = `i1'+1
        }

        if "`I1'" == "" {
                di in re "hausman's test undefined: no overlap between (b1,V1) and (b2,V2)"
                global S_1 .
                global S_2 .
                global S_3
                exit 498
        }

      * create selections from b1, V1, b2, and V2

        tempname nb1 nV1 nb2 nV2 t

        local nI : word count `I1'
        mat `nb1' = J(1,`nI',0)
        mat `nV1' = J(`nI',`nI',0)
        mat `nb2' = J(1,`nI',0)
        mat `nV2' = J(`nI',`nI',0)

        local i 1
        while `i' <= `nI' {
                local ii : word `i' of `I1'
                mat `nb1'[1,`i'] = `b1'[1,`ii']
                local j 1
                while `j' <= `i' {
                        local jj: word `j' of `I1'
                        mat `nV1'[`i',`j'] = `V1'[`ii',`jj']
                        mat `nV1'[`j',`i'] = `V1'[`ii',`jj']
                        local j = `j'+1
                }
                local i = `i'+1
        }

        local i 1
        while `i' <= `nI' {
                local ii : word `i' of `I2'
                mat `nb2'[1,`i'] = `b2'[1,`ii']
                local j 1
                while `j' <= `i' {
                        local jj: word `j' of `I2'
                        mat `nV2'[`i',`j'] = `V2'[`ii',`jj']
                        mat `nV2'[`j',`i'] = `V2'[`ii',`jj']
                        local j = `j'+1
                }
                local i = `i'+1
        }

      * Table with common coefficients of both estimators

        if "`table'" ~= "" {
                if "`estnames'" ~= "" {
                        local e1 : word 1 of `estnames'
                        center 16 `e1'
                        local e1 "$S_1"
                        local e2 : word 2 of `estnames'
                        center 16 `e2'
                        local e2 "$S_1"
                }
                else {
                        local e1 "   Estimator 1  "
                        local e2 "   Estimator 2  "
                }

                di _n in gr "         |    `e1'      `e2'     Coef1-Coef2/ |"
                di    in gr "Variable |    Coef    Std Dev.      Coef    Std Dev.      sqrt(v2-v1) |"
                sepline
                local i 1
                local leq "_"
                while `i' <= `nI' {
                        local vn : word `i' of `names'
                        local ivn = index("`vn'",":")
                        local eq = substr("`vn'",1,`ivn'-1)
                        local vname = substr("`vn'",`ivn'+1,.)
                        if "`eq'" ~= "`leq'" {
                                if `i' ~= 1 { sepline }
                                di in ye "`eq'" in gr _col(10) "|" _col(71) "|"
                                local leq "`eq'"
                        }

                        #del ;
                        scalar `t' = (`nb1'[1,`i']-`nb2'[1,`i']) /
                                sqrt(`nV2'[`i',`i']-`nV1'[`i',`i']) ;
                        local c = 9-length("`vname'") ;
                        di in gr _col(`c') "`vname'" _col(10) "| "
                           in ye %9.0g `nb1'[1,`i'] " " %9.0g sqrt(`nV1'[`i',`i'])
                           "   " %9.0g `nb2'[1,`i'] " " %9.0g sqrt(`nV2'[`i',`i'])
                           _col(61)%7.2f `t' in gr _col(71) "|" ;
                        #del cr

                        local i = `i'+1
                }
                sepline
        }

      * compute H = (nb1-nb2)'(inv(nV2-nV1))*(nb1-nb2)

        tempname nb nV H r

        * invert V2-V1
        mat `nV' = `nV2' - `nV1'
        spdinv `nV'

        * rank of nV
        local i 1
        scalar `r' = 0
        while `i' <= `nI' {
                if `nV'[`i',`i'] ~= 0 { scalar `r' = `r'+1 }
                local i = `i'+1
        }
        if `r' < `nI' {
                di in bl "Warning: V2-V1 is semi-definite"
        }

        mat `nb' = `nb1' - `nb2'
        mat `H'  = `nV' * `nb''
        mat `H'  = `nb' * `H'

        global S_1 = `H'[1,1]
        global S_2 = `r'
        global S_3 "`names'"

      * display statistic

        if "`display'" == "" {
                di _n in gr _col(8) "Hausman's H" _col(24) "= " in ye %9.2f $S_1
                di in gr    _col(8) "Pr H > chi2($S_2) " _col(24) "= " /*
                    */ in ye %9.4f chiprob($S_2,$S_1)
        }
end

* inverse of symmetric (semi)positive-definite matrix
* note that Stata's syminv does not blow up on non-definite matrices
program define spdinv
        local V  "`1'"
        tempname EigVal EigVec
        capt matrix symeig `EigVec' `EigVal' = `V'
        local rc = _rc
        if `rc' > 0 {
                global S_1 .
                global S_2 .
                global S_3
                if `rc' == 505 {
                        di in re "The difference V2-V1 is not symmetric"
                        di in re "hausman was probably invoked with invalid VCEs"
                        exit 505
                }
                else    error `rc'
        }

      * if smallest eigenvalue is negative, V2-V2 is not spd
        if `EigVal'[1,colsof(`EigVal')] < 0 {
                global S_1 .
                global S_2 .
                global S_3
                di in re "The difference V2-V1 is not (semi-)positive definite"
                di in re "Hausman's test is not well-defined"
                exit 498
        }

      * form inverse the lazy way
        mat `V'= syminv(`V')
end

* "center m str" centers str in a field of width m
* the result is returned in S_1
program define center
        local m "`1'"
        mac shift
        local s = substr("`*'", 1, `m')
        local l = length("`s'")
        if `l' == `m' {
                global S_1 "`s'"
        }
        else {
                local b = int((`m'-`l')/2)
                local e = `m'-`l'-`b'
                global S_1 : display _dup(`b') " " "`s'" _dup(`e') " "
        }
end

program define sepline
        di in gr "---------+" _dup(60) "-" "+"
end

exit