R version 4.2.0 (2022-04-22) -- "Vigorous Calisthenics"
Copyright (C) 2022 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> if(!requireNamespace("splines") ||
+ !requireNamespace("dynlm") ||
+ !requireNamespace("plm") ||
+ !requireNamespace("systemfit") ||
+ !requireNamespace("nlme")) q()
Loading required namespace: splines
Loading required namespace: dynlm
Loading required namespace: plm
Loading required namespace: systemfit
>
> ###################################################
> ### chunk number 1: setup
> ###################################################
> options(prompt = "R> ", continue = "+ ", width = 64,
+ digits = 4, show.signif.stars = FALSE, useFancyQuotes = FALSE)
R>
R> options(SweaveHooks = list(onefig = function() {par(mfrow = c(1,1))},
+ twofig = function() {par(mfrow = c(1,2))},
+ threefig = function() {par(mfrow = c(1,3))},
+ fourfig = function() {par(mfrow = c(2,2))},
+ sixfig = function() {par(mfrow = c(3,2))}))
R>
R> library("AER")
Loading required package: car
Loading required package: carData
Loading required package: lmtest
Loading required package: zoo
Attaching package: 'zoo'
The following objects are masked from 'package:base':
as.Date, as.Date.numeric
Loading required package: sandwich
Loading required package: survival
R>
R> suppressWarnings(RNGversion("3.5.0"))
R> set.seed(1071)
R>
R>
R> ###################################################
R> ### chunk number 2: data-journals
R> ###################################################
R> data("Journals")
R> journals <- Journals[, c("subs", "price")]
R> journals$citeprice <- Journals$price/Journals$citations
R> summary(journals)
subs price citeprice
Min. : 2 Min. : 20 Min. : 0.005
1st Qu.: 52 1st Qu.: 134 1st Qu.: 0.464
Median : 122 Median : 282 Median : 1.321
Mean : 197 Mean : 418 Mean : 2.548
3rd Qu.: 268 3rd Qu.: 541 3rd Qu.: 3.440
Max. :1098 Max. :2120 Max. :24.459
R>
R>
R> ###################################################
R> ### chunk number 3: linreg-plot eval=FALSE
R> ###################################################
R> ## plot(log(subs) ~ log(citeprice), data = journals)
R> ## jour_lm <- lm(log(subs) ~ log(citeprice), data = journals)
R> ## abline(jour_lm)
R>
R>
R> ###################################################
R> ### chunk number 4: linreg-plot1
R> ###################################################
R> plot(log(subs) ~ log(citeprice), data = journals)
R> jour_lm <- lm(log(subs) ~ log(citeprice), data = journals)
R> abline(jour_lm)
R>
R>
R> ###################################################
R> ### chunk number 5: linreg-class
R> ###################################################
R> class(jour_lm)
[1] "lm"
R>
R>
R> ###################################################
R> ### chunk number 6: linreg-names
R> ###################################################
R> names(jour_lm)
[1] "coefficients" "residuals" "effects"
[4] "rank" "fitted.values" "assign"
[7] "qr" "df.residual" "xlevels"
[10] "call" "terms" "model"
R>
R>
R> ###################################################
R> ### chunk number 7: linreg-summary
R> ###################################################
R> summary(jour_lm)
Call:
lm(formula = log(subs) ~ log(citeprice), data = journals)
Residuals:
Min 1Q Median 3Q Max
-2.7248 -0.5361 0.0372 0.4662 1.8481
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.7662 0.0559 85.2 <2e-16
log(citeprice) -0.5331 0.0356 -15.0 <2e-16
Residual standard error: 0.75 on 178 degrees of freedom
Multiple R-squared: 0.557, Adjusted R-squared: 0.555
F-statistic: 224 on 1 and 178 DF, p-value: <2e-16
R>
R>
R> ###################################################
R> ### chunk number 8: linreg-summary
R> ###################################################
R> jour_slm <- summary(jour_lm)
R> class(jour_slm)
[1] "summary.lm"
R> names(jour_slm)
[1] "call" "terms" "residuals"
[4] "coefficients" "aliased" "sigma"
[7] "df" "r.squared" "adj.r.squared"
[10] "fstatistic" "cov.unscaled"
R>
R>
R> ###################################################
R> ### chunk number 9: linreg-coef
R> ###################################################
R> jour_slm$coefficients
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.7662 0.05591 85.25 2.954e-146
log(citeprice) -0.5331 0.03561 -14.97 2.564e-33
R>
R>
R> ###################################################
R> ### chunk number 10: linreg-anova
R> ###################################################
R> anova(jour_lm)
Analysis of Variance Table
Response: log(subs)
Df Sum Sq Mean Sq F value Pr(>F)
log(citeprice) 1 126 125.9 224 <2e-16
Residuals 178 100 0.6
R>
R>
R> ###################################################
R> ### chunk number 11: journals-coef
R> ###################################################
R> coef(jour_lm)
(Intercept) log(citeprice)
4.7662 -0.5331
R>
R>
R> ###################################################
R> ### chunk number 12: journals-confint
R> ###################################################
R> confint(jour_lm, level = 0.95)
2.5 % 97.5 %
(Intercept) 4.6559 4.8765
log(citeprice) -0.6033 -0.4628
R>
R>
R> ###################################################
R> ### chunk number 13: journals-predict
R> ###################################################
R> predict(jour_lm, newdata = data.frame(citeprice = 2.11),
+ interval = "confidence")
fit lwr upr
1 4.368 4.247 4.489
R> predict(jour_lm, newdata = data.frame(citeprice = 2.11),
+ interval = "prediction")
fit lwr upr
1 4.368 2.884 5.853
R>
R>
R> ###################################################
R> ### chunk number 14: predict-plot eval=FALSE
R> ###################################################
R> ## lciteprice <- seq(from = -6, to = 4, by = 0.25)
R> ## jour_pred <- predict(jour_lm, interval = "prediction",
R> ## newdata = data.frame(citeprice = exp(lciteprice)))
R> ## plot(log(subs) ~ log(citeprice), data = journals)
R> ## lines(jour_pred[, 1] ~ lciteprice, col = 1)
R> ## lines(jour_pred[, 2] ~ lciteprice, col = 1, lty = 2)
R> ## lines(jour_pred[, 3] ~ lciteprice, col = 1, lty = 2)
R>
R>
R> ###################################################
R> ### chunk number 15: predict-plot1
R> ###################################################
R> lciteprice <- seq(from = -6, to = 4, by = 0.25)
R> jour_pred <- predict(jour_lm, interval = "prediction",
+ newdata = data.frame(citeprice = exp(lciteprice)))
R> plot(log(subs) ~ log(citeprice), data = journals)
R> lines(jour_pred[, 1] ~ lciteprice, col = 1)
R> lines(jour_pred[, 2] ~ lciteprice, col = 1, lty = 2)
R> lines(jour_pred[, 3] ~ lciteprice, col = 1, lty = 2)
R>
R>
R> ###################################################
R> ### chunk number 16: journals-plot eval=FALSE
R> ###################################################
R> ## par(mfrow = c(2, 2))
R> ## plot(jour_lm)
R> ## par(mfrow = c(1, 1))
R>
R>
R> ###################################################
R> ### chunk number 17: journals-plot1
R> ###################################################
R> par(mfrow = c(2, 2))
R> plot(jour_lm)
R> par(mfrow = c(1, 1))
R>
R>
R> ###################################################
R> ### chunk number 18: journal-lht
R> ###################################################
R> linearHypothesis(jour_lm, "log(citeprice) = -0.5")
Linear hypothesis test
Hypothesis:
log(citeprice) = - 0.5
Model 1: restricted model
Model 2: log(subs) ~ log(citeprice)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 179 100
2 178 100 1 0.484 0.86 0.35
R>
R>
R> ###################################################
R> ### chunk number 19: CPS-data
R> ###################################################
R> data("CPS1988")
R> summary(CPS1988)
wage education experience ethnicity
Min. : 50 Min. : 0.0 Min. :-4.0 cauc:25923
1st Qu.: 309 1st Qu.:12.0 1st Qu.: 8.0 afam: 2232
Median : 522 Median :12.0 Median :16.0
Mean : 604 Mean :13.1 Mean :18.2
3rd Qu.: 783 3rd Qu.:15.0 3rd Qu.:27.0
Max. :18777 Max. :18.0 Max. :63.0
smsa region parttime
no : 7223 northeast:6441 no :25631
yes:20932 midwest :6863 yes: 2524
south :8760
west :6091
R>
R>
R> ###################################################
R> ### chunk number 20: CPS-base
R> ###################################################
R> cps_lm <- lm(log(wage) ~ experience + I(experience^2) +
+ education + ethnicity, data = CPS1988)
R>
R>
R> ###################################################
R> ### chunk number 21: CPS-visualization-unused eval=FALSE
R> ###################################################
R> ## ex <- 0:56
R> ## ed <- with(CPS1988, tapply(education,
R> ## list(ethnicity, experience), mean))[, as.character(ex)]
R> ## fm <- cps_lm
R> ## wago <- predict(fm, newdata = data.frame(experience = ex,
R> ## ethnicity = "cauc", education = as.numeric(ed["cauc",])))
R> ## wagb <- predict(fm, newdata = data.frame(experience = ex,
R> ## ethnicity = "afam", education = as.numeric(ed["afam",])))
R> ## plot(log(wage) ~ experience, data = CPS1988, pch = ".",
R> ## col = as.numeric(ethnicity))
R> ## lines(ex, wago)
R> ## lines(ex, wagb, col = 2)
R>
R>
R> ###################################################
R> ### chunk number 22: CPS-summary
R> ###################################################
R> summary(cps_lm)
Call:
lm(formula = log(wage) ~ experience + I(experience^2) + education +
ethnicity, data = CPS1988)
Residuals:
Min 1Q Median 3Q Max
-2.943 -0.316 0.058 0.376 4.383
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.321395 0.019174 225.4 <2e-16
experience 0.077473 0.000880 88.0 <2e-16
I(experience^2) -0.001316 0.000019 -69.3 <2e-16
education 0.085673 0.001272 67.3 <2e-16
ethnicityafam -0.243364 0.012918 -18.8 <2e-16
Residual standard error: 0.584 on 28150 degrees of freedom
Multiple R-squared: 0.335, Adjusted R-squared: 0.335
F-statistic: 3.54e+03 on 4 and 28150 DF, p-value: <2e-16
R>
R>
R> ###################################################
R> ### chunk number 23: CPS-noeth
R> ###################################################
R> cps_noeth <- lm(log(wage) ~ experience + I(experience^2) +
+ education, data = CPS1988)
R> anova(cps_noeth, cps_lm)
Analysis of Variance Table
Model 1: log(wage) ~ experience + I(experience^2) + education
Model 2: log(wage) ~ experience + I(experience^2) + education + ethnicity
Res.Df RSS Df Sum of Sq F Pr(>F)
1 28151 9720
2 28150 9599 1 121 355 <2e-16
R>
R>
R> ###################################################
R> ### chunk number 24: CPS-anova
R> ###################################################
R> anova(cps_lm)
Analysis of Variance Table
Response: log(wage)
Df Sum Sq Mean Sq F value Pr(>F)
experience 1 840 840 2462 <2e-16
I(experience^2) 1 2249 2249 6597 <2e-16
education 1 1620 1620 4750 <2e-16
ethnicity 1 121 121 355 <2e-16
Residuals 28150 9599 0
R>
R>
R> ###################################################
R> ### chunk number 25: CPS-noeth2 eval=FALSE
R> ###################################################
R> ## cps_noeth <- update(cps_lm, formula = . ~ . - ethnicity)
R>
R>
R> ###################################################
R> ### chunk number 26: CPS-waldtest
R> ###################################################
R> waldtest(cps_lm, . ~ . - ethnicity)
Wald test
Model 1: log(wage) ~ experience + I(experience^2) + education + ethnicity
Model 2: log(wage) ~ experience + I(experience^2) + education
Res.Df Df F Pr(>F)
1 28150
2 28151 -1 355 <2e-16
R>
R>
R> ###################################################
R> ### chunk number 27: CPS-spline
R> ###################################################
R> library("splines")
R> cps_plm <- lm(log(wage) ~ bs(experience, df = 5) +
+ education + ethnicity, data = CPS1988)
R>
R>
R> ###################################################
R> ### chunk number 28: CPS-spline-summary eval=FALSE
R> ###################################################
R> ## summary(cps_plm)
R>
R>
R> ###################################################
R> ### chunk number 29: CPS-BIC
R> ###################################################
R> cps_bs <- lapply(3:10, function(i) lm(log(wage) ~
+ bs(experience, df = i) + education + ethnicity,
+ data = CPS1988))
R> structure(sapply(cps_bs, AIC, k = log(nrow(CPS1988))),
+ .Names = 3:10)
3 4 5 6 7 8 9 10
49205 48836 48794 48795 48801 48797 48799 48802
R>
R>
R> ###################################################
R> ### chunk number 30: plm-plot eval=FALSE
R> ###################################################
R> ## cps <- data.frame(experience = -2:60, education =
R> ## with(CPS1988, mean(education[ethnicity == "cauc"])),
R> ## ethnicity = "cauc")
R> ## cps$yhat1 <- predict(cps_lm, newdata = cps)
R> ## cps$yhat2 <- predict(cps_plm, newdata = cps)
R> ##
R> ## plot(log(wage) ~ jitter(experience, factor = 3), pch = 19,
R> ## col = rgb(0.5, 0.5, 0.5, alpha = 0.02), data = CPS1988)
R> ## lines(yhat1 ~ experience, data = cps, lty = 2)
R> ## lines(yhat2 ~ experience, data = cps)
R> ## legend("topleft", c("quadratic", "spline"), lty = c(2,1),
R> ## bty = "n")
R>
R>
R> ###################################################
R> ### chunk number 31: plm-plot1
R> ###################################################
R> cps <- data.frame(experience = -2:60, education =
+ with(CPS1988, mean(education[ethnicity == "cauc"])),
+ ethnicity = "cauc")
R> cps$yhat1 <- predict(cps_lm, newdata = cps)
R> cps$yhat2 <- predict(cps_plm, newdata = cps)
R>
R> plot(log(wage) ~ jitter(experience, factor = 3), pch = 19,
+ col = rgb(0.5, 0.5, 0.5, alpha = 0.02), data = CPS1988)
R> lines(yhat1 ~ experience, data = cps, lty = 2)
R> lines(yhat2 ~ experience, data = cps)
R> legend("topleft", c("quadratic", "spline"), lty = c(2,1),
+ bty = "n")
R>
R>
R> ###################################################
R> ### chunk number 32: CPS-int
R> ###################################################
R> cps_int <- lm(log(wage) ~ experience + I(experience^2) +
+ education * ethnicity, data = CPS1988)
R> coeftest(cps_int)
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.313059 0.019590 220.17 <2e-16
experience 0.077520 0.000880 88.06 <2e-16
I(experience^2) -0.001318 0.000019 -69.34 <2e-16
education 0.086312 0.001309 65.94 <2e-16
ethnicityafam -0.123887 0.059026 -2.10 0.036
education:ethnicityafam -0.009648 0.004651 -2.07 0.038
R>
R>
R> ###################################################
R> ### chunk number 33: CPS-int2 eval=FALSE
R> ###################################################
R> ## cps_int <- lm(log(wage) ~ experience + I(experience^2) +
R> ## education + ethnicity + education:ethnicity,
R> ## data = CPS1988)
R>
R>
R> ###################################################
R> ### chunk number 34: CPS-sep
R> ###################################################
R> cps_sep <- lm(log(wage) ~ ethnicity /
+ (experience + I(experience^2) + education) - 1,
+ data = CPS1988)
R>
R>
R> ###################################################
R> ### chunk number 35: CPS-sep-coef
R> ###################################################
R> cps_sep_cf <- matrix(coef(cps_sep), nrow = 2)
R> rownames(cps_sep_cf) <- levels(CPS1988$ethnicity)
R> colnames(cps_sep_cf) <- names(coef(cps_lm))[1:4]
R> cps_sep_cf
(Intercept) experience I(experience^2) education
cauc 4.310 0.07923 -0.0013597 0.08575
afam 4.159 0.06190 -0.0009415 0.08654
R>
R>
R> ###################################################
R> ### chunk number 36: CPS-sep-anova
R> ###################################################
R> anova(cps_sep, cps_lm)
Analysis of Variance Table
Model 1: log(wage) ~ ethnicity/(experience + I(experience^2) + education) -
1
Model 2: log(wage) ~ experience + I(experience^2) + education + ethnicity
Res.Df RSS Df Sum of Sq F Pr(>F)
1 28147 9582
2 28150 9599 -3 -16.8 16.5 1.1e-10
R>
R>
R> ###################################################
R> ### chunk number 37: CPS-sep-visualization-unused eval=FALSE
R> ###################################################
R> ## ex <- 0:56
R> ## ed <- with(CPS1988, tapply(education, list(ethnicity,
R> ## experience), mean))[, as.character(ex)]
R> ## fm <- cps_lm
R> ## wago <- predict(fm, newdata = data.frame(experience = ex,
R> ## ethnicity = "cauc", education = as.numeric(ed["cauc",])))
R> ## wagb <- predict(fm, newdata = data.frame(experience = ex,
R> ## ethnicity = "afam", education = as.numeric(ed["afam",])))
R> ## plot(log(wage) ~ jitter(experience, factor = 2),
R> ## data = CPS1988, pch = ".", col = as.numeric(ethnicity))
R> ##
R> ##
R> ## plot(log(wage) ~ as.factor(experience), data = CPS1988,
R> ## pch = ".")
R> ## lines(ex, wago, lwd = 2)
R> ## lines(ex, wagb, col = 2, lwd = 2)
R> ## fm <- cps_sep
R> ## wago <- predict(fm, newdata = data.frame(experience = ex,
R> ## ethnicity = "cauc", education = as.numeric(ed["cauc",])))
R> ## wagb <- predict(fm, newdata = data.frame(experience = ex,
R> ## ethnicity = "afam", education = as.numeric(ed["afam",])))
R> ## lines(ex, wago, lty = 2, lwd = 2)
R> ## lines(ex, wagb, col = 2, lty = 2, lwd = 2)
R>
R>
R> ###################################################
R> ### chunk number 38: CPS-region
R> ###################################################
R> CPS1988$region <- relevel(CPS1988$region, ref = "south")
R> cps_region <- lm(log(wage) ~ ethnicity + education +
+ experience + I(experience^2) + region, data = CPS1988)
R> coef(cps_region)
(Intercept) ethnicityafam education experience
4.283606 -0.225679 0.084672 0.077656
I(experience^2) regionnortheast regionmidwest regionwest
-0.001323 0.131920 0.043789 0.040327
R>
R>
R> ###################################################
R> ### chunk number 39: wls1
R> ###################################################
R> jour_wls1 <- lm(log(subs) ~ log(citeprice), data = journals,
+ weights = 1/citeprice^2)
R>
R>
R> ###################################################
R> ### chunk number 40: wls2
R> ###################################################
R> jour_wls2 <- lm(log(subs) ~ log(citeprice), data = journals,
+ weights = 1/citeprice)
R>
R>
R> ###################################################
R> ### chunk number 41: journals-wls1 eval=FALSE
R> ###################################################
R> ## plot(log(subs) ~ log(citeprice), data = journals)
R> ## abline(jour_lm)
R> ## abline(jour_wls1, lwd = 2, lty = 2)
R> ## abline(jour_wls2, lwd = 2, lty = 3)
R> ## legend("bottomleft", c("OLS", "WLS1", "WLS2"),
R> ## lty = 1:3, lwd = 2, bty = "n")
R>
R>
R> ###################################################
R> ### chunk number 42: journals-wls11
R> ###################################################
R> plot(log(subs) ~ log(citeprice), data = journals)
R> abline(jour_lm)
R> abline(jour_wls1, lwd = 2, lty = 2)
R> abline(jour_wls2, lwd = 2, lty = 3)
R> legend("bottomleft", c("OLS", "WLS1", "WLS2"),
+ lty = 1:3, lwd = 2, bty = "n")
R>
R>
R> ###################################################
R> ### chunk number 43: fgls1
R> ###################################################
R> auxreg <- lm(log(residuals(jour_lm)^2) ~ log(citeprice),
+ data = journals)
R> jour_fgls1 <- lm(log(subs) ~ log(citeprice),
+ weights = 1/exp(fitted(auxreg)), data = journals)
R>
R>
R> ###################################################
R> ### chunk number 44: fgls2
R> ###################################################
R> gamma2i <- coef(auxreg)[2]
R> gamma2 <- 0
R> while(abs((gamma2i - gamma2)/gamma2) > 1e-7) {
+ gamma2 <- gamma2i
+ fglsi <- lm(log(subs) ~ log(citeprice), data = journals,
+ weights = 1/citeprice^gamma2)
+ gamma2i <- coef(lm(log(residuals(fglsi)^2) ~
+ log(citeprice), data = journals))[2]
+ }
R> jour_fgls2 <- lm(log(subs) ~ log(citeprice), data = journals,
+ weights = 1/citeprice^gamma2)
R>
R>
R> ###################################################
R> ### chunk number 45: fgls2-coef
R> ###################################################
R> coef(jour_fgls2)
(Intercept) log(citeprice)
4.7758 -0.5008
R>
R>
R> ###################################################
R> ### chunk number 46: journals-fgls
R> ###################################################
R> plot(log(subs) ~ log(citeprice), data = journals)
R> abline(jour_lm)
R> abline(jour_fgls2, lty = 2, lwd = 2)
R>
R>
R> ###################################################
R> ### chunk number 47: usmacro-plot eval=FALSE
R> ###################################################
R> ## data("USMacroG")
R> ## plot(USMacroG[, c("dpi", "consumption")], lty = c(3, 1),
R> ## plot.type = "single", ylab = "")
R> ## legend("topleft", legend = c("income", "consumption"),
R> ## lty = c(3, 1), bty = "n")
R>
R>
R> ###################################################
R> ### chunk number 48: usmacro-plot1
R> ###################################################
R> data("USMacroG")
R> plot(USMacroG[, c("dpi", "consumption")], lty = c(3, 1),
+ plot.type = "single", ylab = "")
R> legend("topleft", legend = c("income", "consumption"),
+ lty = c(3, 1), bty = "n")
R>
R>
R> ###################################################
R> ### chunk number 49: usmacro-fit
R> ###################################################
R> library("dynlm")
R> cons_lm1 <- dynlm(consumption ~ dpi + L(dpi), data = USMacroG)
R> cons_lm2 <- dynlm(consumption ~ dpi + L(consumption),
+ data = USMacroG)
R>
R>
R> ###################################################
R> ### chunk number 50: usmacro-summary1
R> ###################################################
R> summary(cons_lm1)
Time series regression with "ts" data:
Start = 1950(2), End = 2000(4)
Call:
dynlm(formula = consumption ~ dpi + L(dpi), data = USMacroG)
Residuals:
Min 1Q Median 3Q Max
-190.0 -56.7 1.6 49.9 323.9
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -81.0796 14.5081 -5.59 7.4e-08
dpi 0.8912 0.2063 4.32 2.4e-05
L(dpi) 0.0309 0.2075 0.15 0.88
Residual standard error: 87.6 on 200 degrees of freedom
Multiple R-squared: 0.996, Adjusted R-squared: 0.996
F-statistic: 2.79e+04 on 2 and 200 DF, p-value: <2e-16
R>
R>
R> ###################################################
R> ### chunk number 51: usmacro-summary2
R> ###################################################
R> summary(cons_lm2)
Time series regression with "ts" data:
Start = 1950(2), End = 2000(4)
Call:
dynlm(formula = consumption ~ dpi + L(consumption), data = USMacroG)
Residuals:
Min 1Q Median 3Q Max
-101.30 -9.67 1.14 12.69 45.32
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.53522 3.84517 0.14 0.89
dpi -0.00406 0.01663 -0.24 0.81
L(consumption) 1.01311 0.01816 55.79 <2e-16
Residual standard error: 21.5 on 200 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 4.63e+05 on 2 and 200 DF, p-value: <2e-16
R>
R>
R> ###################################################
R> ### chunk number 52: dynlm-plot eval=FALSE
R> ###################################################
R> ## plot(merge(as.zoo(USMacroG[,"consumption"]), fitted(cons_lm1),
R> ## fitted(cons_lm2), 0, residuals(cons_lm1),
R> ## residuals(cons_lm2)), screens = rep(1:2, c(3, 3)),
R> ## lty = rep(1:3, 2), ylab = c("Fitted values", "Residuals"),
R> ## xlab = "Time", main = "")
R> ## legend(0.05, 0.95, c("observed", "cons_lm1", "cons_lm2"),
R> ## lty = 1:3, bty = "n")
R>
R>
R> ###################################################
R> ### chunk number 53: dynlm-plot1
R> ###################################################
R> plot(merge(as.zoo(USMacroG[,"consumption"]), fitted(cons_lm1),
+ fitted(cons_lm2), 0, residuals(cons_lm1),
+ residuals(cons_lm2)), screens = rep(1:2, c(3, 3)),
+ lty = rep(1:3, 2), ylab = c("Fitted values", "Residuals"),
+ xlab = "Time", main = "")
R> legend(0.05, 0.95, c("observed", "cons_lm1", "cons_lm2"),
+ lty = 1:3, bty = "n")
R>
R>
R> ###################################################
R> ### chunk number 54: encompassing1
R> ###################################################
R> cons_lmE <- dynlm(consumption ~ dpi + L(dpi) +
+ L(consumption), data = USMacroG)
R>
R>
R> ###################################################
R> ### chunk number 55: encompassing2
R> ###################################################
R> anova(cons_lm1, cons_lmE, cons_lm2)
Analysis of Variance Table
Model 1: consumption ~ dpi + L(dpi)
Model 2: consumption ~ dpi + L(dpi) + L(consumption)
Model 3: consumption ~ dpi + L(consumption)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 200 1534001
2 199 73550 1 1460451 3951.4 < 2e-16
3 200 92644 -1 -19094 51.7 1.3e-11
R>
R>
R> ###################################################
R> ### chunk number 56: encompassing3
R> ###################################################
R> encomptest(cons_lm1, cons_lm2)
Encompassing test
Model 1: consumption ~ dpi + L(dpi)
Model 2: consumption ~ dpi + L(consumption)
Model E: consumption ~ dpi + L(dpi) + L(consumption)
Res.Df Df F Pr(>F)
M1 vs. ME 199 -1 3951.4 < 2e-16
M2 vs. ME 199 -1 51.7 1.3e-11
R>
R>
R> ###################################################
R> ### chunk number 57: pdata.frame
R> ###################################################
R> data("Grunfeld", package = "AER")
R> library("plm")
R> gr <- subset(Grunfeld, firm %in% c("General Electric",
+ "General Motors", "IBM"))
R> pgr <- pdata.frame(gr, index = c("firm", "year"))
R>
R>
R> ###################################################
R> ### chunk number 58: plm-pool
R> ###################################################
R> gr_pool <- plm(invest ~ value + capital, data = pgr,
+ model = "pooling")
R>
R>
R> ###################################################
R> ### chunk number 59: plm-FE
R> ###################################################
R> gr_fe <- plm(invest ~ value + capital, data = pgr,
+ model = "within")
R> summary(gr_fe)
Oneway (individual) effect Within Model
Call:
plm(formula = invest ~ value + capital, data = pgr, model = "within")
Balanced Panel: n = 3, T = 20, N = 60
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-167.33 -26.14 2.09 26.84 201.68
Coefficients:
Estimate Std. Error t-value Pr(>|t|)
value 0.1049 0.0163 6.42 3.3e-08
capital 0.3453 0.0244 14.16 < 2e-16
Total Sum of Squares: 1890000
Residual Sum of Squares: 244000
R-Squared: 0.871
Adj. R-Squared: 0.861
F-statistic: 185.407 on 2 and 55 DF, p-value: <2e-16
R>
R>
R> ###################################################
R> ### chunk number 60: plm-pFtest
R> ###################################################
R> pFtest(gr_fe, gr_pool)
F test for individual effects
data: invest ~ value + capital
F = 57, df1 = 2, df2 = 55, p-value = 4e-14
alternative hypothesis: significant effects
R>
R>
R> ###################################################
R> ### chunk number 61: plm-RE
R> ###################################################
R> gr_re <- plm(invest ~ value + capital, data = pgr,
+ model = "random", random.method = "walhus")
R> summary(gr_re)
Oneway (individual) effect Random Effect Model
(Wallace-Hussain's transformation)
Call:
plm(formula = invest ~ value + capital, data = pgr, model = "random",
random.method = "walhus")
Balanced Panel: n = 3, T = 20, N = 60
Effects:
var std.dev share
idiosyncratic 4389.3 66.3 0.35
individual 8079.7 89.9 0.65
theta: 0.837
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-187.40 -32.92 6.96 31.43 210.20
Coefficients:
Estimate Std. Error z-value Pr(>|z|)
(Intercept) -109.9766 61.7014 -1.78 0.075
value 0.1043 0.0150 6.95 3.6e-12
capital 0.3448 0.0245 14.06 < 2e-16
Total Sum of Squares: 1990000
Residual Sum of Squares: 258000
R-Squared: 0.87
Adj. R-Squared: 0.866
Chisq: 383.089 on 2 DF, p-value: <2e-16
R>
R>
R> ###################################################
R> ### chunk number 62: plm-plmtest
R> ###################################################
R> plmtest(gr_pool)
Lagrange Multiplier Test - (Honda) for balanced panels
data: invest ~ value + capital
normal = 15, p-value <2e-16
alternative hypothesis: significant effects
R>
R>
R> ###################################################
R> ### chunk number 63: plm-phtest
R> ###################################################
R> phtest(gr_re, gr_fe)
Hausman Test
data: invest ~ value + capital
chisq = 0.04, df = 2, p-value = 1
alternative hypothesis: one model is inconsistent
R>
R>
R> ###################################################
R> ### chunk number 64: EmplUK-data
R> ###################################################
R> data("EmplUK", package = "plm")
R>
R>
R> ###################################################
R> ### chunk number 65: plm-AB
R> ###################################################
R> empl_ab <- pgmm(log(emp) ~ lag(log(emp), 1:2) + lag(log(wage), 0:1) +
+ log(capital) + lag(log(output), 0:1) | lag(log(emp), 2:99),
+ data = EmplUK, index = c("firm", "year"),
+ effect = "twoways", model = "twosteps")
R>
R>
R> ###################################################
R> ### chunk number 66: plm-AB-summary
R> ###################################################
R> summary(empl_ab, robust = FALSE)
Twoways effects Two-steps model Difference GMM
Call:
pgmm(formula = log(emp) ~ lag(log(emp), 1:2) + lag(log(wage),
0:1) + log(capital) + lag(log(output), 0:1) | lag(log(emp),
2:99), data = EmplUK, effect = "twoways", model = "twosteps",
index = c("firm", "year"))
Unbalanced Panel: n = 140, T = 7-9, N = 1031
Number of Observations Used: 611
Residuals:
Min. 1st Qu. Median Mean 3rd Qu. Max.
-0.6191 -0.0256 0.0000 -0.0001 0.0332 0.6410
Coefficients:
Estimate Std. Error z-value Pr(>|z|)
lag(log(emp), 1:2)1 0.4742 0.0853 5.56 2.7e-08
lag(log(emp), 1:2)2 -0.0530 0.0273 -1.94 0.05222
lag(log(wage), 0:1)0 -0.5132 0.0493 -10.40 < 2e-16
lag(log(wage), 0:1)1 0.2246 0.0801 2.81 0.00502
log(capital) 0.2927 0.0395 7.42 1.2e-13
lag(log(output), 0:1)0 0.6098 0.1085 5.62 1.9e-08
lag(log(output), 0:1)1 -0.4464 0.1248 -3.58 0.00035
Sargan test: chisq(25) = 30.11 (p-value = 0.22)
Autocorrelation test (1): normal = -2.428 (p-value = 0.0152)
Autocorrelation test (2): normal = -0.3325 (p-value = 0.739)
Wald test for coefficients: chisq(7) = 372 (p-value = <2e-16)
Wald test for time dummies: chisq(6) = 26.9 (p-value = 0.000151)
R>
R>
R> ###################################################
R> ### chunk number 67: systemfit
R> ###################################################
R> library("systemfit")
Loading required package: Matrix
Please cite the 'systemfit' package as:
Arne Henningsen and Jeff D. Hamann (2007). systemfit: A Package for Estimating Systems of Simultaneous Equations in R. Journal of Statistical Software 23(4), 1-40. http://www.jstatsoft.org/v23/i04/.
If you have questions, suggestions, or comments regarding the 'systemfit' package, please use a forum or 'tracker' at systemfit's R-Forge site:
https://r-forge.r-project.org/projects/systemfit/
R> gr2 <- subset(Grunfeld, firm %in% c("Chrysler", "IBM"))
R> pgr2 <- pdata.frame(gr2, c("firm", "year"))
R>
R>
R> ###################################################
R> ### chunk number 68: SUR
R> ###################################################
R> gr_sur <- systemfit(invest ~ value + capital,
+ method = "SUR", data = pgr2)
R> summary(gr_sur, residCov = FALSE, equations = FALSE)
systemfit results
method: SUR
N DF SSR detRCov OLS-R2 McElroy-R2
system 40 34 4114 11022 0.929 0.927
N DF SSR MSE RMSE R2 Adj R2
Chrysler 20 17 3002 176.6 13.29 0.913 0.903
IBM 20 17 1112 65.4 8.09 0.952 0.946
Coefficients:
Estimate Std. Error t value Pr(>|t|)
Chrysler_(Intercept) -5.7031 13.2774 -0.43 0.67293
Chrysler_value 0.0780 0.0196 3.98 0.00096
Chrysler_capital 0.3115 0.0287 10.85 4.6e-09
IBM_(Intercept) -8.0908 4.5216 -1.79 0.09139
IBM_value 0.1272 0.0306 4.16 0.00066
IBM_capital 0.0966 0.0983 0.98 0.33951
R>
R>
R> ###################################################
R> ### chunk number 69: nlme eval=FALSE
R> ###################################################
R> ## library("nlme")
R> ## g1 <- subset(Grunfeld, firm == "Westinghouse")
R> ## gls(invest ~ value + capital, data = g1, correlation = corAR1())
R>
R>
R>
> proc.time()
user system elapsed
3.003 0.104 3.099
Generated by dwww version 1.15 on Fri Jun 21 07:53:04 CEST 2024.