* HETEROSKEDASTIC CONSISTENT COVARIANCE MATRICES
sample 1 17
read year consume income price
1923  99.2  96.7 101.0
1924  99.0  98.1 100.1
1925 100.0 100.0 100.0
1926 111.6 104.9  90.6
1927 122.2 104.9  86.5
1928 117.6 109.5  89.7
1929 121.1 110.8  90.6
1930 136.0 112.3  82.8
1931 154.2 109.3  70.1
1932 153.6 105.3  65.4
1933 158.5 101.7  61.3
1934 140.6  95.4  62.5
1935 136.2  96.4  63.6
1936 168.0  97.6  52.6
1937 154.3 102.4  59.7 
1938 149.0 101.6  59.5
1939 165.5 103.8  61.3
* Run ols, save residuals and standard errors
ols consume income price / resid=u stderr=ose 
gen1 n=$n
gen1 df=$df
* Square the residuals 
genr u2=u**2
* Copy the independent variables into the X matrix
genr one=1
copy income price one x
matrix xx=inv(x'x)
* Create White's (1980) covariance matrix adjusted for heteroskedasticity.
* (automatically calculated with the HETCOV option on the OLS command).
matrix s0=x'diag(u2)*x/n
matrix hc=n*xx*s0*xx
* Get the corrected standard errors 
matrix hse=sqrt(diag(hc))

* HC1 - The Hinkley method of estimation
matrix hc1=(n/df)*hc
matrix hse1=sqrt(diag(hc1))
* HC2 - The Horn and Duncan estimate
matrix ktt=diag(x*xx*x')
matrix sig2=u2/(1-ktt)
matrix hc2=xx*x'diag(sig2)*x*xx
matrix hse2=sqrt(diag(hc2))
* HC3 - The MacKinnon and White (1985) Jackknife estimate.
matrix ustar=u/(1-ktt)
matrix om=diag(ustar**2)
matrix hc3=((n-1)/n)*xx*(x'om*x-(1/n)*(x'ustar*ustar'x))*xx
matrix hse3=sqrt(diag(hc3))

* Print the different standard error estimates
sample 1 3
namefmt(3x,6(2x,a8)
format(6f10.3)
print ose hse hse1 hse2 hse3 / format
* Print the 4 covariance matrix estimates
print hc hc1 hc2 hc3

* The Cragg (1983) estimator using x**2 as auxiliary variables
sample 1 17
genr income2=income**2
genr price2=price**2
copy income price one income2 price2 q
matrix qq=q*inv(q'diag(u2)*q)*q'
* The coefficient vector ba is Cragg's formula (13), p. 753
matrix ba=inv(x'qq*x)*x'qq*consume
* The covariance matrix is Cragg's formula (14), p. 754 
matrix vba=inv(x'qq*x)
print ba vba

* The Newey-West (1987) variance estimator for autocorrelated errors.
* Reference: Greene (2000), pp. 537-538.
gen1 L=2
ols consume income price / autcov=L pcov
set nodoecho
do #=1,L
gen1 beg=#+1
gen1 w=1-(#/(L+1))
do %=beg,n
matrix omega=u(%)*u(%-#)*x(%,0)'*x(%-#,0)/n
matrix s0=s0 + w*(omega+omega')
endo
endo
matrix autcov=n*xx*s0*xx
print autcov
STOP