rm(list = ls())

n <- 100
t <- 5
nt <- n*t
k <- 1

sim <- 1000

beta <- 1
rho <- .8 # correlation b/t alpha and x
amean <- 0

I.t <- diag(1,t,t)
J.t <- matrix(1,t,t)
qmat <- I.t - 1/t*J.t
I.nt <- diag(1,nt,nt) 
I.n <- diag(1,n,n)
qmatall <- I.nt - I.n %x% (1/t*J.t)
tmeansmat1 <- I.n %x% (1/t*matrix(1,1,t)) # matrix to get N vector of time means
b.wvec <- rep(0,sim)
b.glsvec <- rep(0,sim)

for(i in 1:sim){

 xraw <- as.matrix(rnorm(nt,k))
 a.i <- rnorm(n,amean,1)
 araw <- kronecker(a.i,matrix(1,t,1))
 cormat <- (1-rho)*diag(1,(k+1),(k+1))+rho*matrix(1,(k+1),(k+1))
 cormat.chold <- chol(cormat)
 axraw <- cbind(araw,xraw)
 ax <- axraw%*%cormat.chold
 a <- as.matrix(ax[,1])
 x <- as.matrix(ax[,2:(1+ncol(xraw))])
 
 cor.ax <- cor(x, a)
 cat("\nCorrelation b/t alpha and x:", cor.ax, "\n")
 
 y <- beta * x + a + rnorm(nt)
 
# Fixed effects
 
 inxx.w <- solve(t(x) %*% qmatall %*% x)
 b.w <- inxx.w %*% t(x) %*% qmatall %*%y
 y.tilde <- as.matrix(qmatall %*% y)
 x.tilde <- as.matrix(qmatall %*% x)
 u.tilde <- as.matrix(y.tilde - (x.tilde %*% b.w))
 sigmasq.u <- (t(u.tilde) %*% u.tilde)/(n*(t-1)-k)
 se.w <- sqrt(diag(as.numeric(sigmasq.u) * inxx.w))

 b.wvec[i] <- b.w
  
 cat("iteration:",i,"\n")
 
}

cat("mean of Within beta:", mean(b.wvec), "\n")