Linear regression with Newey and West standard errors

[Christopher Koedding]

Hey guys,

I would like to set up a linear regression in Mata. Most importantly, I want to use heteroskedasticity-and-autocorrelation-consistent standard errors such as Newey and West standard errors instead of the regular ones.

With regular standard errors this is something like:

mata:
mata clear
mata set matastrict on
real rowvector myreg(real matrix Xall)
{
real colvector y, b, Xy
real matrix X, XX
y = Xall[.,1] // dependent var is first column of Xall
X = Xall[.,2::cols(Xall)] // the remaining cols are the independent variables
X = X,J(rows(X),1,1) // add a constant

XX = quadcross(X, X) // linear regression
Xy = quadcross(X, y)
b = invsym(XX) * Xy

return(rows(X), b')
}
end

Does anybody know how to set this up with Newey and West standard errors instead?

Help is highly appreciated.

Best,
Christopher

[Christophe Kolodziejczyk]

The formulas for the NW covariance matrix is given in the manual under the command newey. Here is one solution in Mata. I have used a structure, but it is not strictly necessary (and maybe an overkill). I just wanted to minimize the number of arguments I had to pass to the two functions I have created to ease the computations. The program should check that the number of lags is less than the number of observations. Note that I modified a bit your myreg() function.

use "http://www.stata-press.com/data/r12/idle2.dta" , clear

mata:
mata clear
mata set matastrict on


struct NW {
        
        real colvector y, e
        real matrix X, V
        real scalar m, N, k
        
        real rowvector b

}


real matrix computeLagCross(struct NW scalar nw, real scalar l)
{
        real matrix omega, x, xl
        real colvector e, el
        
        e  = nw.e[l+1..nw.N]
        el = nw.e[1..nw.N-l]
        
        x  = nw.X[l+1..nw.N,]
        xl = nw.X[1..nw.N-l,]
        
        omega = quadcross(x, 1, e:*el, xl, 1) + quadcross(xl, 1, e:*el, x, 1)
        
        return(omega)
}

real matrix newey(struct NW scalar nw)
{
        real scalar N, k, l, scaleFact
        real matrix omega, omega0
        
        scaleFact = nw.N/(nw.N-nw.k)
        omega = scaleFact*quadcross(nw.X, 1, nw.e:^2, nw.X, 1)

        if (nw.m>0)
        {
                for (l = 1 ; l <= nw.m ; l++)
                {
                        omega = omega + scaleFact*(1-l/(nw.m+1))*computeLagCross(nw,l)
                }
            
        }
        
        return(omega)


}

real rowvector myreg(real matrix Xall, real scalar m)
{
        real colvector y, b, Xy, e
        real matrix X, XX
        
        struct NW scalar nw
        
        nw.y = Xall[.,1] // dependent var is first column of Xall
        nw.X = Xall[.,2::cols(Xall)] // the remaining cols are the independent variables
        // X = X,J(rows(X),1,1) // add a constant

        XX = quadcross(nw.X, 1, nw.X, 1) // linear regression
        Xy = quadcross(nw.X, 1, nw.y, 0)
        nw.b = qrsolve(XX,Xy) // invsym(XX) * Xy
        
        nw.b
        nw.N = rows(nw.X)
        nw.k = cols(nw.X)+1
        nw.m = m
 
        nw.e = nw.y-(nw.X,J(rows(nw.X),1,1))*nw.b
        
        nw.V = invsym(XX)*newey(nw)*invsym(XX)
        st_matrix("NeweyWest",nw.V)
        
        return((nw.N, nw.b'))
}

X = st_data(.,"usr idle")
myreg(X,1)

end

tsset time
newey usr idle , lag(1)
mat li e(V)
mat li NeweyWest

[Christopher Koedding]

Christophe, thank you very much indeed!

When applying the code to my dataset, however, I get an error message. I did the following:

* TSSET the data
tsset permno D

Now I used your code to set up the regression in Mata:

(...)
X = st_data(.,"ret_rf CRSP_MKT_RF CRSP_MKT_RF_1")
myreg(X,1)

end

Then I want to execute rolling regressions in Stata using the -rangestat- command and referring to the regression set up in Mata.

* Use RANGESTAT to run rolling regressions of excess stock returns on excess market returns
* Regressions with a constant over a rolling window of 60 periods by permno
rangestat (myreg) ret_rf CRSP_MKT_RF CRSP_MKT_RF_1, by(permno) interval(D -59 0) casewise

This yields the error r(3001):

myreg(): 3001 expected 2 arguments but received 1
do_user_stats(): - function returned error
<istmt>: - function returned error


Do I have to adjust the Mata code to account for the fact that I want to use two instead of just one independent variable?
Or did I mess up something else here?

[Christophe Kolodziejczyk]

Apparently rangestat only accepts Mata function with only one argument (the matrix X of variables). So myreg has to have only one arguments and the number of lags for computing the NW covariance matrix has to be passed to the function in another way. One solution could be to use st_numscalar(). It's not ideal, but it should work, though I haven't tried.

real rowvector myreg(real matrix Xall)
{

....
nw.m = st_numscalar("m")
if (nw.m==.) nw.m = 0 
....


}

scalar m = 2
rangestat ....

[Christopher Koedding]

Christophe, thanks a lot! It really seems as if -rangestat- only works with one argument.

Would you be so kind to post the complete code (starting from mata: and ending with rangestat (myreg))? I'm still very new to Mata and do not want to make a mistake here.

Wishing you a merry Christmas!
Best,
Christopher