##########################################################################
# draw[,1]  - Price of good 1                                            #
# draw[,2]  - Quantity of good 1                                         #
# draw[,3]  - Price of good 2                                            #
# draw[,4]  - Income                                                     #
##########################################################################

vec <- function(z) { 
   dimz <- dim(z)
   zret <- t(z[1,])
   for ( i in 2:dimz[1] ) {
     zret <- cbind(zret,t(z[i,])) 
     }
   t(zret)
}

y     <- vec(draw[,1:2])
x     <- cbind(matrix(1,nrow=nrow(draw),ncol=1),draw[,3],draw[,4])%x%diag(1,2,2)
biols <- solve(t(x)%*%x)%*%(t(x)%*%y)
 
erhat <- y-x%*%biols

errhat <- matrix(0,nrow=nrow(draw),ncol=2)
for (i in 1:nrow(draw)) {
  errhat[i,1] <- erhat[2*(i-1)+1,]
  errhat[i,2] <- erhat[2*(i-1)+2,] }

print(errhat)