rm(list = ls())

datamat <- matrix(scan("hrdata.txt"),byrow=T,ncol=15)

y <- datamat[,3]
x <- datamat[,7:15]

# fields are 
# 1) idorigin 2) year 3) ainew 4) sdnew 5) ailag 6) sdlag 7) polrt 
# 8) lpop 9) perchpop 10) pcgnp 11) perchpcg 12) left 13) mil2 
# 14) iwarcow2 15) cwarcow;

n <- 110
t <- 17
nt <- n*t

# define useful matrices 

I.t <- diag(1,t,t)

e <- matrix(1,t,1)

J.t <- e %*% t(e)

J.t <- matrix(1,t,t)

qmat <- I.t - 1/t*J.t

I.nt <- diag(1,n*t,n*t) 

I.n <- diag(1,n,n)

qmatall <- I.nt - I.n %x% (1/t*J.t)

xc <- cbind(rep(1,nrow(x)),x)

k <- ncol(x)

inxx.w <- solve(t(x) %*% qmatall %*% x)

b.w <- inxx.w %*% t(x) %*% qmatall %*% y

# FE Standard Errors

y.tilde <- as.matrix(qmatall %*% y)
x.tilde <- as.matrix(qmatall %*% x)

u.tilde <- as.matrix(y.tilde - (x.tilde %*% b.w))
ssq.w <- (t(u.tilde) %*% u.tilde)/(nt-n-k)
se.w <- sqrt(diag(as.numeric(ssq.w) * inxx.w))

fevarnames <- as.matrix(c("polrt", "lpop", "perchpop", "pcgnp", "perchpcg", "left", "mil2", "iwarcow2", "cwarcow"))

print("FE estimates")
print(cbind(fevarnames,b.w,se.w),digits=3)