Dear R Users,
please find attached what I believe to be the solution to my problem. Note
that I am still not 100% sure if my approach really does what it is intended
to do and if it is applicable to my case at all.
Any comment or correction is highly appreciated.
Best wishes,
Philipp
##################################################################
##################################################################
##################################################################
#packages
library(lmtest)
#Data
y<-rnorm(100)
x1<-rnorm(100)
x2<-x1+rnorm(100)/5
x3<-rnorm(100)
#Simulated multicollinearity in the data
plot(x1,x2)
cor(x1,x2)
#I wish to estimate y=c+x1+x2+x3+u but cor(x1,x2) is high. Therefore twostep
procedure: (1) x2~c1+x1+u1 (2) y=c2+x1+u1+x3+u2.
res1<-lm(x2~x1)
coeftest(res1,vcov.=vcov(res1))
vcov(res1)
#I am able to calculate this manually.
X1<-as.matrix(data.frame(A=rep(1,100),B=x1),stringsAsFactors=FALSE)
betas<-solve(crossprod(X1,X1))%*%t(X1)%*%x2
round(coef(res1),digits=10)==round(betas,digits=10)
betas
u1<-resid(res1)
n1<-length(y)
k1<-ncol(X1)
vcov1 = 1/(n1-k1) * as.numeric(t(u1)%*%u1) * solve(t(X1)%*%X1) #= 1/T
vcov1
#Now I do my second regression:
Y<-as.matrix(y)
X2<-as.matrix(data.frame(A=1,B=x1,C=u1,D=x3))
res2<-lm(y~x1+u1+x3)
coeftest(res2,vcov.=vcov(res2))
vcov(res2)
#This is the (biased) variance-covariance matrix:
u2<-resid(res2)
n2<-length(y)
k2<-ncol(X2)
vcov2 = 1/(n2-k2) * as.numeric(t(u2)%*%u2) * solve(t(X2)%*%X2) #= 1/T
vcov2
#In order to implement the Shanken(1992) adjustment, I need to calculate the
variance-covariance matrix manually in the following way: (See Chochrane
(2005) Asset Pricing Ch.12):
var_coef<-function(x){
x_mm<-model.matrix(x)
x_s2<-summary(x)$sigma^2
solve(t(x_mm)%*%x_mm)%*%t(x_mm)%*%(x_s2*diag(nobs(x)))%*%x_mm%*%solve(t(x_mm)%*%x_mm)
}
vcov_manual<-var_coef(x=res1)
vcov1
round(vcov1,digits=10)==round(vcov_manual,digits=10)
#I calculate the covariance matrix of the residuals cov(u,u').
cov_resid<-function(x){
x_mm<-model.matrix(x)
x_s2<-summary(x)$sigma^2
(diag(nobs(x))-x_mm%*%solve(t(x_mm)%*%x_mm)%*%t(x_mm))%*%(x_s2*diag(nobs(x)))%*%t(diag(nobs(x))-x_mm%*%solve(t(x_mm)%*%x_mm)%*%t(x_mm))
}
cov_resid<-cov_resid(x=res1)
cov_resid[95:100,95:100] #to show only a part
#Another way to go to:
H1<-X1%*%solve(t(X1)%*%X1)%*%t(X1)
((diag(100)-H1)*summary(res1)$sigma^2)[95:100,95:100]
#Seems to be correct
round(((diag(100)-H1)*summary(res1)$sigma^2)[95:100,95:100],digits=10)==round(cov_resid[95:100,95:100],digits=10)
#Now I include the Shanken adjustment:
...*(1+t(lambdas)%*%solve(sig_f)%*%lambdas)...
var_coef_adj<-function(x_1,x_2){
lambdas<-coef(x_2)
x_mm<-model.matrix(x_2)
x_s<-((diag(100)-H1)*summary(x_1)$sigma^2)
x_s_f<-vcov(x_2)
1/1*(solve(t(x_mm)%*%x_mm)%*%t(x_mm)%*%(x_s)%*%x_mm%*%solve(t(x_mm)%*%x_mm)*as.numeric(1+t(lambdas)%*%solve(x_s_f)%*%lambdas)+x_s_f)
}
vcov2_adj<-var_coef_adj(x_1=res1, x_2=res2)
var_coef_adj(x_1=res1, x_2=res2)
coeftest(res2)
coeftest(res2,vcov.=vcov2_adj)
-----
____________________________________
EBS Universitaet fuer Wirtschaft und Recht
FARE Department
Wiesbaden/ Germany
http://www.ebs.edu/index.php?id=finacc&L=0
--
View this message in context:
http://r.789695.n4.nabble.com/Calculate-Adjusted-vcov-Matrix-acc-to-Shanken-1992-Generated-Regressor-Problem-tp4681404p4681631.html
Sent from the R help mailing list archive at Nabble.com.
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
