> mdata <-read.table("ch2001.csv", header=T, sep=",")
> year<-mdata[,1]
> ch  <-mdata[,2]
> ydh <-mdata[,3]
> wh  <-mdata[,4]
> m1 <-lm(ch~ydh+wh)

> #library(foreign)
> library(lmtest)
> # Huber-White robust regression
> coeftest(m1, vcov = vcovHC(m1, "HC1")) 

t test of coefficients:

              Estimate Std. Error t value  Pr(>|t|)    
(Intercept) 6.0631e+04 3.3130e+04  1.8301  0.092174 .  
ydh         5.2795e-01 1.3717e-01  3.8490  0.002314 ** 
wh          4.3395e-02 6.8344e-03  6.3495 3.667e-05 ***
---
Signif. codes:  0 �e***�f 0.001 �e**�f 0.01 �e*�f 0.05 �e.�f 0.1 �e �f 1

> # Newey-West HAC estimator for heteroskedasticity and autocorrelations
> library(sandwich)
> coeftest(m1, vcov = NeweyWest(m1, lag=1, prewhite=FALSE, adjust=TRUE)) 

t test of coefficients:

              Estimate Std. Error t value  Pr(>|t|)    
(Intercept) 6.0631e+04 4.1557e+04  1.4590   0.17024    
ydh         5.2795e-01 1.7399e-01  3.0344   0.01038 *  
wh          4.3395e-02 7.3630e-03  5.8937 7.324e-05 ***
---
Signif. codes:  0 �e***�f 0.001 �e**�f 0.01 �e*�f 0.05 �e.�f 0.1 

> # Prais-Winsten for autocorrelations (not for heteroskedasticity)
> library(prais)
> sample <-data.frame(ch,ydh,wh)
> prais.winsten(ch~ydh+wh,data=sample)
[[1]]

Call:
lm(formula = fo)

Residuals:
    Min      1Q  Median      3Q     Max 
-2383.6 -1438.9  -976.6   483.6  6143.1 

Coefficients:
           Estimate Std. Error t value Pr(>|t|)    
Intercept 6.174e+04  5.098e+04   1.211 0.249193    
ydh       5.227e-01  2.122e-01   2.463 0.029883 *  
wh        4.365e-02  8.638e-03   5.054 0.000283 ***
---
Signif. codes:  0 �e***�f 0.001 �e**�f 0.01 �e*�f 0.05 �e.�f 0.1 �e �f 1

Residual standard error: 2636 on 12 degrees of freedom
Multiple R-squared:  0.9999,    Adjusted R-squared:  0.9999 
F-statistic: 5.836e+04 on 3 and 12 DF,  p-value: < 2.2e-16



******************************************
Prais-Winsten for autocorrelations (Stata)
******************************************
. prais ch ydh wh, rhotype(regress)

Iteration 0:  rho = 0.0000
Iteration 1:  rho = -0.0217
Iteration 2:  rho = -0.0245
Iteration 3:  rho = -0.0248
Iteration 4:  rho = -0.0248
Iteration 5:  rho = -0.0249
Iteration 6:  rho = -0.0249

Prais-Winsten AR(1) regression -- iterated estimates

      Source |       SS           df       MS      Number of obs   =        15
-------------+----------------------------------   F(2, 12)        =     87.77
       Model |  1.2201e+09         2   610056334   Prob > F        =    0.0000
    Residual |  83404542.5        12  6950378.54   R-squared       =    0.9360
-------------+----------------------------------   Adj R-squared   =    0.9254
       Total |  1.3035e+09        14  93108372.1   Root MSE        =    2636.4

------------------------------------------------------------------------------
          ch |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         ydh |   .5227037   .2122373     2.46   0.030     .0602784     .985129
          wh |   .0436551   .0086377     5.05   0.000     .0248353     .062475
       _cons |   61742.52    50983.2     1.21   0.249    -49340.33    172825.4
-------------+----------------------------------------------------------------
         rho |  -.0248518
------------------------------------------------------------------------------
Durbin-Watson statistic (original)    2.013586
Durbin-Watson statistic (transformed) 1.989451

.