1###################################################
2### chunk number 1: setup
3###################################################
4options(prompt = "R> ", continue = "+  ", width = 64,
5  digits = 4, show.signif.stars = FALSE, useFancyQuotes = FALSE)
6
7options(SweaveHooks = list(onefig =   function() {par(mfrow = c(1,1))},
8                           twofig =   function() {par(mfrow = c(1,2))},
9                           threefig = function() {par(mfrow = c(1,3))},
10                           fourfig =  function() {par(mfrow = c(2,2))},
11                           sixfig =   function() {par(mfrow = c(3,2))}))
12
13library("AER")
14
15suppressWarnings(RNGversion("3.5.0"))
16set.seed(1071)
17
18
19###################################################
20### chunk number 2: ps-summary
21###################################################
22data("PublicSchools")
23summary(PublicSchools)
24
25
26###################################################
27### chunk number 3: ps-plot eval=FALSE
28###################################################
29## ps <- na.omit(PublicSchools)
30## ps$Income <- ps$Income / 10000
31## plot(Expenditure ~ Income, data = ps, ylim = c(230, 830))
32## ps_lm <- lm(Expenditure ~ Income, data = ps)
33## abline(ps_lm)
34## id <- c(2, 24, 48)
35## text(ps[id, 2:1], rownames(ps)[id], pos = 1, xpd = TRUE)
36
37
38###################################################
39### chunk number 4: ps-plot1
40###################################################
41ps <- na.omit(PublicSchools)
42ps$Income <- ps$Income / 10000
43plot(Expenditure ~ Income, data = ps, ylim = c(230, 830))
44ps_lm <- lm(Expenditure ~ Income, data = ps)
45abline(ps_lm)
46id <- c(2, 24, 48)
47text(ps[id, 2:1], rownames(ps)[id], pos = 1, xpd = TRUE)
48
49
50###################################################
51### chunk number 5: ps-lmplot eval=FALSE
52###################################################
53## plot(ps_lm, which = 1:6)
54
55
56###################################################
57### chunk number 6: ps-lmplot1
58###################################################
59plot(ps_lm, which = 1:6)
60
61
62###################################################
63### chunk number 7: ps-hatvalues eval=FALSE
64###################################################
65## ps_hat <- hatvalues(ps_lm)
66## plot(ps_hat)
67## abline(h = c(1, 3) * mean(ps_hat), col = 2)
68## id <- which(ps_hat > 3 * mean(ps_hat))
69## text(id, ps_hat[id], rownames(ps)[id], pos = 1, xpd = TRUE)
70
71
72###################################################
73### chunk number 8: ps-hatvalues1
74###################################################
75ps_hat <- hatvalues(ps_lm)
76plot(ps_hat)
77abline(h = c(1, 3) * mean(ps_hat), col = 2)
78id <- which(ps_hat > 3 * mean(ps_hat))
79text(id, ps_hat[id], rownames(ps)[id], pos = 1, xpd = TRUE)
80
81
82###################################################
83### chunk number 9: influence-measures1 eval=FALSE
84###################################################
85## influence.measures(ps_lm)
86
87
88###################################################
89### chunk number 10: which-hatvalues
90###################################################
91which(ps_hat > 3 * mean(ps_hat))
92
93
94###################################################
95### chunk number 11: influence-measures2
96###################################################
97summary(influence.measures(ps_lm))
98
99
100###################################################
101### chunk number 12: ps-noinf eval=FALSE
102###################################################
103## plot(Expenditure ~ Income, data = ps, ylim = c(230, 830))
104## abline(ps_lm)
105## id <- which(apply(influence.measures(ps_lm)$is.inf, 1, any))
106## text(ps[id, 2:1], rownames(ps)[id], pos = 1, xpd = TRUE)
107## ps_noinf <- lm(Expenditure ~ Income, data = ps[-id,])
108## abline(ps_noinf, lty = 2)
109
110
111###################################################
112### chunk number 13: ps-noinf1
113###################################################
114plot(Expenditure ~ Income, data = ps, ylim = c(230, 830))
115abline(ps_lm)
116id <- which(apply(influence.measures(ps_lm)$is.inf, 1, any))
117text(ps[id, 2:1], rownames(ps)[id], pos = 1, xpd = TRUE)
118ps_noinf <- lm(Expenditure ~ Income, data = ps[-id,])
119abline(ps_noinf, lty = 2)
120
121
122###################################################
123### chunk number 14: journals-age
124###################################################
125data("Journals")
126journals <- Journals[, c("subs", "price")]
127journals$citeprice <- Journals$price/Journals$citations
128journals$age <- 2000 - Journals$foundingyear
129
130
131###################################################
132### chunk number 15: journals-lm
133###################################################
134jour_lm <- lm(log(subs) ~ log(citeprice), data = journals)
135
136
137###################################################
138### chunk number 16: bptest1
139###################################################
140bptest(jour_lm)
141
142
143###################################################
144### chunk number 17: bptest2
145###################################################
146bptest(jour_lm, ~ log(citeprice) + I(log(citeprice)^2),
147  data = journals)
148
149
150###################################################
151### chunk number 18: gqtest
152###################################################
153gqtest(jour_lm, order.by = ~ citeprice, data = journals)
154
155
156###################################################
157### chunk number 19: resettest
158###################################################
159resettest(jour_lm)
160
161
162###################################################
163### chunk number 20: raintest
164###################################################
165raintest(jour_lm, order.by = ~ age, data = journals)
166
167
168###################################################
169### chunk number 21: harvtest
170###################################################
171harvtest(jour_lm, order.by = ~ age, data = journals)
172
173
174###################################################
175### chunk number 22:
176###################################################
177library("dynlm")
178
179
180###################################################
181### chunk number 23: usmacro-dynlm
182###################################################
183data("USMacroG")
184consump1 <- dynlm(consumption ~ dpi + L(dpi),
185  data = USMacroG)
186
187
188###################################################
189### chunk number 24: dwtest
190###################################################
191dwtest(consump1)
192
193
194###################################################
195### chunk number 25: Box-test
196###################################################
197Box.test(residuals(consump1), type = "Ljung-Box")
198
199
200###################################################
201### chunk number 26: bgtest
202###################################################
203bgtest(consump1)
204
205
206###################################################
207### chunk number 27: vcov
208###################################################
209vcov(jour_lm)
210vcovHC(jour_lm)
211
212
213###################################################
214### chunk number 28: coeftest
215###################################################
216coeftest(jour_lm, vcov = vcovHC)
217
218
219###################################################
220### chunk number 29: sandwiches
221###################################################
222t(sapply(c("const", "HC0", "HC1", "HC2", "HC3", "HC4"),
223  function(x) sqrt(diag(vcovHC(jour_lm, type = x)))))
224
225
226###################################################
227### chunk number 30: ps-anova
228###################################################
229ps_lm <- lm(Expenditure ~ Income, data = ps)
230ps_lm2 <- lm(Expenditure ~ Income + I(Income^2), data = ps)
231anova(ps_lm, ps_lm2)
232
233
234###################################################
235### chunk number 31: ps-waldtest
236###################################################
237waldtest(ps_lm, ps_lm2, vcov = vcovHC(ps_lm2, type = "HC4"))
238
239
240###################################################
241### chunk number 32: vcovHAC
242###################################################
243rbind(SE = sqrt(diag(vcov(consump1))),
244  QS = sqrt(diag(kernHAC(consump1))),
245  NW = sqrt(diag(NeweyWest(consump1))))
246
247
248###################################################
249### chunk number 33: solow-lm
250###################################################
251data("OECDGrowth")
252solow_lm <- lm(log(gdp85/gdp60) ~ log(gdp60) +
253  log(invest) + log(popgrowth + .05), data = OECDGrowth)
254summary(solow_lm)
255
256
257###################################################
258### chunk number 34: solow-plot eval=FALSE
259###################################################
260## plot(solow_lm)
261
262
263###################################################
264### chunk number 35: solow-lts
265###################################################
266library("MASS")
267solow_lts <- lqs(log(gdp85/gdp60) ~ log(gdp60) +
268  log(invest) + log(popgrowth + .05), data = OECDGrowth,
269  psamp = 13, nsamp = "exact")
270
271
272###################################################
273### chunk number 36: solow-smallresid
274###################################################
275smallresid <- which(
276  abs(residuals(solow_lts)/solow_lts$scale[2]) <= 2.5)
277
278
279###################################################
280### chunk number 37: solow-nohighlev
281###################################################
282X <- model.matrix(solow_lm)[,-1]
283Xcv <- cov.rob(X, nsamp = "exact")
284nohighlev <- which(
285  sqrt(mahalanobis(X, Xcv$center, Xcv$cov)) <= 2.5)
286
287
288###################################################
289### chunk number 38: solow-goodobs
290###################################################
291goodobs <- unique(c(smallresid, nohighlev))
292
293
294###################################################
295### chunk number 39: solow-badobs
296###################################################
297rownames(OECDGrowth)[-goodobs]
298
299
300###################################################
301### chunk number 40: solow-rob
302###################################################
303solow_rob <- update(solow_lm, subset = goodobs)
304summary(solow_rob)
305
306
307###################################################
308### chunk number 41: quantreg
309###################################################
310library("quantreg")
311
312
313###################################################
314### chunk number 42: cps-lad
315###################################################
316library("quantreg")
317data("CPS1988")
318cps_f <- log(wage) ~ experience + I(experience^2) + education
319cps_lad <- rq(cps_f, data = CPS1988)
320summary(cps_lad)
321
322
323###################################################
324### chunk number 43: cps-rq
325###################################################
326cps_rq <- rq(cps_f, tau = c(0.25, 0.75), data = CPS1988)
327summary(cps_rq)
328
329
330###################################################
331### chunk number 44: cps-rqs
332###################################################
333cps_rq25 <- rq(cps_f, tau = 0.25, data = CPS1988)
334cps_rq75 <- rq(cps_f, tau = 0.75, data = CPS1988)
335anova(cps_rq25, cps_rq75)
336
337
338###################################################
339### chunk number 45: cps-rq-anova
340###################################################
341anova(cps_rq25, cps_rq75, joint = FALSE)
342
343
344###################################################
345### chunk number 46: rqbig
346###################################################
347cps_rqbig <- rq(cps_f, tau = seq(0.05, 0.95, by = 0.05),
348  data = CPS1988)
349cps_rqbigs <- summary(cps_rqbig)
350
351
352###################################################
353### chunk number 47: rqbig-plot eval=FALSE
354###################################################
355## plot(cps_rqbigs)
356
357
358###################################################
359### chunk number 48: rqbig-plot1
360###################################################
361plot(cps_rqbigs)
362
363
364