## Exercise Chapter 3
## Problem 3
## Background
install.packages(“AER”)
suppressMessages(library(“AER”))
## attach the data-set Journals to the current R-session
data(“Journals”, package = “AER”)
## ?Journals # Check the help file # 一组数据，美国图书馆经济期刊的订阅量
## Select variable “subs” and “price”
journals <- Journals[, c("subs", "price")] # sub为订阅量，price为订阅价格；方括号为select；逗号前是行，逗号后是列，这里选中列
## Define variable 'journals-price per citation'
journals$citeprice <- Journals$price/Journals$citations # 定义列；$用来引用Journals中的数据
## Define variable 'journals-age'
journals$age <- 2020 - Journals$foundingyear # 定义列；foundingyear为期刊成立年份
## Check variable name in 'journals'
names(journals)
## [1] "subs" "price" "citeprice" "age"
## Estimate the coefficients beta_1 and beta_2 of the linear regression model: log(Y) = beta_1 + beta_2*log(X) + eps
## log(Y) = log(subs) and log(X) = log(citeprice)
## Problem 3a 画residuals - fitted value图看是否有heteroscedastic error-term variances(异方差)
## Solution
jour_lm <- lm(log(subs) ~ log(citeprice), data = journals) # lm 用于拟合线性模型，格式 lm(Y~X1+...+Xk, data), Y~X1+...+Xk 为通过X预测Y，其中+用于分隔预测变量
## Diagnostic plot residuals against fitted values
## plot(y = resid(jour_lm), x = fitted(jour_lm)) # 朴素画法
## Or 高级画法
plot(jour_lm, which = 1) # which是图的类型
## 图中点的分布明显左边更高，右边更低，说明不均匀，有different variance in the error term，所以说是heteroscedastic
## Problem 3b Estimate the standard error of the OLS estimator beta_hat_2
## Solution
## HC robust variance estimation
library("sandwich")
## Robust estimation of the variance of \hat\beta
Var_hat_beta_HC3 <- sandwich::vcovHC(jour_lm, type = "HC3") # sandwich::vcovHC为Heteroscedasticity-Consistent Covariance Matrix Estimation；HC3最常用，作业都用了HC3
# type为estimation type
## Robust standard error of \hat\beta_2
sqrt(diag(Var_hat_beta_HC3)[2]) # sqrt开方；diag(x)用x中元素作为主对角元素创造对角矩阵；方括号选出2代表2. variance component
## log(citeprice)
## 0.03447364
## Comparison with the classic standard error estimation
sqrt(diag(vcov(jour_lm))[2]) # Calculate Variance-Covariance Matrix for a Fitted Model Object
## log(citeprice)
## 0.0356132