############################ # # Lecture 4: Statistical Properties of OLS # POLS 509: The Linear Model # Dr. Justin Esarey # ############################ # set random seed set.seed(123456) # what does the bias in ADL models look like? # do a simulation to see b.hat.store<-matrix(data=NA, ncol=2, nrow=5000) colnames(b.hat.store)<-c("constant", "y lag") pb<-txtProgressBar(min=0, max=5000, style=3) for(j in 1:5000){ setTxtProgressBar(pb, value=j) # create a simple ADL model y<-runif(1, min=0, max=10) for(i in 1:30){ y[i+1]<-.435+.859*y[i]+rnorm(1,mean=0, sd=1) } # run a regression of y_t on y_t-1 y.std<-y[2:31] y.lag<-y[1:30] X<-cbind(1,y.lag) b.hat<-qr.solve(t(X)%*%X)%*%t(X)%*%y.std b.hat.store[j,]<-t(b.hat) } mean(b.hat.store[,1])-.435 mean(b.hat.store[,2])-.859 plot(density(b.hat.store[,1]), main="Bias in Est. Constant for ADL Model") abline(v=.435, lty=2) # this is the true value abline(v=mean(b.hat.store[,1]), lty=2, lwd=3) legend("topright", lty=c(2,2), lwd=c(1,3), legend=c("True Value", "Mean Estimate")) plot(density(b.hat.store[,2]), main="Bias in Est. Lag Coef. for ADL Model") abline(v=.859, lty=2) # this is the true value abline(v=mean(b.hat.store[,2]), lty=2, lwd=3) legend("topleft", lty=c(2,2), lwd=c(1,3), legend=c("True Value", "Mean Estimate")) # example: find a VCV matrix # create some data X<-matrix(data=runif(100, min=0, max=10), ncol=2) X<-cbind(1,X) b<-c(2.8, 1.3, 6.5) y<-X%*%b+rnorm(50, mean=0, sd=2) # run a regression b.hat<-qr.solve(t(X)%*%X)%*%t(X)%*%y y.hat<-X%*%b.hat u<-y-X%*%b.hat plot(y~y.hat) abline(c(0,1),lty=1) # get VCV matrix uTu<-(t(u)%*%u) vcv<-as.vector(uTu)*(1/(50-3))*qr.solve(t(X)%*%X) summary(lm(y~X-1)) sqrt(diag(vcv)) vcv vcov(lm(y~X-1))