pdf("nonnorm_e.pdf") beta <- matrix(c(50,10),ncol=1) X0 <- matrix(rep(1,22),ncol=1) X1 <- matrix(rep(seq(0,10,1),2),ncol=1) X <- cbind(X0,X1) n <- nrow(X) EY <- X %*% beta XX <- t(X)%*%X XXI <- solve(t(X)%*%X) XXIX <- solve(t(X)%*%X) %*% t(X) XX XXI nsim <- 10000 tn_2 <- qt(.975,n-2) #### Non-normal Errors (3*t(3)) #### betahat <- matrix(rep(0,(2*nsim)),ncol=2) s2 <- matrix(rep(0,nsim),ncol=1) s2beta <- matrix(rep(0,(4*nsim)),ncol=4) beta1lo <- matrix(rep(0,nsim),ncol=1) beta1hi <- matrix(rep(0,nsim),ncol=1) covrate <- 0 for (i in 1:nsim) { eps <- 3.0*rt(n,3) Y <- EY + eps bsim <- XXIX %*% Y s2sim <- (t(Y-X%*%bsim) %*% (Y-X%*%bsim))/(n-2) s2bsim <- s2sim[1,1]* XXI beta1lo[i] <- bsim[2,1] - tn_2*sqrt(s2bsim[2,2]) beta1hi[i] <- bsim[2,1] + tn_2*sqrt(s2bsim[2,2]) if ((beta1lo[i] <= beta[2,1]) & (beta1hi[i] >= beta[2,1])) covrate=covrate+1 betahat[i,] <- t(bsim) s2[i] <- s2sim for (i1 in 1:2) { for (i2 in 1:2) { s2beta[i,(((i1-1)*2)+i2)] <- s2bsim[i1,i2] }} } mean(betahat[,2]) sd(betahat[,2]) min(betahat[,2]) max(betahat[,2]) mean(s2) var(s2) mean(s2beta[,4]) mean(beta1hi-beta1lo) covrate <- covrate/nsim covrate brk <- seq(-10.0,30.0,0.01) h1 <- hist(betahat[,2], plot=F) range(0, h1\$density, dnorm(10,10,.35032)) ylim <- range(0, h1\$density, dnorm(10,10,.35032)) hist(betahat[,2], breaks=brk, freq=F, xlim=c(9.0,11.0), ylim=c(0,1.5)) curve(dnorm(x,10,0.35032), add=T) ### Heteroskedastic Errors e ~ N(0,((X+.5))^2) #### betahat <- matrix(rep(0,(2*nsim)),ncol=2) s2 <- matrix(rep(0,nsim),ncol=1) s2beta <- matrix(rep(0,(4*nsim)),ncol=4) beta1lo <- matrix(rep(0,nsim),ncol=1) beta1hi <- matrix(rep(0,nsim),ncol=1) covrate <- 0 for (i in 1:nsim) { eps <- rnorm(22) Y <- EY + eps*(X[,2]+0.5) bsim <- XXIX %*% Y s2sim <- (t(Y-X%*%bsim) %*% (Y-X%*%bsim))/20 s2bsim <- s2sim[1,1]* XXI beta1lo[i] <- bsim[2,1] - tn_2*sqrt(s2bsim[2,2]) beta1hi[i] <- bsim[2,1] + tn_2*sqrt(s2bsim[2,2]) if ((beta1lo[i] <= beta[2,1]) & (beta1hi[i] >= beta[2,1])) covrate=covrate+1 betahat[i,] <- t(bsim) s2[i] <- s2sim for (i1 in 1:2) { for (i2 in 1:2) { s2beta[i,(((i1-1)*2)+i2)] <- s2bsim[i1,i2] }} } mean(betahat[,2]) sd(betahat[,2]) mean(s2) var(s2) mean(s2beta[,4]) mean(beta1hi-beta1lo) covrate <- covrate/nsim covrate brk <- seq(-10.0,30.0,0.01) h1 <- hist(betahat[,2], plot=F) range(0, h1\$density, dnorm(10,10,0.4763)) ylim <- range(0, h1\$density, dnorm(10,10,.4763)) hist(betahat[,2], breaks=brk, freq=F, xlim=c(8.0,12.0), ylim=c(0,1.0)) curve(dnorm(x,10,0.4763), add=T) ### Correlated Errors e_t = 0.5*e_(t-1) + v_t #### betahat <- matrix(rep(0,(2*nsim)),ncol=2) s2 <- matrix(rep(0,nsim),ncol=1) s2beta <- matrix(rep(0,(4*nsim)),ncol=4) beta1lo <- matrix(rep(0,nsim),ncol=1) beta1hi <- matrix(rep(0,nsim),ncol=1) covrate <- 0 sigma <- 3 rho <- 0.5 for (i in 1:nsim) { eps <- rep(0,n) v <- rnorm(n-1,mean=0,sd=sigma) eps[1] <- rnorm(1,mean=0,sd=(sigma/sqrt(1-rho^2))) for (i1 in 2:n) { eps[i1] <- rho*eps[i1-1] + v[i1-1] } Y <- EY + eps bsim <- XXIX %*% Y s2sim <- (t(Y-X%*%bsim) %*% (Y-X%*%bsim))/20 s2bsim <- s2sim[1,1]* XXI beta1lo[i] <- bsim[2,1] - tn_2*sqrt(s2bsim[2,2]) beta1hi[i] <- bsim[2,1] + tn_2*sqrt(s2bsim[2,2]) if ((beta1lo[i] <= beta[2,1]) & (beta1hi[i] >= beta[2,1])) covrate=covrate+1 betahat[i,] <- t(bsim) s2[i] <- s2sim for (i1 in 1:2) { for (i2 in 1:2) { s2beta[i,(((i1-1)*2)+i2)] <- s2bsim[i1,i2] }} } mean(betahat[,2]) sd(betahat[,2]) mean(s2) var(s2) mean(s2beta[,4]) mean(beta1hi-beta1lo) covrate <- covrate/nsim covrate brk <- seq(-10.0,30.0,0.01) h1 <- hist(betahat[,2], plot=F) range(0, h1\$density, dnorm(10,10,.035032)) ylim <- range(0, h1\$density, dnorm(10,10,.2629)) hist(betahat[,2], breaks=brk, freq=F, xlim=c(9.0,11.0), ylim=c(0,1.5)) curve(dnorm(x,10,0.2629), add=T) ### Measurement Error- Observed X, True value=Z="Old" X (0,1,...,10)x2, X=Z+U U~N(0,0.1667^2) ### betahat <- matrix(rep(0,(2*nsim)),ncol=2) s2 <- matrix(rep(0,nsim),ncol=1) s2beta <- matrix(rep(0,(4*nsim)),ncol=4) beta1lo <- matrix(rep(0,nsim),ncol=1) beta1hi <- matrix(rep(0,nsim),ncol=1) covrate <- 0 sigma_u <- 1 sigma_e <- 3 XE0 <- rep(1,22) for (i in 1:nsim) { U <- rnorm(22) XE1 <- X[,2]+sigma_u*U XE <- cbind(XE0,XE1) eps <- rnorm(22) Y <- EY + sigma_e*eps XX <- t(XE)%*%XE XXI <- solve(t(XE)%*%XE) XXIX <- solve(t(XE)%*%XE) %*% t(XE) bsim <- XXIX %*% Y s2sim <- (t(Y-XE%*%bsim) %*% (Y-XE%*%bsim))/20 s2bsim <- s2sim[1,1]* XXI beta1lo[i] <- bsim[2,1] - tn_2*sqrt(s2bsim[2,2]) beta1hi[i] <- bsim[2,1] + tn_2*sqrt(s2bsim[2,2]) if ((beta1lo[i] <= beta[2,1]) & (beta1hi[i] >= beta[2,1])) covrate=covrate+1 betahat[i,] <- t(bsim) s2[i] <- s2sim for (i1 in 1:2) { for (i2 in 1:2) { s2beta[i,(((i1-1)*2)+i2)] <- s2bsim[i1,i2] }} } mean(betahat[,2]) sd(betahat[,2]) mean(s2) var(s2) mean(s2beta[,4]) mean(beta1hi-beta1lo) covrate <- covrate/nsim covrate brk <- seq(-10.0,30.0,0.01) h1 <- hist(betahat[,2], plot=F) range(0, h1\$density, dnorm(10,10,0.20226)) ylim <- range(0, h1\$density, dnorm(10,10,.20226)) hist(betahat[,2], breaks=brk, freq=F, xlim=c(7.0,13.0), ylim=c(0,2.0)) curve(dnorm(x,10,0.20226), add=T) dev.off()