Hello all,
I'm attempting to run a SUR model and adjust for heteroskedasticity and
autocorrelation using var.hac, which implements the Newey-West
Heteroskedasticity and Autocorrelation Consistent covariance matrix. I'm
running Splus 7.0 with Finmetrics 2.0.1 on a Windows XP Pro.
I figured out how to estimate the variance-covariance matrix of a SUR using
matrix algebra, and I think all I need to do is swap my var() statement with
the var.hac() statement (see below). However, that's not how it works when
implementing an OLS with a Newey-West correction, so I'm a bit concerned.
Any advice would be hugely appreciated!! An example demonstrating my code
and my concerns is below.
Cheers,
Andy
# simulate a multivariate data set.
set.seed(10)
x1 <- abs(rnorm(100))
x2 <- x1^2
sim.dat <- data.frame(x1 = x1, x2 = x2, y1 = 1 + 2 * x1 + rnorm(100), y2 = 4
+ 5 * x2 + rnorm(100))
# estimate an SUR model
testsur = SUR(list(y1 ~ x1, y2 ~ x2), data = sim.dat, df = 3)
# Calculate the variance-covariance matrix by hand using Equation 10.10 in
Zivot and Wang (2006)
# re-arrange the data
xx1 = cbind(1, x1, 0, 0)
xx2 = cbind(0, 0, 1, x2)
xx = rbind(xx1, xx2)
# calculate sigma using var() and correcting for the default degrees of
freedom (eg, Equation 10.9 in Z & W) - this is where I think I'd swap in
var.hac(), but let var.hac() do the degrees of freedom adjustment
my.sigma = (var(residuals(testsur)) * 99)/100
# estimate the assymptotic variance-covariance
my.v = kronecker(solve(my.sigma), diag(100))
avarB = solve(crossprod(xx, my.v) %*% xx)
# Compare 10.10 with SUR output - not exact, but close
avarB
x2
0.0265833644 -0.021461723 0.001960952 -0.0008997789
-0.0214617234 0.026941248 -0.001113889 0.0011295070
0.0019609523 -0.001113889 0.015912270 -0.0049318926
x2 -0.0008997789 0.001129507 -0.004931893 0.0050010425
vcov(testsur)
y1.(Intercept) y1.x1 y2.(Intercept) y2.x2
y1.(Intercept) 0.0265833667 -0.021461730 0.001960581 -0.0008996088
y1.x1 -0.0214617297 0.026941256 -0.001113679 0.0011292934
y2.(Intercept) 0.0019605807 -0.001113679 0.015912246 -0.0049318869
y2.x2 -0.0008996088 0.001129293 -0.004931887 0.0050010368
# But simply swapping var() and var.hac() is not how it works in OLS
LS.fit = OLS(y1 ~ x1, data = sim.dat)
my.x1 = cbind(1, x1)
# Calculate the variance-covariance matrix by hand based on equation 6.3
from Zivot and Wang (2006) and adjust for the default degrees of freedom in
var()
xpx = solve(crossprod(my.x1))
ols.cov = (var(residuals(LS.fit)) * xpx * 99)/98
# Compare 6.3 with the OLS output with no correction
ols.cov
x1
0.02715283 -0.02193355
x1 -0.02193355 0.02753354
vcov(LS.fit)
(Intercept) x1
(Intercept) 0.02715283 -0.02193355
x1 -0.02193355 0.02753354
# But to get the robust variance-covariance matrix we use equation 6.20 in
Z&W (as also demonstrated in the var.hac() help)
xpx.rob = solve(crossprod(my.x1)/nrow(my.x1))
hh = my.x1 * LS.fit$residuals
ols.cov.rob = xpx.rob %*% var.hac(hh, window = "bartlett", df = 2) %*%
xpx.rob/100
# Compare var.hac approach with OLS with Newey-West Correction
> sqrt(diag(ols.cov.rob))
[1] 0.1578415 0.1485211
> summary(LS.fit, correction = "nw")
....
Coefficients:
Value Std. Error t value Pr(>|t|)
(Intercept) 0.9827 0.1578 6.2259 0.0000
x1 1.9583 0.1485 13.1855 0.0000
# Which is NOT the same as just inserting var.hac() into equation 6.3
ols.cov.silly = var.hac(residuals(LS.fit), window = "bartlett", df = 2) *
xpx
> sqrt(diag(ols.cov.silly))
[1] 0.1428708 0.1438690
> sqrt(diag(ols.cov.rob))
[1] 0.1578415 0.1485211
*****************************************************
Andrew B. Cooper
Associate Professor
School of Resource and Environmental Management
Simon Fraser University
TASC 1 Building, Room 8405
8888 University
Burnaby, British Columbia, CANADA V5A 1S6
*****************************************************
Office: 778-782-3954
FAX: 778-782-4968
----------------------------------------------------------
Web Page: http://www.rem.sfu.ca/faculty/cooper.htm
Lab Web Page: http://www.rem.sfu.ca/fishgrp/
----------------------------------------------------------
|