1### R code from vignette source 'sandwich.Rnw' 2 3################################################### 4### code chunk number 1: preliminaries 5################################################### 6library("zoo") 7library("sandwich") 8library("strucchange") 9library("lmtest") 10options(prompt = "R> ", continue = "+ ") 11 12 13################################################### 14### code chunk number 2: hac-kweights 15################################################### 16curve(kweights(x, kernel = "Quadratic", normalize = TRUE), 17 from = 0, to = 3.2, xlab = "x", ylab = "K(x)") 18curve(kweights(x, kernel = "Bartlett", normalize = TRUE), 19 from = 0, to = 3.2, col = 2, add = TRUE) 20curve(kweights(x, kernel = "Parzen", normalize = TRUE), 21 from = 0, to = 3.2, col = 3, add = TRUE) 22curve(kweights(x, kernel = "Tukey", normalize = TRUE), 23 from = 0, to = 3.2, col = 4, add = TRUE) 24lines(c(0, 0.5), c(1, 1), col = 6) 25lines(c(0.5, 0.5), c(1, 0), lty = 3, col = 6) 26lines(c(0.5, 3.2), c(0, 0), col = 6) 27curve(kweights(x, kernel = "Quadratic", normalize = TRUE), 28 from = 0, to = 3.2, col = 1, add = TRUE) 29 30text(0.5, 0.98, "Truncated", pos = 4) 31text(0.8, kweights(0.8, "Bartlett", normalize = TRUE), "Bartlett", pos = 4) 32text(1.35, kweights(1.4, "Quadratic", normalize = TRUE), "Quadratic Spectral", pos = 2) 33text(1.15, 0.29, "Parzen", pos = 4) 34arrows(1.17, 0.29, 1, kweights(1, "Parzen", normalize = TRUE), length = 0.1) 35text(1.3, 0.2, "Tukey-Hanning", pos = 4) 36arrows(1.32, 0.2, 1.1, kweights(1.1, "Tukey", normalize = TRUE), length = 0.1) 37 38 39################################################### 40### code chunk number 3: loadlibs1 41################################################### 42library("sandwich") 43library("lmtest") 44 45 46################################################### 47### code chunk number 4: hc-data 48################################################### 49data("PublicSchools") 50ps <- na.omit(PublicSchools) 51ps$Income <- ps$Income * 0.0001 52 53 54################################################### 55### code chunk number 5: hc-model 56################################################### 57fm.ps <- lm(Expenditure ~ Income + I(Income^2), data = ps) 58 59 60################################################### 61### code chunk number 6: hc-test1 62################################################### 63coeftest(fm.ps, df = Inf, vcov = vcovHC(fm.ps, type = "HC0")) 64 65 66################################################### 67### code chunk number 7: hc-test2 68################################################### 69coeftest(fm.ps, df = Inf, vcov = vcovHC(fm.ps, type = "HC4")) 70 71 72################################################### 73### code chunk number 8: hc-plot 74################################################### 75plot(Expenditure ~ Income, data = ps, 76 xlab = "per capita income", 77 ylab = "per capita spending on public schools") 78inc <- seq(0.5, 1.2, by = 0.001) 79lines(inc, predict(fm.ps, data.frame(Income = inc)), col = 4, lty = 2) 80fm.ps2 <- lm(Expenditure ~ Income, data = ps) 81abline(fm.ps2, col = 4) 82text(ps[2,2], ps[2,1], rownames(ps)[2], pos = 2) 83 84 85################################################### 86### code chunk number 9: hac-data 87################################################### 88data("Investment") 89 90 91################################################### 92### code chunk number 10: hac-model 93################################################### 94fm.inv <- lm(RealInv ~ RealGNP + RealInt, data = Investment) 95 96 97################################################### 98### code chunk number 11: hac-test1 99################################################### 100coeftest(fm.inv, df = Inf, vcov = NeweyWest(fm.inv, lag = 4, prewhite = FALSE)) 101 102 103################################################### 104### code chunk number 12: hac-test2 105################################################### 106coeftest(fm.inv, df = Inf, vcov = NeweyWest) 107 108 109################################################### 110### code chunk number 13: hac-test3 111################################################### 112parzenHAC <- function(x, ...) kernHAC(x, kernel = "Parzen", prewhite = 2, 113 adjust = FALSE, bw = bwNeweyWest, ...) 114coeftest(fm.inv, df = Inf, vcov = parzenHAC) 115 116 117################################################### 118### code chunk number 14: hac-plot 119################################################### 120library("scatterplot3d") 121s3d <- scatterplot3d(Investment[,c(5,7,6)], 122 type = "b", angle = 65, scale.y = 1, pch = 16) 123s3d$plane3d(fm.inv, lty.box = "solid", col = 4) 124 125 126################################################### 127### code chunk number 15: loadlibs2 128################################################### 129library("strucchange") 130data("RealInt") 131 132 133################################################### 134### code chunk number 16: sc-ocus 135################################################### 136ocus <- gefp(RealInt ~ 1, fit = lm, vcov = kernHAC) 137 138 139################################################### 140### code chunk number 17: sc-bp 141################################################### 142bp <- breakpoints(RealInt ~ 1) 143confint(bp, vcov = kernHAC) 144 145 146################################################### 147### code chunk number 18: sc-plot 148################################################### 149par(mfrow = c(1, 2)) 150plot(ocus, aggregate = FALSE, main = "") 151plot(RealInt, ylab = "Real interest rate") 152lines(ts(fitted(bp), start = start(RealInt), freq = 4), col = 4) 153lines(confint(bp, vcov = kernHAC)) 154 155 156################################################### 157### code chunk number 19: sandwich.Rnw:786-787 158################################################### 159options(prompt = " ") 160 161 162################################################### 163### code chunk number 20: sandwich.Rnw:805-807 (eval = FALSE) 164################################################### 165## library("sandwich") 166## library("lmtest") 167## library("strucchange") 168 169 170################################################### 171### code chunk number 21: sandwich.Rnw:814-815 (eval = FALSE) 172################################################### 173## data("PublicSchools") 174## ps <- na.omit(PublicSchools) 175## ps$Income <- ps$Income * 0.0001 176 177 178################################################### 179### code chunk number 22: sandwich.Rnw:819-820 (eval = FALSE) 180################################################### 181## fm.ps <- lm(Expenditure ~ Income + I(Income^2), data = ps) 182 183 184################################################### 185### code chunk number 23: sandwich.Rnw:824-829 (eval = FALSE) 186################################################### 187## sqrt(diag(vcov(fm.ps))) 188## sqrt(diag(vcovHC(fm.ps, type = "const"))) 189## sqrt(diag(vcovHC(fm.ps, type = "HC0"))) 190## sqrt(diag(vcovHC(fm.ps, type = "HC3"))) 191## sqrt(diag(vcovHC(fm.ps, type = "HC4"))) 192 193 194################################################### 195### code chunk number 24: sandwich.Rnw:833-835 (eval = FALSE) 196################################################### 197## coeftest(fm.ps, df = Inf, vcov = vcovHC(fm.ps, type = "HC0")) 198## coeftest(fm.ps, df = Inf, vcov = vcovHC(fm.ps, type = "HC4")) 199 200 201################################################### 202### code chunk number 25: sandwich.Rnw:855-856 (eval = FALSE) 203################################################### 204## data("Investment") 205 206 207################################################### 208### code chunk number 26: sandwich.Rnw:860-861 (eval = FALSE) 209################################################### 210## fm.inv <- lm(RealInv ~ RealGNP + RealInt, data = Investment) 211 212 213################################################### 214### code chunk number 27: sandwich.Rnw:879-881 (eval = FALSE) 215################################################### 216## plot(Investment[, "RealInv"], type = "b", pch = 19, ylab = "Real investment") 217## lines(ts(fitted(fm.inv), start = 1964), col = 4) 218 219 220################################################### 221### code chunk number 28: sandwich.Rnw:897-898 (eval = FALSE) 222################################################### 223## data("RealInt") 224 225 226################################################### 227### code chunk number 29: sandwich.Rnw:902-905 (eval = FALSE) 228################################################### 229## ocus <- gefp(RealInt ~ 1, fit = lm, vcov = kernHAC) 230## plot(ocus, aggregate = FALSE) 231## sctest(ocus) 232 233 234################################################### 235### code chunk number 30: sandwich.Rnw:909-912 (eval = FALSE) 236################################################### 237## fs <- Fstats(RealInt ~ 1, vcov = kernHAC) 238## plot(fs) 239## sctest(fs) 240 241 242################################################### 243### code chunk number 31: sandwich.Rnw:917-919 (eval = FALSE) 244################################################### 245## bp <- breakpoints(RealInt ~ 1) 246## confint(bp, vcov = kernHAC) 247## plot(bp) 248 249 250################################################### 251### code chunk number 32: sandwich.Rnw:923-926 (eval = FALSE) 252################################################### 253## plot(RealInt, ylab = "Real interest rate") 254## lines(ts(fitted(bp), start = start(RealInt), freq = 4), col = 4) 255## lines(confint(bp, vcov = kernHAC)) 256 257 258