###############################################################################
# AEB 6572 - Econometric Methods II                                           #
# Charles B. Moss - cbmoss@ufl.edu                                            #
# R Regression Program GIVE03-19-12.R                                         #
###############################################################################
# This program demonstrates a system instrumental variable approach to        #
# motivate generalize method of moments estimation.                           #
###############################################################################
# Psuedo - data                                                               #
# dta[,1] - Observation Number                                                #
# dta[,2] - p2 (alternative output for supply equation)                       #
# dta[,3] - w1 (input price for supply equation)                              #
# dta[,4] - w2 (input price for supply equation)                              #
# dta[,5] - y (consumer income for demand equation)                           #
# dta[,6] - p3 (price of alternative good for demand equation)                #
# dta[,7] - p1 Endogenous price of good                                       #
# dta[,8] - q1 Endogenous quantity of good                                    #
###############################################################################

dta <- read.csv("3SLS03-13-12.csv")

###############################################################################
# Construct system x and y matrices                                           #
# x[,1] - beta_{10} Constant for supply equation                              #
# x[,2] - gamma_{11} Effect of price (p1)                                     #
# x[,3] - beta_{11} Effect of alternative output price (p2)                   #
# x[,4] - beta_{12} Effect of input price (w1)                                #
# x[,5] - beta_{13} Effect of input price (w2)                                #
# x[,6] - beta_{20} Constant for the demand equation                          #
# x[,7] - gamma_{12} Effect of price (p1)                                     #
# x[,8] - beta_{24} Effect of income (y)                                      #
# x[,9] - beta_{25} Effect of altenative consumption price (p3)               #
###############################################################################

x <- matrix(0,nrow=2*nrow(dta),ncol=9)
y <- matrix(0,nrow=2*nrow(dta),ncol=1)
for (i in 1:nrow(dta)) {
###############################################################################
# Data for the supply equation - 2*(i-1)+1                                    #
###############################################################################
  x[2*(i-1)+1,1] <- 1
  x[2*(i-1)+1,2] <- dta[i,7]
  x[2*(i-1)+1,3] <- dta[i,2]
  x[2*(i-1)+1,4] <- dta[i,3]
  x[2*(i-1)+1,5] <- dta[i,4]
###############################################################################
# Data for the demand equation - 2*(i-1)+2                                    #
###############################################################################
  x[2*(i-1)+2,6] <- 1
  x[2*(i-1)+2,7] <- dta[i,7]
  x[2*(i-1)+2,8] <- dta[i,5]
  x[2*(i-1)+2,9] <- dta[i,6]
###############################################################################
# Data for dependent variable                                                 #
###############################################################################
  y[2*(i-1)+1] <- dta[i,8]
  y[2*(i-1)+2] <- dta[i,8]
}

###############################################################################
# Instrumental Variables                                                      #
# z[,1]  - constant                                                           #
# z[,2]  - p2 (alternative output for supply equation)                        #
# z[,3]  - w1 (input price for supply equation)                               #
# z[,4]  - w2 (input price for supply equation)                               #
# z[,5]  - y (consumer income for demand equation)                            #
# z[,6]  - p3 (price of alternative good for demand equation)                 #
###############################################################################
# For the equivalent to the 2SLS estimator, we know the W_T matrix            #
###############################################################################

z <- cbind(matrix(1,nrow=nrow(dta),ncol=1),as.matrix(dta[,2:6]))
zvar <- t(z)%*%z/nrow(z)

###############################################################################
# For 2SLS I set the variance matrix to an identity times the variance of the #
#  instrument variance                                                        #
###############################################################################

what <- diag(1,nrow=2,ncol=2)%x%solve(zvar)
zz <- z %x% diag(1,nrow=2,ncol=2)
b2sls <- solve(t(x)%*%zz%*%what%*%t(zz)%*%x)%*%(t(x)%*%zz%*%what%*%t(zz)%*%y)
print(b2sls)

###############################################################################
# Next, I can estimate the variance for each equation                         #
###############################################################################

errhat <- y - x%*%b2sls
err    <- matrix(0,nrow=nrow(dta),ncol=2)
for (i in 1:nrow(dta)) { err[i,] <- cbind(errhat[2*(i-1)+1],errhat[2*(i-1)+2])}
vhat <- t(err)%*%err/nrow(err)

print(vhat)
what <- solve(vhat)%x%solve(zvar)
b3sls <- solve(t(x)%*%zz%*%what%*%t(zz)%*%x)%*%(t(x)%*%zz%*%what%*%t(zz)%*%y)

print(cbind(b2sls,b3sls))

###############################################################################
# Now, let's extend the estimator by adding some more instruments             #
###############################################################################

z      <- cbind(z,z[,2]*z[,2],z[,2]*z[,3],z[,2]*z[,4],
                z[,3]*z[,3],z[,4]*z[,4],z[,5]*z[,5])
zvar   <- t(z)%*%z/nrow(z)
what   <- diag(1,nrow=2,ncol=2)%x%solve(zvar)
zz     <- z %x% diag(1,nrow=2,ncol=2)
bgmm1  <- solve(t(x)%*%zz%*%what%*%t(zz)%*%x)%*%(t(x)%*%zz%*%what%*%t(zz)%*%y)

errhat <- y - x%*%bgmm1
err    <- matrix(0,nrow=nrow(dta),ncol=2)
for (i in 1:nrow(dta)) { err[i,] <- cbind(errhat[2*(i-1)+1],errhat[2*(i-1)+2])}
vhat   <- t(err)%*%err/nrow(err)
what   <- solve(vhat)%x%solve(zvar)
bgmm2  <- solve(t(x)%*%zz%*%what%*%t(zz)%*%x)%*%(t(x)%*%zz%*%what%*%t(zz)%*%y)


print(cbind(b2sls,b3sls,bgmm1,bgmm2))