bn <-c("CRIM","ZN","INDUS","CHAS","NOX","RM","AGE","DIS","RAD","TAX",
        "PTRATIO","B","LSTAT","MV")
bh <- read.table("bh.dat",col.names=bn)

n <- nrow(bh)

ones <- mat.or.vec(n,1) + 1 
ones <- as.data.frame(ones)
colnames(ones) <- "(Intercept)"
bh <- cbind(bh,ones)

bh["MV"] <- 1000.0*bh["MV"]
bh["B"] <- bh["B"]/1000.0
bh["NOX"] <- 10.0*bh["NOX"]
bh["LSTAT"] <- bh["LSTAT"]/100.0

#b0 <- sqrt(bh["B"]) + 0.63 
#colnames(b0) <- "B0"
#bh <- cbind(bh,b0)

nox2 <- bh["NOX"]^2
colnames(nox2) <- "NOX2"
bh <- cbind(bh,nox2)

rm2 <- bh["RM"]^2
colnames(rm2) <- "RM2"
bh <- cbind(bh,rm2)

logMV <- log(bh["MV"])
colnames(logMV) <- "logMV"
bh <- cbind(bh,logMV)

logDIS <- log(bh["DIS"])
colnames(logDIS) <- "logDIS"
bh <- cbind(bh,logDIS)

logRAD <- log(bh["RAD"])
colnames(logRAD) <- "logRAD"
bh <- cbind(bh,logRAD)

logLSTAT <- log(bh["LSTAT"])
colnames(logLSTAT) <- "logLSTAT"
bh <- cbind(bh,logLSTAT)

write.table(bh,file="bh_reg1.csv",sep=",",row.names=F,col.names=T)
write.table(bh,file="bh_reg1.dat",sep=" ",row.names=F,col.names=F)

hedonic <- lm("logMV ~ RM2+AGE+logDIS+logRAD+TAX+PTRATIO+B+logLSTAT+CRIM+ZN+INDUS+CHAS+NOX2",bh) 

print("mean(bh) = ")
print(mean(bh))

print("sd(bh) = ")
print(sd(bh))

print("hedonic = ")
print(hedonic)
summary(hedonic)

b <- hedonic$coefficients

logMVhat <- b["(Intercept)"]*bh["(Intercept)"] 
logMVhat <- logMVhat + b["RM2"]*bh["RM2"] 
logMVhat <- logMVhat + b["AGE"]*bh["AGE"] 
logMVhat <- logMVhat + b["logDIS"]*bh["logDIS"] 
logMVhat <- logMVhat + b["logRAD"]*bh["logRAD"] 
logMVhat <- logMVhat + b["TAX"]*bh["TAX"] 
logMVhat <- logMVhat + b["PTRATIO"]*bh["PTRATIO"] 
logMVhat <- logMVhat + b["B"]*bh["B"] 
logMVhat <- logMVhat + b["logLSTAT"]*bh["logLSTAT"] 
logMVhat <- logMVhat + b["CRIM"]*bh["CRIM"] 
logMVhat <- logMVhat + b["ZN"]*bh["ZN"] 
logMVhat <- logMVhat + b["INDUS"]*bh["INDUS"] 
logMVhat <- logMVhat + b["CHAS"]*bh["CHAS"] 
logMVhat <- logMVhat + b["NOX2"]*bh["NOX2"] 

colnames(logMVhat) <- "logMVhat"
bh <- cbind(bh,logMVhat)

MVhat <- exp(logMVhat)
colnames(MVhat) <- "MVhat"
bh <- cbind(bh,MVhat)

What <- -MVhat*(2*b["NOX2"]*bh["NOX"])
colnames(What) <- "What"
bh <- cbind(bh,What)

logWhat <- log(What)
colnames(logWhat) <- "logWhat"
bh <- cbind(bh,logWhat)

logNOX <- log(bh["NOX"])
colnames(logNOX) <- "logNOX"
bh <- cbind(bh,logNOX)

ps.options(horizontal=FALSE,onefile=FALSE);
ps.options(paper="special",pagecentre=FALSE);
ps.options(width=8.0,height=4.0);

postscript(file="bh01.eps")
par(mfrow=c(1,2),mar=c(4,4,1.5,2)+0.1) # mar=c(bottom,left,top,right)- 

idx17000 <- (bh[,"MV"] <= 17000)
idx23000 <- (bh[,"MV"] > 17000) & (bh[,"MV"] <= 23000)
idx30000 <- (bh[,"MV"] >= 30000)

plot(x=bh[,"NOX"],y=bh[,"What"],type="n",
  ylab="Marginal Willingness to Pay, $",
  xlab="Nitrous Oxide, pphm")
  points(x=bh[idx17000,"NOX"],y=bh[idx17000,"What"],pch='o',col="green")
  points(x=bh[idx23000,"NOX"],y=bh[idx23000,"What"],pch='o',col="green")
  points(x=bh[idx30000,"NOX"],y=bh[idx30000,"What"],pch='o',col="red")

plot(x=bh[,"MV"],y=bh[,"What"],type="n",
  ylab="Marginal Willingness to Pay, $",
  xlab="Median House Value, $")
  points(x=bh[,"MV"],y=bh[,"What"],pch='o')

dev.off()

write.table(bh,file="bh_reg2.csv",sep=",",row.names=F,col.names=T)
write.table(bh,file="bh_reg2.dat",sep=" ",row.names=F,col.names=F)

#demand <- lm("What ~ NOX + MV",bh) 
demand <- lm("logWhat ~ logNOX + logMV",bh) 

print("demand =")
print(demand)
summary(demand)

c <- demand$coefficients

wtp <- seq(1,8.5,0.25)
m <- length(wtp)
wtp <- as.data.frame(wtp)
colnames(wtp) <- "NOX"
logwtp <- log(wtp)
colnames(logwtp) <- "logNOX"
wtp <- cbind(wtp,logwtp)

ones <- mat.or.vec(m,1) + 1 
ones <- as.data.frame(ones)
colnames(ones) <- "(Intercept)"
wtp <- cbind(wtp,ones)

inc17000  <- log(17000.0)*ones
inc23000 <- log(23000.0)*ones
inc30000 <- log(30000.0)*ones

logWhat <- c["(Intercept)"]*wtp["(Intercept)"] 
logWhat <- logWhat + c["logNOX"]*wtp["logNOX"] 

logW17000 <- logWhat + c["logMV"]*inc17000
logW23000 <- logWhat + c["logMV"]*inc23000
logW30000 <- logWhat + c["logMV"]*inc30000

W17000 <- exp(logW17000)
W23000 <- exp(logW23000)
W30000 <- exp(logW30000)

colnames(W17000) <- "W17000"
colnames(W23000) <- "W23000"
colnames(W30000) <- "W30000"

wtp <- cbind(wtp,W17000,W23000,W30000)

ps.options(horizontal=FALSE,onefile=FALSE);
ps.options(paper="special",pagecentre=FALSE);
ps.options(width=8.0,height=5.0);

postscript(file="bh02.eps")
par(mfrow=c(1,1),mar=c(4,4,1.5,2)+0.1) # mar=c(bottom,left,top,right)

xlim <- c(1,10)
ylim <- c(min(wtp["W17000"]),max(wtp["W30000"]))

plot(x=xlim,y=ylim,type="n",
  ylab="Marginal Willingness to Pay, $",
  xlab="Nitrous Oxide, pphm")
  lines(x=wtp[,"NOX"],y=wtp[,"W17000"],lty="solid",col="blue")
  text(x=c(wtp[m,"NOX"]+0.75),y=c(wtp[m,"W17000"]),label="MV=17,000")
  lines(x=wtp[,"NOX"],y=wtp[,"W23000"],lty="solid",col="green")
  text(x=c(wtp[m,"NOX"]+0.75),y=c(wtp[m,"W23000"]),label="MV=23,000")
  lines(x=wtp[,"NOX"],y=wtp[,"W30000"],lty="solid",col="red")
  text(x=c(wtp[m,"NOX"]+0.75),y=c(wtp[m,"W30000"]),label="MV=30,000")

dev.off()