Announcement

Collapse
No announcement yet.
X
  • Filter
  • Time
  • Show
Clear All
new posts

  • Linear regression with Newey and West standard errors

    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:
    Code:
    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
    Last edited by Christopher Koedding; 23 Dec 2016, 06:31.

  • #2
    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.

    Code:
    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

    Comment


    • #3
      Christophe, thank you very much indeed!

      When applying the code to my dataset, however, I get an error message. I did the following:
      Code:
      * TSSET the data
      tsset permno D
      Now I used your code to set up the regression in Mata:
      Code:
      (...)
      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.
      Code:
      * 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?

      Comment


      • #4
        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.


        Code:
        real rowvector myreg(real matrix Xall)
        {
        
        ....
        nw.m = st_numscalar("m")
        if (nw.m==.) nw.m = 0 
        ....
        
        
        }
        
        scalar m = 2
        rangestat ....

        Comment


        • #5
          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
          Last edited by Christopher Koedding; 24 Dec 2016, 03:31.

          Comment

          Working...
          X