#Set up
	rm(list=ls(all=TRUE))
	library(ggplot2)
	library(quantreg)
	library(gplots)
	library(xlsReadWrite)
	theme_set(theme_bw())
	setwd("I:/Project/2017 LD GHG Rule/Footprint_analysis/Outputs")
	fleetorig <- xlsReadWrite::read.xls("I:/Project/2017 LD GHG Rule/Footprint_analysis/fleet/2008_20120509_fleet.xls")

#Drop Electrics and Diesels -Leave Hybrids per NHTSA NPRM
fleet<-fleetorig[fleetorig$Engine.Cycle %in% c("H","G","R"),]

##Calculate Weights - Don't remove vehicles after this, as these are based on total volumes
	a<-subset(fleet,NHTSA.Defined.New.NHTSA.Car.Truck == "C")
	car.vol<-sum(a$Volume.2021)
	b<-subset(fleet,NHTSA.Defined.New.NHTSA.Car.Truck == "T")
	truck.vol<-sum(b$Volume.2021)
	fleet$Weight<-ifelse(fleet$NHTSA.Defined.New.NHTSA.Car.Truck=="C",fleet$Volume.2021/car.vol,fleet$Volume.2021/truck.vol)
	
annot<-function(a)
	{paste(signif(coef(a)[[1]],digits=3),"+",signif(coef(a)[[2]],digits=3),"x : R2=",signif(summary(a)[[8]],digits=3),sep = "")}
	
annotrq<-function(a)
	{paste(signif(coef(a)[[1]],digits=3),"+",signif(coef(a)[[2]],digits=3),"x",sep = "")}
	

analysis<-function(x) {
	#Create Fleet
	
	fleet1<-subset(fleet,NHTSA.Defined.New.NHTSA.Car.Truck == x)
		
	#Diagnostics
	#HP/WT ~ Footprint
	lm.1<-lm(fleet1$Horsepower/fleet1$Curb.Weight ~ fleet1$Footprint..PU.Average, weights=fleet1$Weight)
	#HP/WT ~ Footprint Squared
	lm.2<-lm(fleet1$Horsepower/fleet1$Curb.Weight ~ poly(fleet1$Footprint..PU.Average,2,raw=TRUE), weights=fleet1$Weight)
	#WT/FP ~ FP.
	lm.12<-lm(I(fleet1$Curb.Weight/fleet1$Footprint..PU.Average) ~ fleet1$Footprint..PU.Average, weights = fleet1$Weight)
	#HP ~ WT 
	lm.3<-lm(fleet1$Horsepower~poly(fleet1$Curb.Weight,1, raw=TRUE), weights=fleet1$Weight)
	#HP ~WT Squared
	lm.4<-lm(fleet1$Horsepower~poly(fleet1$Curb.Weight,2, raw=TRUE), weights=fleet1$Weight)
	#CO2 ~ Footprint
	lm.5<-lm(fleet1$CO2.Package.1 ~ fleet1$Footprint..PU.Average, weights=fleet1$Weight)
	#CO2 ~ HP/WT + WT
	lm.8<-lm(fleet1$CO2.Package.1 ~ I(fleet1$Horsepower/fleet1$Curb.Weight)+I(fleet1$Curb.Weight),weights = fleet1$Weight)

	#Create CO2 w/ Adjustments
	#HP/WT Adjusted CO2
	fleet1$hp.wt.avg<-sum(fleet1$Horsepower/fleet1$Curb.Weight*fleet1$Weight)
	fleet1$hp.wt.CO2<-(fleet1$CO2.Package.1-(fleet1$Horsepower/fleet1$Curb.Weight-fleet1$hp.wt.avg)*coef(lm.8)[2])
	#WT/FP adjusted CO2
	fleet1$crb.wt.fp.avg<-sum(fleet1$Curb.Weight/fleet1$Footprint..PU.Average*fleet1$Weight)
	fleet1$crbwt.norm.CO2<-(fleet1$CO2.Package.1-coef(lm.8)[3]*(fleet1$Curb.Weight-fleet1$crb.wt.fp.avg*fleet1$Footprint..PU.Average))
	#HP/WT and WT/FP Adjusted CO2
	fleet1$hpwt.crbwt.norm.CO2<-(fleet1$hp.wt.CO2-coef(lm.8)[3]*(fleet1$Curb.Weight-fleet1$crb.wt.fp.avg*fleet1$Footprint..PU.Average))

	#Return Answers
	ans<-list(fleet1,lm.1,lm.2,lm.12,lm.3,lm.4,lm.5,lm.8)
	invisible(return(ans))
}

#Function for the regressions needed on each of the CO2 datasets.
regressions<-function(dataset){
	datasetb<-data.frame(dataset[1])
	#List the datasets here
	CO2set<-data.frame(datasetb$CO2, datasetb$No.Tech.CO2, datasetb$CO2.Package.1, datasetb$hp.wt.CO2, datasetb$crbwt.norm.CO2, datasetb$hpwt.crbwt.norm.CO2)
		runr<-function(dsb){
			weighted.ols	<-lm(dsb~datasetb$Footprint..PU.Average,weights = datasetb$Weight)
			unweighted.ols	<-lm(dsb~datasetb$Footprint..PU.Average,weights = NULL)
			weighted.mad	<-rq(dsb~datasetb$Footprint..PU.Average, tau=.5, weights = datasetb$Weight)
			unweighted.mad	<-rq(dsb~datasetb$Footprint..PU.Average, tau=.5,weights = NULL)
		invisible(return(regr<-list(weighted.ols,unweighted.ols,weighted.mad,unweighted.mad)))
		}
	invisible(return(j<-lapply(CO2set,runr)))
	}

#Create the Datasets, and excute the regressions, then merge the car and truck fleets
cars.egr<-analysis("C")
car.regressions<-regressions(cars.egr)
trucks.egr<-analysis("T")
truck.regressions<-regressions(trucks.egr)
fleet4<-data.frame(rbind(cars.egr[[1]], trucks.egr[[1]]))

	
car_base_equation<-annot(cars.egr[[8]])
truck_base_equation<-annot(trucks.egr[[8]])
		

summarize_regressions<-function(name,number){		
	name<-data.frame(rbind(
		cbind(name,"weighted OLS",annot(car.regressions[[number]][[1]]),annot(truck.regressions[[number]][[1]])),
		cbind(name,"unweighted OLS",annot(car.regressions[[number]][[2]]),annot(truck.regressions[[number]][[2]])),
		cbind(name,"weighted MAD",annotrq(car.regressions[[number]][[3]]),annotrq(truck.regressions[[number]][[3]])),
		cbind(name,"unweighted MAD",annotrq(car.regressions[[number]][[4]]),annotrq(truck.regressions[[number]][[4]]))
		))
	
	return(name)
	}
	
all_regressions<-rbind(	
summarize_regressions("2008",1),
summarize_regressions("No.Tech",2),
summarize_regressions("Tech.Normalized",3),
summarize_regressions("HP/WT-Tech-Normalized",4),
summarize_regressions("WT/FP-Tech-Normalized",5),
summarize_regressions("HP/WT&WT/FP-Tech-Normalized",6)
)
colnames(all_regressions)<-c("Dataset","Regression","Car Equation", "Truck Equation")
all_regressions


dchart<-function(tit,jj,mult){ 
	#Base CO2
	fleet4$RCO2<-fleet4[,jj]
	base.chart<-ggplot(data=fleet4, aes(y=RCO2, x=Footprint..PU.Average,  size=Volume.2021, weight=Weight))+
		geom_point(aes(shape=1))+
		labs(x="Footprint", y="CO2", size="2021 Projected Sales")+
		opts(title=paste(tit,"CO2 v. Footprint"))+
		facet_wrap(~NHTSA.Defined.New.NHTSA.Car.Truck)+
		#scale_y_continuous(limits=c(0,800))+
		scale_x_continuous(limits=c(35,80))+
		#weighted lines
		geom_smooth(data = subset(fleet4, NHTSA.Defined.New.NHTSA.Car.Truck =="C"),method=lm, formula= y~ poly(x,1, raw=TRUE), se=FALSE) + 
		geom_smooth(data = subset(fleet4, NHTSA.Defined.New.NHTSA.Car.Truck =="T"),method=lm, formula= y~ poly(x,1, raw=TRUE), se=FALSE) +
		geom_quantile(data = subset(fleet4, NHTSA.Defined.New.NHTSA.Car.Truck =="C"),quantiles=.5, colour="red")+
		geom_quantile(data = subset(fleet4, NHTSA.Defined.New.NHTSA.Car.Truck =="T"),quantiles=.5, colour="red")+
		#unweighted lines
		geom_smooth(data = subset(fleet4, NHTSA.Defined.New.NHTSA.Car.Truck =="C"),method=lm, formula= y~ poly(x,1, raw=TRUE), se=FALSE, weight=1,linetype=5) + 
		geom_smooth(data = subset(fleet4, NHTSA.Defined.New.NHTSA.Car.Truck =="T"),method=lm, formula= y~ poly(x,1, raw=TRUE), se=FALSE,weight=1,linetype=5) +
		geom_quantile(data = subset(fleet4, NHTSA.Defined.New.NHTSA.Car.Truck =="C"),quantiles=.5, colour="red",weight=1,linetype=5)+
		geom_quantile(data = subset(fleet4, NHTSA.Defined.New.NHTSA.Car.Truck =="T"),quantiles=.5, colour="red",weight=1,linetype=5)+
		#no legend
		opts(legend.position = "none")
		#add equations
		#annotate(geom="text",x=40,y=700, colour="black", data = subset(fleet4, NHTSA.Defined.New.NHTSA.Car.Truck =="C"),
		#			label= paste(all_regressions[1+(mult*4),2],all_regressions[1+(mult*4),3]) , size=2.75,hjust=0, vjust=0 )
	return(base.chart)
}



	
	physics<-function(){
#Plot of WT/ v. FP - Weighted OLS
p2<-ggplot(data=fleet, aes(y=I(fleet$Curb.Weight), x=fleet$Footprint..PU.Average, size=Volume.2021,  weight=fleet$Weight, shape=1))+
	geom_point(aes(shape=1))+
	labs(x="Footprint", y="WT", colour="C/T", size="2021 Sales")+
	geom_smooth(method=lm, formula= y~ poly(x,1, raw=TRUE), se=FALSE)+
	geom_smooth(span=1,colour="red")+
	facet_wrap(~NHTSA.Defined.New.NHTSA.Car.Truck)+
	opts(title="WT v. FP - Weighted OLS and Loess Fit")+
	opts(legend.position = "none")

#Plot of WT/FP v. FP - Weighted OLS
p1<-ggplot(data=fleet, aes(y=I(fleet$Curb.Weight/fleet$Footprint..PU.Average), x=fleet$Footprint..PU.Average, size=Volume.2021,  weight=fleet$Weight))+
	geom_point(aes(shape=1))+
	labs(x="Footprint", y="WT/FP", colour="C/T", size="2021 Sales")+
	geom_smooth(method=lm, formula= y~ poly(x,1, raw=TRUE), se=FALSE)+
	geom_smooth(span=1,colour="red")+
	facet_wrap(~NHTSA.Defined.New.NHTSA.Car.Truck)+
	opts(title="WT/FP v. FP - Weighted OLS and Loess Fit")+
	opts(legend.position = "none")
	
# Plot of HP/WT v. FP - Weighted OLS
p4<-ggplot(data=fleet, aes(y=I(fleet$Horsepower/fleet$Curb.Weight), x=fleet$Footprint..PU.Average, size=Volume.2021, weight=fleet$Weight))+
	geom_point(aes(shape=1))+
	geom_smooth(method=lm, formula= y~ poly(x,1, raw=TRUE), se=FALSE)+
	facet_wrap(~NHTSA.Defined.New.NHTSA.Car.Truck)+
	labs(x="Footprint", y="HP/WT", colour="C/T", size="2021 Sales")+
	opts(title="HP/WT v. FP - Weighted OLS")+
	opts(legend.position = "none")

return(list(p1,p2,p4))
}

#Footprint Histogram	
fphist<-ggplot(data=fleet,aes(x=Footprint..PU.Average, weight=Volume.2021))+
		geom_bar(color="black", fill="white", binwidth=1)+
		facet_wrap(~NHTSA.Defined.New.NHTSA.Car.Truck,ncol=1)+
		labs(x="Footprint",y="MY 2021 Projected Sales")+
		opts(title=paste("Footprint Distribution by Car and Truck"))+
		geom_vline(xintercept=41, linetype=2)+
		geom_vline(xintercept=56, linetype=2)+
		geom_vline(xintercept=74, linetype=6)+
		geom_vline(xintercept=41)
		


tester3<-function(df){lm(Curb.Weight ~ Footprint..PU.Average, weights=Volume.2021, data=df)}
phys1a<-dlply(fleet4,.(NHTSA.Defined.New.NHTSA.Car.Truck),tester3)[[1]]
phys1b<-dlply(fleet4,.(NHTSA.Defined.New.NHTSA.Car.Truck),tester3)[[2]]
tester2<-function(df){lm(Curb.Weight/Footprint..PU.Average ~ Footprint..PU.Average, weights=Volume.2021, data=df)}
phys2a<-dlply(fleet4,.(NHTSA.Defined.New.NHTSA.Car.Truck),tester2)[[1]]
phys2b<-dlply(fleet4,.(NHTSA.Defined.New.NHTSA.Car.Truck),tester2)[[2]]
tester<-function(df){lm(Horsepower/Curb.Weight ~ Footprint..PU.Average, weights=Volume.2021, data=df)}
phys3a<-dlply(fleet4,.(NHTSA.Defined.New.NHTSA.Car.Truck),tester)[[1]]
phys3b<-dlply(fleet4,.(NHTSA.Defined.New.NHTSA.Car.Truck),tester)[[2]]

lapply(list(phys1a,phys1b,phys2a,phys2b,phys3a,phys3b),summary)

pdf(paste("2008_curve_analysis-",Sys.Date( ),".pdf"))
textplot(all_regressions,cex=0.5,halign=c("left"))
fphist
dchart("2008","CO2",0)	
dchart("No.Tech","No.Tech.CO2",1)
dchart("Tech-Normalized","CO2.Package.1",2)
dchart("HP/WT-Tech-Normalized","hp.wt.CO2",3)
dchart("WT/FP-Tech-Normalized","crbwt.norm.CO2",4)
dchart("HP/WT&WT/FP-Tech-Normalized","hpwt.crbwt.norm.CO2",5)
physics()
dev.off()