dta <- read.table("AggregAg.dta")
x <- as.matrix(dta[,5:9])
z <- as.matrix(cbind(dta[,5:9],dta[,5:9]*dta[,5:9]))
y <- as.matrix(dta[,1])

w <- diag(x=1,nrow=10,ncol=10)

fr <- function(b){
       t(y-x%*%b)%*%z%*%w%*%(t(z)%*%(y-x%*%b)) }

gfr <- function(b){
      -2*t(y-x%*%b)%*%z%*%w%*%t(z)%*%x }

res<-optim(c(1.0,1.0,1.0,1.0,1.0),fr,gfr,method="BFGS")

sse <- as.numeric(t(y-x%*%res$par)%*%(y-x%*%res$par)/nrow(y))
shat <- (sse/nrow(z))*(t(z)%*%z)
mhat <- nrow(z)*solve((t(x)%*%z)%*%w%*%(t(z)%*%x))%*%(t(x)%*%z)%*%w
vhat <- mhat %*% shat %*% t(mhat)

print(cbind(res$par,sqrt(diag(vhat)/nrow(x)),res$par/sqrt(diag(vhat)/nrow(x)),
            1-pnorm(abs(res$par/sqrt(diag(vhat)/nrow(x))))))

lm.rhs1 <- lm(y ~ 0 + x)
print(summary.lm(lm.rhs1))

w <- solve(shat)
res<-optim(c(1.0,1.0,1.0,1.0,1.0),fr,gfr,method="BFGS")
sse <- as.numeric(t(y-x%*%res$par)%*%(y-x%*%res$par)/nrow(y))
shat <- (sse/nrow(z))*(t(z)%*%z)
mhat <- nrow(z)*solve((t(x)%*%z)%*%w%*%(t(z)%*%x))%*%(t(x)%*%z)%*%w
vhat <- mhat %*% shat %*% t(mhat)

print(cbind(res$par,sqrt(diag(vhat)/nrow(x)),res$par/sqrt(diag(vhat)/nrow(x)),
            1-pnorm(abs(res$par/sqrt(diag(vhat)/nrow(x))))))

jt <- 1/nrow(z)*t(y-x%*%res$par)%*%z%*%solve(shat)%*%t(z)%*%(y-x%*%res$par)
print(jt)
print(pchisq(jt,5))