nba1415 <- read.csv("http://www.stat.ufl.edu/~winner/data/nbaodds201415.csv", header=T) attach(nba1415); names(nba1415) (N.pop <- length(TotalPts[OT==0])) Y <- TotalPts[OT==0]/100 X0 <- rep(1,N.pop) X1a <- OvrUndr[OT==0]/100 X1 <- X1a - mean(X1a) X <- cbind(X0,X1) (beta <- solve(t(X)%*%X) %*% (t(X)%*%Y)) E.Y <- X %*% beta eps <- Y - E.Y (sigma2 <- sum(eps^2)/N.pop) n.sample <- 25 (M_XX <- (1/N.pop) * (t(X)%*%X)) epsepst <- diag(diag(eps %*% t(eps))) (M_XOmX <- (1/N.pop) * (t(X)%*% epsepst %*% X)) (V.n12.bols.beta <- solve(M_XX) %*% M_XOmX %*% solve(M_XX)) (V.bols.n <- V.n12.bols.beta/n.sample) par(mfrow=c(2,1)) plot(X1a,Y,main="Total Points (Y) vs Over/Under (X) - NBA 2014/15", xlab="Over/Under",ylab="Total Points") abline(lm(Y ~ X1a)) plot(E.Y,eps,main="Errors versus Expected Value") hist(X1a,breaks=30,main="Over/Under") hist(eps,breaks=30,main=expression(paste(epsilon))) mean(X1); sd(X1); mean(eps); sd(eps); cor(X1,eps) par(mfrow=c(1,1)) set.seed(135678) n.sim <- 100000 beta.hat <- matrix(rep(0,2*n.sim),ncol=2) s2 <- numeric(n.sim) xpx <- matrix(rep(0,3*n.sim),ncol=3) xpe <- matrix(rep(0,2*n.sim),ncol=2) for (i in 1:n.sim) { nba.sample <- sample(1:N.pop,n.sample,replace=F) X.s <- X[nba.sample,] Y.s <- Y[nba.sample] eps.s <- eps[nba.sample] XPX.s <- t(X.s) %*% X.s beta.hat.s <- solve(XPX.s) %*% t(X.s) %*% Y.s beta.hat[i,] <- t(beta.hat.s) e.s <- Y.s - X.s%*%beta.hat.s s2[i] <- (t(e.s) %*% (e.s))/(n.sample-2) XPe.s <- t(X.s) %*% eps.s xpx[i,1] <- XPX.s[1,1] xpx[i,2] <- XPX.s[1,2] xpx[i,3] <- XPX.s[2,2] xpe[i,1] <- XPe.s[1,1] xpe[i,2] <- XPe.s[2,1] } mean(s2) (E.xpx <- matrix(c(mean(xpx[,1]),mean(xpx[,2]),mean(xpx[,2]),mean(xpx[,3])), byrow=T,ncol=2)) mean(beta.hat[,1]); mean(beta.hat[,2]) var(beta.hat[,1]); var(beta.hat[,2]); cov(beta.hat[,1],beta.hat[,2]) cor(beta.hat[,1],beta.hat[,2]) (V.beta.hat <- sigma2 * solve(E.xpx)) xpe <-(1/sqrt(n.sample))*xpe (mean(xpe[,1])) (mean(xpe[,2])) (var(xpe[,1])) (var(xpe[,2])) (cov(xpe[,1],xpe[,2])) par(mfrow=c(2,1)) hist(beta.hat[,1],breaks=50,main="Histogram of OLS Estmates of B0") hist(beta.hat[,2],breaks=50,main="Histogram of OLS Estmates of B1") par(mfrow=c(1,1)) plot(beta.hat[,1],beta.hat[,2],pch=16,cex=.4, main="Scatterplot of Estimates of B1 vs B0",xlab="B0",ylab="B1")