# dta <- read.table("PseudoLand.dta")
# dta <- read.table("AlabamaLand.dta")
# dta <- read.table("ArkansasLand.dta")
# dta <- read.table("FloridaLand.dta")
# dta <- read.table("GeorgiaLand.dta")
# dta <- read.table("KentuckyLand.dta")
# dta <- read.table("Louisiana.dta")
# dta <- read.table("MarylandLand.dta")
dta[,1:3] <- dta[,1:3]/100
dv <- as.matrix(dta[,1])
x  <- as.matrix(cbind(dta[,2:4]))
tt <- c(1:nrow(x))
z  <- as.matrix(cbind(dta[,2:4],dta[,2:3]*dta[,2:3],tt,1))
z[,4] <- z[,4]/10
# z  <- as.matrix(cbind(dta[,2:4],dta[,2:4]*dta[,2:4]))
# z  <- as.matrix(dta[,2:7])
n  <- nrow(dta)
m  <- ncol(z)
w  <- diag(x=1,nrow=m,ncol=m)
w  <- solve(t(z) %*% z/nrow(z))

df <- function(b) {
      bb <- (1+b[4]*x[,3])*(1+b[4]*x[,3])
      g1 <- -1.0
      g2 <- -b[4]*x[,3]*x[,1]/(1+b[4]*x[,3])
      g3 <- x[,2]/(1+b[4]*x[,3])
      g4 <- -b[2]*x[,3]*x[,1]/(1+b[4]*x[,3]) + b[2]*b[4]*x[,3]*x[,3]*x[,1]/bb - b[3]*x[,3]*x[,2]/bb
      as.matrix(cbind(g1,g2,g3,g4))}

landval <- function(b) {
          as.matrix(dv-b[1]-b[2]*(b[4]*x[,3])*x[,1]/(1+(b[4]*x[,3]))+b[3]*x[,2]/(1+(b[4]*x[,3]))) }

lsq <- function(b) {
        as.numeric(t(landval(b)) %*% landval(b) )}

fr <- function(b) {
        landv <- landval(b)
        t(landv) %*% z %*% w %*% t(z) %*% landv/n }

b0 <- cbind(1.0,1.0,1.0,0.008)
res.lsq <- optim(b0,lsq)
print("Least Squares")
print(res.lsq$par)

print("Original GMM")
res.gmm <- optim(b0,fr)
print(res.gmm$par)
g    <- z * (landval(res.gmm$par) %*% matrix(1,nrow=1,ncol=m))
shat <- t(g) %*% g/nrow(z)

for (i in 1:20) {
    w <- solve(shat)
    res.gmm <- optim(b0,fr)
    g     <- z * (landval(res.gmm$par) %*% matrix(1,nrow=1,ncol=m))
    shato <- shat
    shat  <- t(g) %*% g/nrow(z) }

print("Final GMM")
print(res.gmm$par)
print(shat-shato)
print(sqrt(sum((shat-shato)*(shat-shato))))

gg   <- t(z) %*% df(res.gmm$par)/nrow(z)
mt   <- solve(t(gg) %*% w %*% gg) %*% t(gg) %*% w
vt   <- mt %*% shat %*% t(mt)/nrow(z)
print(cbind(t(res.gmm$par),sqrt(diag(vt)),t(res.gmm$par)/sqrt(diag(vt))))