# QUANTILE REGRESSION AND MCMILLEN'S LWR QUANTILE MODEL ############################
------------------------------------------------------------------------------------
# Set working directory to where data file is located
# setwd("~/.......")
  
library(car)
library(spdep)
library(lmtest)
library(sandwich)
library(rgdal)
library(RColorBrewer)
library(maptools)
library(lattice)
library(quantreg)
library(sphet)
library(GWmodel)
library(spatialprobit)
library(splm)
library(quantreg)
library(sm)
library(locfit)
library(RANN)
# install.packages("~/McSpatial_2.0.tar.gz") from working directory
library(McSpatial)

# 0. INTRODUCTION ##############
--------------------------------

# Read the data
cityhouse <- read.delim("cityhouse.txt")
View(cityhouse)
attach(cityhouse)


# 1. OLS AND STANDARD QUANTILE MODELS ##############################################
------------------------------------------------------------------------------------

# 1.1. Data treatment (if necessary) ###############################################
     
# 1.1.1. Create a dataframe (complete database) ####################################

qrdb <- data.frame("LPRICE"=cityhouse$LPRICE, "FL_1"=cityhouse$FL_1,
  "ATTIC"=cityhouse$ATTIC,   "HOUSE"=cityhouse$HOUSE,
  "DUPLEX"=cityhouse$DUPLEX, "BEDSIT"=cityhouse$BEDSIT,
  "LM2"=cityhouse$LM2, "REFORM"=cityhouse$REFORM, "NEW"=cityhouse$NEW,
  "CONT"=cityhouse$CONT, "NOISE"=cityhouse$NOISE,
  "SERVICE"=cityhouse$SERVICE, "DISCEN"=cityhouse$DISCEN,
  "MINMETRO"=cityhouse$MINMETRO, "MINHUB"=cityhouse$MINHUB,
  "MINM30"=cityhouse$MINM30, "MINM40"=cityhouse$MINM40,
  "MINAIRP"=cityhouse$MINAIRP, "DISPARK"=cityhouse$DISPARK,
  "BEXPT"=cityhouse$BEXPT, "FSOCIO"=cityhouse$FSOCIO,
  "EDUCK"=cityhouse$EDUCK, "PERIPH"=cityhouse$PERIPH,
  "XLONG"=cityhouse$XLONG, "YLAT"=cityhouse$XLAT,
  "POPDENS"=cityhouse$POPDENS,"PPOP2Y"=cityhouse$PPOP2Y,
  "PPOP5Y"=cityhouse$PPOP5Y, "PPOP10Y"=cityhouse$PPOP10Y,
  "PD11001"=cityhouse$PD11001, "NIGHT"=cityhouse$NIGHT,
  "RESNOISE"=cityhouse$RESNOISE, "RESCONT"=cityhouse$RESCONT,
  "SERVK"=cityhouse$SERVK, "PPOP10YK"=cityhouse$PPOP10YK,
  "RESMEASUR"=cityhouse$RESMEASUR, "DIFMEASUR"=cityhouse$DIFMEASUR,
  "MINEXIT2"=cityhouse$MINEXIT2, "MEDEXIT2"=cityhouse$MEDEXIT2,
  "MINMETR2"=cityhouse$MINMETR2, "MINHUB0"=cityhouse$MINHUBO,
  "MINHUB20"=cityhouse$MINHUB2O,"MINMETRP"=cityhouse$MINMETRP,
  "MINHUB2"=cityhouse$MINHUB2, "MINHUB3"=cityhouse$MINHUB3,
  "MINAIRP2"=cityhouse$MINAIRP2, "MINM40P"=cityhouse$MINM40P,
  "MINAIRP2"=cityhouse$MINAIRP2, "MINAIRP3"=cityhouse$MINAIRP3,
  "MINAIRBF"=cityhouse$MINAIRBF)

# 1.1.2. Create a model formula "form" #############################################

form <- LPRICE~FL_1+ATTIC+HOUSE+DUPLEX+BEDSIT+LM2+REFORM+NEW+
  DISCEN+MINAIRP+MINMETRO+MINMETR2+MINHUB+MINHUB2+MINM30+
  MINEXIT+BEXPT+SERVK+EDUCK+PPOP10Y+CONT+NOISE+PERIPH

# 1.2. OLS Basic model #############################################################

ols <- lm(form)
summary(ols)

     
# 1.3. Standard QR #################################################################

# 1.3.1. For individual quantiles ##################################################

qr10 <- rq(form, tau=.10, data=qrdb)
summary(qr10)
errqr10 <- residuals(qr10)
errqr10.tab <- errqr10 ; write.table(errqr10.tab, "errqr10.tab.txt", sep="\t") 

qr25 <- rq(form, tau=.25, data=qrdb)
summary(qr25)
errqr25 <- residuals(qr25)
errqr25.tab <- errqr25 ; write.table(errqr25.tab, "errqr25.tab.txt", sep="\t") 

qr50 <- rq(form, tau=.50, data=qrdb)
summary(qr50)
errqr50 <- residuals(qr50)
errqr50.tab <- errqr50 ; write.table(errqr50.tab, "errqr50.tab.txt", sep="\t") 

qr75 <- rq(form, tau=.75, data=qrdb)
summary(qr75)
errqr75 <- residuals(qr75)
errqr75.tab <- errqr75 ; write.table(errqr75.tab, "errqr75.tab.txt", sep="\t") 

qr90 <- rq(form, tau=.90, data=qrdb)
summary(qr90)
errqr90 <- residuals(qr90)
errqr90.tab <- errqr90 ; write.table(errqr90.tab, "errqr90.tab.txt", sep="\t") 

# 1.3.2 For multiple quantiles ######################################################

# Main results
param_qr <- rq(form,tau=seq(.05,.95,.05), data=qrdb)
summary(param_qr, se="nid")

# Coefficient summaries
param_qr_2 <- qregbmat(form,taumat=seq(.10,.90,.10),
                graphb=FALSE, graph.factor=TRUE,data=qrdb)
summary(param_qr_2, se="nid")   

# Coefficient graphs
param_qr <- qregbmat(form,taumat=seq(.01,.99,.01),graphb=TRUE,
  graph.factor=TRUE,data=qrdb)


# 2. LOCALLY WEIGHTED REGRESSION QUANTILE MODEL ####################################
# Conditionally parametric LWR quantile estimation
------------------------------------------------------------------------------------

# 2.1. Model estimation ############################################################
# Euclidean distance - tricube kernel - 30% window size

e <-system.time(     
  lwr_cpar_et30 <- qregcpar(LPRICE~FL_1+ATTIC+HOUSE+DUPLEX+BEDSIT+
          LM2+REFORM+NEW+CONT+NOISE+EDUCK+SERVK+PPOP10YK,
          nonpar=~XLONG+YLAT,
          taumat=c(.10,.20,.30,.40,.50,.60,.70,.80,.90), kern="tcub",
          window=.30, distance="EUCLID", data=qrdb))

coef <- as.data.frame(lwr_cpar_et30$xcoef)
coef.intercept <- coef[,1:9] ; write.table(coef.intercept, "coef.intercept.txt", sep="\t") 
coef.FL_1 <- coef[,10:18] ; write.table(coef.FL_1, "coef.FL_1.txt", sep="\t") 
coef.ATTIC <- coef[,19:27] ; write.table(coef.ATTIC, "coef.ATTIC.txt", sep="\t") 
coef.HOUSE <- coef[,28:36] ; write.table(coef.HOUSE, "coef.HOUSE.txt", sep="\t") 
coef.DUPLEX <- coef[,37:45] ; write.table(coef.DUPLEX, "coef.DUPLEX.txt", sep="\t") 
coef.BEDSIT <- coef[,46:54] ; write.table(coef.BEDSIT, "coef.BEDSIT.txt", sep="\t") 
coef.LM2 <- coef[,55:63] ; write.table(coef.LM2, "coef.LM2.txt", sep="\t") 
coef.REFORM <- coef[,64:72] ; write.table(coef.REFORM, "coef.REFORM.txt", sep="\t") 
coef.NEW <- coef[,73:81] ; write.table(coef.NEW, "coef.NEW.txt", sep="\t") 
coef.CONT <- coef[,82:90] ; write.table(coef.CONT, "coef.CONT.txt", sep="\t") 
coef.NOISE <- coef[,91:99] ; write.table(coef.NOISE, "coef.NOISE.txt", sep="\t") 
coef.EDUCK <- coef[,100:108] ; write.table(coef.EDUCK, "coef.EDUCK.txt", sep="\t") 
coef.SERVK <- coef[,109:117] ; write.table(coef.SERVK, "coef.SERVK.txt", sep="\t") 
coef.PPOP10YK <- coef[,118:126] ; write.table(coef.PPOP10YK, "coef.PPOP10YK.txt", sep="\t") 

# 2.2. Kernel Density Plots #######################################################

# 2.2.1. For CONT #################################################################

# Stack the columns for quantiles:
names(coef.CONT)<-c("Q10","Q20","Q30","Q40","Q50","Q60","Q70","Q80","Q90")
vec.CONT <- stack(coef.CONT)

# Set the corresponding categories for each graph:
low <- subset(vec.CONT,vec.CONT$ind %in% c("Q10","Q20","Q30","Q40"))
med <- subset(vec.CONT,vec.CONT$ind %in% c("Q50"))
high <- subset(vec.CONT,vec.CONT$ind %in% c("Q60","Q70","Q80","Q90"))
all <- subset(vec.CONT,vec.CONT$ind %in% c("Q10","Q30","Q50","Q70", "Q90"))

# Plot for all the categories (only 5 categories)
all$ind <- factor(all$ind)
sm.density.compare(all$values, all$ind,
          xlab="CONT coefficient for all prices")
colfill<-c(2:(2+length(levels(all$ind))))
legend(locator(1), levels(all$ind), fill=colfill) 

# Plot for low
low$ind <- factor(low$ind)
sm.density.compare(low$values, low$ind,
          xlab="Air pollution coefficient for low prices")
colfill<-c(2:(2+length(levels(low$ind))))
legend(locator(1), levels(low$ind), fill=colfill) 

# Plot for median
d <- density(med$values, xlab="Air pollution coefficient for the median")
plot(d)

# Plot for high
high$ind <- factor(high$ind)
sm.density.compare(high$values, high$ind, xlab="Air pollution coefficient for high prices")
colfill<-c(2:(2+length(levels(high$ind))))
legend(locator(1), levels(high$ind), fill=colfill) 

# 2.2.2. For NOISE #################################################################

# Stack the columns for quantiles:
names(coef.NOISE)<-c("Q10","Q20","Q30","Q40","Q50","Q60","Q70","Q80","Q90")
vec.NOISE <- stack(coef.NOISE)

# Set the corresponding categories for each graph:
low <- subset(vec.NOISE,vec.NOISE$ind %in% c("Q10","Q20","Q30","Q40"))
med <- subset(vec.NOISE,vec.NOISE$ind %in% c("Q50"))
high <- subset(vec.NOISE,vec.NOISE$ind %in% c("Q60","Q70","Q80","Q90"))
all <- subset(vec.NOISE,vec.NOISE$ind %in% c("Q10","Q30","Q50","Q70", "Q90"))

# Plot for all the categories (only 5 categories)
all$ind <- factor(all$ind)
sm.density.compare(all$values, all$ind,
          xlab="Noise coefficient for all prices")
colfill<-c(2:(2+length(levels(all$ind))))
legend(locator(1), levels(all$ind), fill=colfill) 

# Plot for low
low$ind <- factor(low$ind)
sm.density.compare(low$values, low$ind,
          xlab="Noise coefficient for low prices")
colfill<-c(2:(2+length(levels(low$ind))))
legend(locator(1), levels(low$ind), fill=colfill) 

# Plot for median
d <- density(med$values, xlab="Noise coefficient for the median")
plot(d)

# Plot for high
high$ind <- factor(high$ind)
sm.density.compare(high$values, high$ind, xlab="Noise coefficient for high prices")
colfill<-c(2:(2+length(levels(high$ind))))
legend(locator(1), levels(high$ind), fill=colfill) 

# 2.2.3. For LM2 ###################################################################

# Stack the columns for quantiles:
names(coef.LM2)<-c("Q10","Q20","Q30","Q40","Q50","Q60","Q70","Q80","Q90")
vec.LM2 <- stack(coef.LM2)

# Set the corresponding categories for each graph:
low <- subset(vec.LM2,vec.LM2$ind %in% c("Q10","Q20","Q30","Q40"))
med <- subset(vec.LM2,vec.LM2$ind %in% c("Q50"))
high <- subset(vec.LM2,vec.LM2$ind %in% c("Q60","Q70","Q80","Q90"))
all <- subset(vec.LM2,vec.LM2$ind %in% c("Q10","Q30","Q50","Q70", "Q90"))

# Plot for all the categories (only 5 categories)
all$ind <- factor(all$ind)
sm.density.compare(all$values, all$ind,
                   xlab="LM2 coefficient for all prices")
colfill<-c(2:(2+length(levels(all$ind))))
legend(locator(1), levels(all$ind), fill=colfill) 

# Plot for low
low$ind <- factor(low$ind)
sm.density.compare(low$values, low$ind,
                   xlab="LM2 coefficient for low prices")
colfill<-c(2:(2+length(levels(low$ind))))
legend(locator(1), levels(low$ind), fill=colfill) 

# Plot for median
d <- density(med$values, xlab="LM2 coefficient for the median")
plot(d)

# Plot for high
high$ind <- factor(high$ind)
sm.density.compare(high$values, high$ind, xlab="LM2 coefficient for high prices")
colfill<-c(2:(2+length(levels(high$ind))))
legend(locator(1), levels(high$ind), fill=colfill)