#####################################

#Set up

•rm(list=ls(all=TRUE))

•library(ggplot2)

•library(quantreg)

•theme_set(theme_bw())

•setwd("I:/Project/2017 LD GHG Rule/Footprint_analysis/Outputs")

•fleetorig <- read.table("I:/Project/2017 LD GHG 

Rule/Footprint_analysis/fleet/20110520_fleet.txt", header=TRUE, sep="\t", 

na.strings="NA", dec=".", strip.white=TRUE)

#Drop Hybrids, Electrics and Diesels

fleet<-fleetorig[fleetorig$Engine.Cycle %in% c("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.v

ol,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))}

•

annotrq<-function(a)

•{paste(signif(coef(a)[[1]],digits=3),"+",signif(coef(a)[[2]],digits=3),"x")}

•

analysis<-function(x,col,title) {

•fleet1<-subset(fleet,NHTSA.Defined.New.NHTSA.Car.Truck == x)

•CO2<-cbind(fleet1$No.Tech.CO2,fleet1$CO2.Package.1,fleet1$CO2.Package.2)

•lm.1<-lm(fleet1$Horsepower/fleet1$Curb.Weight ~ fleet1$Footprint..PU.Average, 

weights=fleet1$Weight)

•lm.2<-lm(fleet1$Horsepower/fleet1$Curb.Weight ~ 

poly(fleet1$Footprint..PU.Average,2,raw=TRUE), weights=fleet1$Weight)

•lm.3<-lm(fleet1$Horsepower~poly(fleet1$Curb.Weight,1, raw=TRUE), weights=fleet1$Weight)

•lm.4<-lm(fleet1$Horsepower~poly(fleet1$Curb.Weight,2, raw=TRUE), weights=fleet1$Weight)

•lm.5<-lm(CO2[,col] ~ fleet1$Footprint..PU.Average, weights=fleet1$Weight)

•

#Rick's method for adjusting hp.wt

•lm.8<-lm(CO2[,col] ~ 

I(fleet1$Horsepower/fleet1$Curb.Weight)+I(fleet1$Curb.Weight),weights = fleet1$Weight)

•#Remove HP.Wt Adjustment,

•#fleet1$hp.wt.avg<-sum(fleet1$Horsepower/fleet1$Curb.Weight*fleet1$Weight)

•fleet1$hp.wt.CO2<-CO2[,col]#-(fleet1$Horsepower/fleet1$Curb.Weight-fleet1$hp.wt.avg)*coef

(lm.8)[2])

•lm.9<-lm(fleet1$hp.wt.CO2 ~ fleet1$Footprint..PU.Average, weights=fleet1$Weight)

•lm.10<-lm(CO2[,col] ~ fleet1$Footprint..PU.Average + 

I(fleet1$Horsepower/fleet1$Curb.Weight)+I(fleet1$Curb.Weight/fleet1$Footprint..PU.Average)

,weights = fleet1$Weight)

•

#Adjustment for HP.WT and FP.WT

•fleet1$crb.wt.fp.avg<-sum(fleet1$Curb.Weight/fleet1$Footprint..PU.Average*fleet1$Weight)

•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)

•lm.11<-lm(fleet1$hpwt.crbwt.norm.CO2~ fleet1$Footprint..PU.Average, 

weights=fleet1$Weight)

#FP/WT v. FP.

•lm.12<-lm(I(fleet1$Curb.Weight/fleet1$Footprint..PU.Average) ~ 

fleet1$Footprint..PU.Average, weights = fleet1$Weight)

•

#CO2 v FP + HP/WT + WT/FP•
•lm.13<-lm(CO2[,col] ~ I(fleet1$Curb.Weight/fleet1$Footprint..PU.Average) + 

fleet1$Footprint..PU.Average + I(fleet1$Horsepower/fleet1$Curb.Weight), weights = 

fleet1$Weight)

•lm.14<-lm(CO2[,col] ~ I(fleet1$Curb.Weight/fleet1$Footprint..PU.Average) + 

poly(fleet1$Footprint..PU.Average,2,raw=TRUE) + I(fleet1$Horsepower/fleet1$Curb.Weight), 

weights = fleet1$Weight)

•lm.15<-lm(CO2[,col] ~ I(fleet1$Curb.Weight/fleet1$Footprint..PU.Average) + 

poly(fleet1$Footprint..PU.Average,3,raw=TRUE) + I(fleet1$Horsepower/fleet1$Curb.Weight), 

weights = fleet1$Weight)

•

•••

•••

•ans<-list(fleet1,lm.1,lm.2,lm.3,lm.4,lm.5,lm.8,lm.9,lm.10,lm.11,lm.12,lm.13,lm.14,lm.15)

•invisible(return(ans))

}

••

mregres<-function(vtype,a){•

fleeta<-subset(fleet4,NHTSA.Defined.New.NHTSA.Car.Truck == vtype)

#####WT/FP

olsa1<-lm( I(fleeta$Curb.Weight/fleeta$Footprint..PU.Average)~ 

fleeta$Footprint..PU.Average,weights = fleeta$Weight)

#####HP/WT

olsb1<-lm( I(fleeta$Horsepower/fleeta$Curb.Weight)~ fleeta$Footprint..PU.Average,weights 

= fleeta$Weight)

return(list(olsa1,olsb1))

}

dchart<-function(dtt,tit,hj){

#Unadjusted CO2

dtt$art<-fleet4$CO2.Package.1

clabels<-list(annot(cars.no.tech[[6]]),annot(cars.egr[[6]]))

tlabels<-list(annot(trucks.no.tech[[6]]),annot(trucks.egr[[6]]))

#Base CO2 - Weighted OLS (faceted B& W)

base.chart<-ggplot(data=dtt, aes(y=art, x=Footprint..PU.Average, size=Volume.2021, 

weight=Weight))+

••geom_point(aes(shape=1))+

••labs(x="Footprint", y="Technology Normalized 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,500))+

••scale_x_continuous(limits=c(35,80))

#Base CO2 - Weighted OLS

base.ols<-ggplot(data=dtt, aes(y=art, x=Footprint..PU.Average, 

color=NHTSA.Defined.New.NHTSA.Car.Truck, size=Volume.2021, weight=Weight))+

••geom_point(aes(shape=1))+

••geom_smooth(method=lm, formula= y~ poly(x,1, raw=TRUE), se=FALSE)+

••labs(x="Footprint", y="Base CO2", colour="C/T")+

••opts(title=paste(tit,"Base CO2 v. Footprint (Wt. OLS)"))+

••coord_cartesian(xlim = c(35,80), ylim=c(0,500))+

••annotate(geom="text",x=50,y=50, colour="black", label= paste("Car:",clabels[[hj]]) , 

size=3 ) +

••annotate(geom="text", x = 50, y= 100, colour="black", 

label=paste("Truck:",tlabels[[hj]]) , size=3 )
#WT/FP Normalized CO2 - OLS

clabels<-list(annot(cars.no.tech[[10]]),annot(cars.egr[[10]]))

tlabels<-list(annot(trucks.no.tech[[10]]),annot(trucks.egr[[10]]))

hp.fp.ols<-•ggplot(data=dtt, aes(y=hpwt.crbwt.norm.CO2, x=Footprint..PU.Average, 

color=NHTSA.Defined.New.NHTSA.Car.Truck, size=Volume.2021, weight=Weight))+

•geom_point(aes(shape=1))+

•geom_smooth(method=lm, formula= y~ poly(x,1, raw=TRUE),se=FALSE)+

•labs(x="Footprint", y="WT/FP Normalized CO2", colour="C/T")+

•opts(title=paste(tit,"WT/FP Normalized CO2 v. Footprint (Wt. OLS)"))•+

•coord_cartesian(xlim = c(35,80), ylim=c(0,500))+

•annotate(geom="text",x=50,y=50, colour="black", label= paste("Car:",clabels[[hj]]) , 

size=3 ) +

•annotate(geom="text", x = 50, y= 100, colour="black", 

label=paste("Trucks:",tlabels[[hj]]) , size=3 )

•

•return(list(base.chart,base.ols,hp.fp.ols))

•

}

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")

#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")

•

# 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")

return(list(p1,p2,p4))

}

cars.no.tech<-analysis("C",1,"No Tech")

trucks.no.tech<-analysis("T",1,"No Tech")

cars.egr<-analysis("C",2,"EGR")

trucks.egr<-analysis("T",2,"EGR")

fleet3<-data.frame(rbind(cars.no.tech[[1]], trucks.no.tech[[1]]))

fleet4<-data.frame(rbind(cars.egr[[1]], trucks.egr[[1]]))
pltcars.no.tech<-mregres("C",1)

plttks.no.tech<-mregres("T",1)

pltcars.egr<-mregres("C",2)

plttks.egr<-mregres("T",2)

pdf(paste("WT.FP_curves_for_analysis-",Sys.Date(),".pdf",sep=""))

dchart(fleet4, "Max ICE Tech -",2)

physics()

fleet$sales.2008<-fleet$Total.Production.Volume

ggplot(data=fleet, aes(y=CO2, x=Footprint..PU.Average, size=sales.2008))+

•geom_point(shape=1)+

•labs(x="Footprint", y="2008 CO2", colour="Car/Truck", size="2008 Sales")+

•opts(title=paste("2008 CO2 v. Footprint"))+

•facet_wrap(~NHTSA.Defined.New.NHTSA.Car.Truck)+

•coord_cartesian(xlim = c(35,80), ylim=c(0,700))

summary(cars.egr[[7]])

summary(trucks.egr[[7]])

sd(cars.egr[[1]]$CO2)

sd(cars.egr[[1]]$CO2.Package.1)

sd(trucks.egr[[1]]$CO2)

sd(trucks.egr[[1]]$CO2.Package.1)

#Get Densities

man.sales<-aggregate((fleet$Volume.2021), by=list(fleet$Our.Class), sum)

colnames(man.sales)<- c("Our.Class", "Man.2021.Sales")

man.footprint<-stack(by(fleet,list(fleet$Our.Class),function(subset)weighted.mean(subset$F

ootprint..PU.Average, subset$Volume.2021), simplify=TRUE))

colnames(man.footprint)<- c("Avg.Footprint","Our.Class")

•tyu2<-merge(man.sales,fleet)•

•tyu2$man.proportion<-tyu2$Volume.2021/tyu2$Man.2021.Sales

•tyu2$wt.fp.proportion<-tyu2$Curb.Weight/tyu2$Footprint..PU.Average*tyu2$man.proportion

•manwtfp<-aggregate((tyu2$wt.fp.proportion), by=list(tyu2$Our.Class), sum)

•colnames(manwtfp) <- c("Our.Class", "Avg.Density")

•manwtfp<- transform(manwtfp, Our.Class=reorder(Our.Class,Avg.Density))

•testmerge<-merge(manwtfp,man.sales)

•manwtfp<-merge(testmerge, man.footprint)

•#wt/fp averages

•

ct<-signif(sum(cars.egr[[1]]$Curb.Weight/cars.egr[[1]]$Footprint..PU.Average*cars.egr[[1]]

$Weight),digits=3)

dt<-signif(sum(trucks.egr[[1]]$Curb.Weight/trucks.egr[[1]]$Footprint..PU.Average*trucks.eg

r[[1]]$Weight), digits=3)

•ggplot(data=manwtfp,aes(x=Avg.Footprint, y=Avg.Density,label=Our.Class))+

••geom_point(size=1)+

••geom_text(hjust=-.05, angle=30, size=3)+

••labs(x="Footprint",y="WT/FP")+

••opts(title=paste("WT/FP by Vehicle Class"))+

••opts(axis.text.x=theme_text(hjust=1))+

••geom_hline(yintercept=ct, linetype=2, colour="red")+

••annotate(geom="text", x = 65, y= 71, colour="red",label=paste("Car fleet avg.:",ct), 

size=3)+

••geom_hline(yintercept=dt, linetype=2, colour="blue")+

••annotate(geom="text", x = 65, y= 83, colour="blue",label=paste("Truck fleet avg.:",dt), 

size=3)+

••scale_x_continuous(limits=c(40,70))+

••scale_y_continuous(limits=c(65,100))

•

#Footprint Histogram•
•ggplot(data=fleet,aes(x=Footprint..PU.Average, weight=Volume.2021))+

••geom_bar(color="darkblue", 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, colour="darkgreen")+

••geom_vline(xintercept=56, linetype=2, colour="darkgreen")+

••geom_vline(xintercept=74, linetype=4, colour="red")+

••geom_vline(xintercept=41, linetype=4, colour="red")

#wt/fp co2

fleetwadj<-rbind(cars.egr[[1]],trucks.egr[[1]])

ggplot(data=fleetwadj, aes(y=hpwt.crbwt.norm.CO2/8887, x=Footprint..PU.Average, 

size=Volume.2021, weight=Weight))+

•geom_point(aes(shape=1))+

•facet_wrap(~NHTSA.Defined.New.NHTSA.Car.Truck)+

•geom_smooth(method=lm, formula= y~ poly(x,1, raw=TRUE), se=FALSE)+

•labs(x="Footprint", y="GPM", colour="C/T")+

•opts(title=paste("WT/FP Normalized Fuel Consumption v. Footprint"))+

•scale_y_continuous(limits=c(0,0.06))•

dev.off()