@fanxy
2016-07-14T20:39:16.000000Z
字数 2077
阅读 2610
樊潇彦
复旦大学经济学院
经济数学
# 准备工作
setwd("D:\\...")
install.packages("gvlma")
install.packages("car")
install.packages("sandwich")
install.packages("lmtest")
install.packages("AER")
Y=c(11,10,11,6,3)
P=c(1,1,0,1,0)
A=c(1,1,1,0,0)
college=data.frame(Y,P,A)
lm(Y~P+A,data=college)
我们假定真实参数 ,并模拟出 和 :
# 生成模拟数据
N=200 # 样本数
beta=c(1, 0.5, 0.1, 1) # 参数
K=length(beta)
x=matrix(rep(NaN,N*K),nrow=N) # N*K 空矩阵
mu_x=c(1,0,0,0); sd_x=c(0,1,2,3) # 均值和标准差
set.seed(1) # 设定随机种子使模拟结果可复制
for (k in 1:K) { # 模拟生成x和y
x[,k]=rnorm(N,mu_x[k],sd_x[k])
}
y=x%*%beta
colnames(x)=c("con","x1","x2","epsilon")
mydata=data.frame(y, x)
# 估计
x=x[,1:3]
beta_bar=solve(t(x)%*%x)%*%t(x)%*%y # 线性最小二乘估计量
myOLS=lm(y~x1+x2, data=mydata) # OLS 回归
summary(myOLS) # 报告OLS回归结果
如果某个解释变量是离散变量(如年龄组),想把它设为哑变量,就用factor()
。此外Kabacoff(2013)列出了 lm()
的其他常用设定,可以用来在解释变量中加入二次项、交互项,以及去掉某些解释变量:
gvlma
包中的 gvlma()
函数可以对模型假定进行综合检验,另外car
包则提供了对多重共线性vif()
、残差自相关durbinWatsonTest()
和异方差ncvTest()
等问题的检验命令。
library(gvlma)
gvlma(myOLS) # 模型假定是否可接受
library(car)
vif(myOLS)>4 # 检验多重共线性问题
durbinWatsonTest(myOLS) # 检验残差自相关问题
ncvTest(myOLS) # 检验异方差问题
spreadLevelPlot(myOLS) # 做图查看
library(sandwich)
sd=sqrt(diag(vcov(myOLS)))
sd_adjust=sqrt(diag(vcovHC(myOLS))) # 如果不仅有异方差,还有残差的自相关,用 vcovHAC() 代替 vcovHC()
data.frame(sd, sd_adjust) # 比较 beta 的标准差
library(lmtest)
coeftest(myOLS,vcov=vcovHC) # 报告调整标准差后的结果
在假定异方差和自相关的具体形式后,可以进行FGLS估计,具体方法是:(1)直接在lm()
中设定相应的weights
参数;(2)MASS
包中的lm.gls()
命令;(3)nlme
包中的gls()
命令。
myOLS
的回归中x1
有内生性,现在用x2
作为x1
的工具变量进行回归:
library(AER)
myiv = ivreg(formula = y~ x1 | x2, data=mydata)
summary(myiv)
上述结果等同于先回归 ,然后用 对 做OLS回归。下面报告的两阶段最小二乘回归(2SLS)的参数估计结果和前面一样,但标准差未经调整:
lm1=lm(x1~ x2, data=mydata)
x1_hat=fitted(lm1)
mydata=data.frame(mydata, x1_hat)
lm2SLS=lm(y~x1_hat, data=mydata)
summary(lm2SLS)