8
$\begingroup$

Could someone please help with the formula for the Newey West standard error of $\beta_1$ (without matrix notation) for the following regression:

$Y_t=\beta_0+\beta_1X_t+\epsilon_t$

where $\epsilon_t$ is believed to be serial correlated and/or heteroskedastic.

$\endgroup$
  • $\begingroup$ Can you be specific about your model and lags? $\endgroup$ – Dimitriy V. Masterov Aug 30 '16 at 18:53
  • $\begingroup$ Hi Dimitriy. I've edited the question to be more specific. Any help would be much appreciated $\endgroup$ – ts_highbury Aug 30 '16 at 19:00
  • $\begingroup$ Are you trying to implement newey west standard errors in a particular programming environment or just looking for the math? $\endgroup$ – Matthew Gunn Aug 30 '16 at 19:14
  • $\begingroup$ I'm just looking for the math! $\endgroup$ – ts_highbury Aug 30 '16 at 19:15
  • 1
    $\begingroup$ The end of stata.com/manuals13/tsnewey.pdf gives the formulas using matrix notation. I advise you learn matrix notation. $\endgroup$ – Matthew Gunn Aug 30 '16 at 19:37
7
$\begingroup$

Take your model of $$y_t=\beta_0+\beta_1x_t+u_t,$$ where $t=1,...,T$. We will assume there are no other regressors and that the serial correlation only lasts up to one period (so shocks do not persist for very long).

To get the Newey-West/HAC standard error of $\beta_1$ that is robust to heteroskedasticity and autocorrelation up to 1 lag, you should:

  1. Estimate the model with OLS, which gives you usual $SE(\beta_1)$, the RMSE $\hat \sigma$, and the residuals $\hat u_1,...,\hat u_T$.
  2. Get the ${\hat r_1,...,\hat r_T}$ residuals from the auxiliary OLS regression of $x_t$ on a constant and calculate $\hat a_t = \hat u_t \cdot \hat r_t$ for $t=1,...,T$. If you had more regressors, you would include them as additional covariates.
  3. Assuming that the serial correlation lasts up to one period, calculate $$\hat v(1)=\sum_{t=1}^T \hat a_t^2+ \sum_{t=2}^T \hat a_t \cdot \hat a_{t-1}.$$ The first term in $\hat v$ is what gets you the het-consistent standard error. The second term is the autocorrelation part. Assuming the autocorrelation is positive, this is why your standard errors blow up: you have less information. If you autocorrelation lasted longer, you would have additional weighted terms for each each lag. You might want to use a finite sample correction multiplier of $\frac{T}{T-k}$ here, where we count the constant as one of the $k$ covariates.
  4. Calculate N-W standard error of $\beta_1$ as $$\left[ \frac{SE(\beta_1)}{\hat \sigma} \right]^2 \cdot \sqrt{ \hat v(1)}.$$

Here's an example with Stata (with the full data shown by the list command in case you want to use other software). This has a slight wrinkle in that Stata uses a default finite sample correction of $\frac{T}{T-k}$ that is not default in other statistics packages and is usually not shown in textbook formulas, though it is a sensible thing to do:

. /* N-W Standard Errors With One Regressor and Serial Correlation That Dies Down After 1 Peri
> od */
. webuse idle2, clear

. tsset time
        time variable:  time, 1 to 30
                delta:  1 unit

. list time usr idle, clean noobs

    time   usr   idle  
       1     0    100  
       2     0    100  
       3     0     97  
       4     1     98  
       5     2     94  
       6     0     98  
       7     2     90  
       8     3     85  
       9     1     68  
      10     2     91  
      11     2     94  
      12     2     89  
      13     1     88  
      14     4     92  
      15     7     74  
      16     7     76  
      17     8     71  
      18     4     78  
      19     5     75  
      20    10     74  
      21    16     65  
      22    12     63  
      23     3     83  
      24     2     60  
      25     3     85  
      26     5     87  
      27     6     83  
      28     5     84  
      29     1     98  
      30     1     98  

. /* Step 1 */
. reg usr idle

      Source |       SS           df       MS      Number of obs   =        30
-------------+----------------------------------   F(1, 28)        =     28.07
       Model |   210.35435         1   210.35435   Prob > F        =    0.0000
    Residual |  209.812316        28  7.49329702   R-squared       =    0.5006
-------------+----------------------------------   Adj R-squared   =    0.4828
       Total |  420.166667        29  14.4885057   Root MSE        =    2.7374

------------------------------------------------------------------------------
         usr |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        idle |  -.2281501   .0430607    -5.30   0.000    -.3163559   -.1399442
       _cons |   23.13483    3.67706     6.29   0.000     15.60271    30.66694
------------------------------------------------------------------------------

. scalar se_beta1 = _se[idle]

. scalar sigmahat    = e(rmse)

. predict double uhat, resid

. /* Step 2 */
. reg idle

      Source |       SS           df       MS      Number of obs   =        30
-------------+----------------------------------   F(0, 29)        =      0.00
       Model |           0         0           .   Prob > F        =         .
    Residual |      4041.2        29  139.351724   R-squared       =    0.0000
-------------+----------------------------------   Adj R-squared   =    0.0000
       Total |      4041.2        29  139.351724   Root MSE        =    11.805

------------------------------------------------------------------------------
        idle |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _cons |       84.6    2.15524    39.25   0.000     80.19204    89.00796
------------------------------------------------------------------------------

. predict double rhat, resid

. gen double ahat = uhat*rhat

. /* Step 3 */
. gen double v = ahat^2

. replace v    = v + ahat*L1.ahat in 2/L
(29 real changes made)

. sum v

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
           v |         30    3070.853    8464.421  -256.7766   33901.46

. scalar v1     = r(sum) 

. scalar v1_fsc = r(sum)*(30/28) // Stata uses a finite sample correction of T/(T-k)

. /* Step 4 */
. di "Usual N-W SE = " sqrt(scalar(v1))*[scalar(se_beta1)/scalar(sigmahat)]^2
Usual N-W SE = .07510689

. di "Stata's N-W SE  = " sqrt(scalar(v1_fsc))*[scalar(se_beta1)/scalar(sigmahat)]^2
Stata's N-W SE  = .07774301

. /* Comapare to newey command */
. newey usr idle, lag(1)

Regression with Newey-West standard errors      Number of obs     =         30
maximum lag: 1                                  F(  1,        28) =       8.61
                                                Prob > F          =     0.0066

------------------------------------------------------------------------------
             |             Newey-West
         usr |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        idle |  -.2281501    .077743    -2.93   0.007    -.3873994   -.0689007
       _cons |   23.13483   7.119611     3.25   0.003     8.550965    37.71869
------------------------------------------------------------------------------

. di _se[idle]
.07774301

As you can see, with the finite sample correction, newey matches what we did by hand at 0.07774301. I am not sure this example contains a whole lot of intuition, but YMMV.

This is based on Wooldridge, Jeffrey M. "A computationally simple heteroskedasticity and serial correlation robust standard error for the linear regression model." Economics Letters 31.3 (1989): 239-243.


Stata Code:

/* N-W Standard Errors With One Regressor and Serial Correlation That Dies Down After 1 Period */
webuse idle2, clear
tsset time
list time usr idle, clean noobs
/* Step 1 */
reg usr idle
scalar se_beta1 = _se[idle]
scalar sigmahat = e(rmse)
predict double uhat, resid
/* Step 2 */
reg idle
predict double rhat, resid
gen double ahat = uhat*rhat
/* Step 3 */
gen double v = ahat^2
replace v    = v + ahat*L1.ahat in 2/L
sum v
scalar v1     = r(sum) 
scalar v1_fsc = r(sum)*(30/28) // Stata uses a finite sample correction of T/(T-k)
/* Step 4 */
di "Usual N-W SE = " sqrt(scalar(v1))*[scalar(se_beta1)/scalar(sigmahat)]^2
di "Stata's N-W SE  = " sqrt(scalar(v1_fsc))*[scalar(se_beta1)/scalar(sigmahat)]^2
/* Compare to newey command */
newey usr idle, lag(1)
di _se[idle]
/* Compare To Smaller Robust Sandwich SE */ 
sum v
scalar v0     = r(sum)*(30/28)
di "Stata's Sandwich SE  = " sqrt(scalar(v0))*[scalar(se_beta1)/scalar(sigmahat)]^2
reg usr idle, robust
di _se[idle]
| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Fantastic! Thank you Dimitriy. A quick question: I have seen papers say things like, "the lag length is 1.5 times the overlap". Any idea what the overlap is? $\endgroup$ – ts_highbury Sep 1 '16 at 4:58
  • $\begingroup$ @ts_highbury I think people often estimate TS models where the outcome is a difference $y_t-y_{t-1}$. For example, with monthly data, one observation might be the change from December to January, and another from January to February, and so on. This creates overlap, since January is a part of both observations, which induces autocorrelation even if there is none to begin with. Another example would be smoothed data, like a moving average. This heuristic tells you how to pick the maximum lag order of autocorrelation to deal with this issue. $\endgroup$ – Dimitriy V. Masterov Sep 1 '16 at 5:41
  • $\begingroup$ However, I am not entirely certain this is correct since this is not really my area. $\endgroup$ – Dimitriy V. Masterov Sep 1 '16 at 5:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.