R version 2.11.1 (2010-05-31)
Copyright (C) 2010 The R Foundation for Statistical Computing
ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> 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) = ")
[1] "mean(bh) = "
> print(mean(bh))
         CRIM            ZN         INDUS          CHAS           NOX 
 3.613524e+00  1.136364e+01  1.113678e+01  6.916996e-02  5.546951e+00 
           RM           AGE           DIS           RAD           TAX 
 6.284634e+00  6.857490e+01  3.795043e+00  9.549407e+00  4.082372e+02 
      PTRATIO             B         LSTAT            MV   (Intercept) 
 1.845553e+01  3.566740e-01  1.265306e-01  2.253281e+04  1.000000e+00 
         NOX2           RM2         logMV        logDIS        logRAD 
 3.210877e+01  3.998932e+01  9.942268e+00  1.188032e+00  1.867661e+00 
     logLSTAT 
-2.234205e+00 
> 
> print("sd(bh) = ")
[1] "sd(bh) = "
> print(sd(bh))
        CRIM           ZN        INDUS         CHAS          NOX           RM 
8.601545e+00 2.332245e+01 6.860353e+00 2.539940e-01 1.158777e+00 7.026171e-01 
         AGE          DIS          RAD          TAX      PTRATIO            B 
2.814886e+01 2.105710e+00 8.707259e+00 1.685371e+02 2.164946e+00 9.129486e-02 
       LSTAT           MV  (Intercept)         NOX2          RM2        logMV 
7.141062e-02 9.197104e+03 0.000000e+00 1.392125e+01 9.079531e+00 4.087569e-01 
      logDIS       logRAD     logLSTAT 
5.395465e-01 8.748333e-01 6.008913e-01 
> 
> print("hedonic = ")
[1] "hedonic = "
> print(hedonic)

Call:
lm(formula = "logMV ~ RM2+AGE+logDIS+logRAD+TAX+PTRATIO+B+logLSTAT+CRIM+ZN+INDUS+CHAS+NOX2",     data = bh)

Coefficients:
(Intercept)          RM2          AGE       logDIS       logRAD          TAX  
  9.756e+00    6.328e-03    9.074e-05   -1.913e-01    9.571e-02   -4.203e-04  
    PTRATIO            B     logLSTAT         CRIM           ZN        INDUS  
 -3.112e-02    3.637e-01   -3.712e-01   -1.186e-02    8.016e-05    2.395e-04  
       CHAS         NOX2  
  9.140e-02   -6.380e-03  

> summary(hedonic)

Call:
lm(formula = "logMV ~ RM2+AGE+logDIS+logRAD+TAX+PTRATIO+B+logLSTAT+CRIM+ZN+INDUS+CHAS+NOX2", 
    data = bh)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.71176 -0.09169 -0.00566  0.09895  0.79780 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  9.756e+00  1.496e-01  65.221  < 2e-16 ***
RM2          6.328e-03  1.312e-03   4.823 1.89e-06 ***
AGE          9.074e-05  5.263e-04   0.172 0.863179    
logDIS      -1.913e-01  3.339e-02  -5.727 1.78e-08 ***
logRAD       9.571e-02  1.913e-02   5.002 7.91e-07 ***
TAX         -4.203e-04  1.227e-04  -3.426 0.000664 ***
PTRATIO     -3.112e-02  5.013e-03  -6.208 1.14e-09 ***
B            3.637e-01  1.031e-01   3.527 0.000460 ***
logLSTAT    -3.712e-01  2.501e-02 -14.841  < 2e-16 ***
CRIM        -1.186e-02  1.245e-03  -9.532  < 2e-16 ***
ZN           8.016e-05  5.056e-04   0.159 0.874105    
INDUS        2.395e-04  2.364e-03   0.101 0.919318    
CHAS         9.140e-02  3.320e-02   2.753 0.006129 ** 
NOX2        -6.380e-03  1.131e-03  -5.639 2.88e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.1825 on 492 degrees of freedom
Multiple R-squared: 0.8059,	Adjusted R-squared: 0.8008 
F-statistic: 157.1 on 13 and 492 DF,  p-value: < 2.2e-16 

> 
> 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()
null device 
          1 
> 
> 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 =")
[1] "demand ="
> print(demand)

Call:
lm(formula = "logWhat ~ logNOX + logMV", data = bh)

Coefficients:
(Intercept)       logNOX        logMV  
    -1.2563       0.7214       0.7351  

> summary(demand)

Call:
lm(formula = "logWhat ~ logNOX + logMV", data = bh)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.63319 -0.08263  0.00447  0.09130  0.56668 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.25628    0.23732  -5.294 1.79e-07 ***
logNOX       0.72136    0.03985  18.100  < 2e-16 ***
logMV        0.73512    0.01964  37.422  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.1546 on 503 degrees of freedom
Multiple R-squared: 0.736,	Adjusted R-squared: 0.7349 
F-statistic: 701.1 on 2 and 503 DF,  p-value: < 2.2e-16 

> 
> 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()
null device 
          1 
> 
> 
> proc.time()
   user  system elapsed 
  0.345   0.024   0.563