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()