Read and Plot the Time Series

Set the working directory appropriately first. Then read / plot the time series.

setwd("C:/Users/lfult/Desktop/Electric Car Simulation/")
mydata=read.csv("R Code/gaskwh.csv")

###############Define the Time Series##################
gas=ts(mydata$centsgal, start=c(1995, 1), frequency=12)
kw=ts(mydata$centskwh, start=c(1995,1), frequency=12)
mv=ts(mydata[,-c(1,2,3)], start=c(1995,1), frequency=12)
#######################################################

###############Describe the Data######################
library(psych)
## Warning: package 'psych' was built under R version 3.5.1
describe(gas)
##    vars   n   mean    sd median trimmed    mad min   max range skew
## X1    1 280 220.88 88.17 222.25  216.01 114.75  94 405.1 311.1 0.32
##    kurtosis   se
## X1    -1.16 5.27
describe(kw)
##    vars   n  mean  sd median trimmed  mad  min  max range skew kurtosis
## X1    1 280 10.28 1.8  10.17   10.24 2.58 7.58 13.3  5.72  0.1    -1.58
##      se
## X1 0.11
myplot=plot(gas, main="Average Retail Prices for Regular Gasoline", type="l", ylab="Cents per Gallon")

myplot
## NULL
myplot2=plot(kw, main="Average Retail Prices per kWh", type="l", ylab="Cents per kWh")

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

###############Some Basic Linear Models################
mylm1=lm(scale(mv[,1])~seq(1:280))
mylm2=lm(scale(mv[,2])~seq(1:280))
#######################################################

##################Basic Decomposition##################
mydec=decompose(kw)
plot(mydec)

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

Correlations

Conduct some basic correlational / autocorrelational analysis

#################Basic Correlation#####################
cor(gas,kw)
## [1] 0.8091864
lm(formula = scale(gas) ~ scale(kw), intercept = 0)
## Warning: In lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) :
##  extra argument 'intercept' will be disregarded
## 
## Call:
## lm(formula = scale(gas) ~ scale(kw), intercept = 0)
## 
## Coefficients:
## (Intercept)    scale(kw)  
##  -4.681e-16    8.092e-01
acf(gas)

pacf(gas)

acf(kw)

pacf(kw)

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

Time Series Functions

Generate some basic time series (ETs & ARIMA). These auto-generated models are useful enough.

#############Auto-Optimized ETS / ARIMA################
library(fpp)
## Warning: package 'fpp' was built under R version 3.5.1
## Loading required package: forecast
## Warning: package 'forecast' was built under R version 3.5.1
## 
## Attaching package: 'forecast'
## The following object is masked _by_ '.GlobalEnv':
## 
##     gas
## Loading required package: fma
## Loading required package: expsmooth
## Loading required package: lmtest
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 3.5.1
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: tseries
## Warning: package 'tseries' was built under R version 3.5.1
#for external accuracy
t1=ts(gas[1:224], start=c(1995,1), frequency=12)
tst1=gas[225:280]
t2=ts(kw[1:224], start=c(1995,1),frequency=12)
tst2=kw[225:280]
mya1=auto.arima(t1)
mye1=ets(t1)
mya2=auto.arima(t2)
mye2=ets(t2)

#total models
myarima=auto.arima(gas)
plot(forecast(myarima))

myets=ets(gas)
plot(forecast(myets))

myarima2=auto.arima(kw)
plot(forecast(myarima2))

myets2=ets(kw)
plot(forecast(myets2))

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

Internal Accuracy Metrics

Generate internal accuracy metrics.

###################Internal Accuaracy###################
myacc1=accuracy(myets)
myacc2=accuracy(myarima)
accgas=rbind(myacc1, myacc2)

myacc3=accuracy(myets2)
myacc4=accuracy(myarima2)
totalacc=rbind(accgas,myacc3,myacc4)
row.names(totalacc)=c("ETS Gas", "ARIMA Gas", "ETS kWh", "ARIMA kWh")
totalacc
##                    ME       RMSE       MAE        MPE      MAPE      MASE
## ETS Gas   0.449730367 13.0960931 8.5381625 0.08202279 3.8348781 0.2346900
## ARIMA Gas 0.666057828 12.6892972 8.8020604 0.21998747 3.8805029 0.2419438
## ETS kWh   0.004297674  0.1334283 0.1018960 0.04328157 0.9895515 0.3551585
## ARIMA kWh 0.003761583  0.1336724 0.1026441 0.04342523 0.9947409 0.3577658
##                   ACF1
## ETS Gas    0.376128593
## ARIMA Gas  0.001732673
## ETS kWh   -0.004674220
## ARIMA kWh  0.009977391
#######################################################

External Accuracy Metrics

#####################External Accuracy####################
f1=forecast(mye1, 56)
f2=forecast(mya1, 56)
f3=forecast(mye2, 56)
f4=forecast(mya2, 56)

a1=accuracy(f1,tst1)
a2=accuracy(f2,tst1)
a3=accuracy(f3,tst2)
a4=accuracy(f4,tst2)
acc=rbind(a1,a2,a3,a4)

row.names(acc)=c("ETS Gas Train", "ETS Gas Test","ARIMA Gas Train", "ARIMA Gas Test", "ETS kWh Train", "ETS kWh Test","ARIMA kWh Train", "ARIMA kWh Test")
acc
##                            ME        RMSE          MAE          MPE
## ETS Gas Train    7.535433e-01  13.1927298   8.31395274   0.20328062
## ETS Gas Test    -8.944382e+01 103.5284078  89.52606588 -38.09531368
## ARIMA Gas Train -2.885121e-01  12.3871079   8.49989969  -0.68235353
## ARIMA Gas Test  -1.279332e+02 143.6820508 127.93321549 -53.92168776
## ETS kWh Train    6.386101e-03   0.1203615   0.09446788   0.06878985
## ETS kWh Test     1.052992e-01   0.2891325   0.25130092   0.85808158
## ARIMA kWh Train  8.642786e-03   0.1251915   0.09562729   0.09089753
## ARIMA kWh Test  -7.997010e-02   0.2550320   0.20436609  -0.62560787
##                       MAPE       MASE         ACF1
## ETS Gas Train    3.8805277  0.8886175  0.374794641
## ETS Gas Test    38.1201899  9.5687848           NA
## ARIMA Gas Train  3.9606118  0.9084920 -0.011340723
## ARIMA Gas Test  53.9216878 13.6738435           NA
## ETS kWh Train    0.9706435  0.5061590  0.045742707
## ETS kWh Test     1.9971164  1.3464705           NA
## ARIMA kWh Train  0.9812281  0.5123711 -0.005756583
## ARIMA kWh Test   1.6108186  1.0949937           NA
#######################################################

Forecasts

Generate forecasts based on best-performing models (ETS for both kWh and Gas by MASE & MAPE)

##################Forecasts####################
(par(mfrow=c(2,1)))
## $mfrow
## [1] 1 1
plot(forecast(myets2, 96), main="Cents / kWh 8-Year Forecast")
plot(forecast(myets,96), main="Cents / Gallon 8-Year Forecast")

f_kw=forecast(myets2,96)$mean
f_gas=forecast(myets,96)$mean

par(mfrow=c(2,1))
plot(f_kw, main="Cents / kWh")
plot(f_gas, main="Cents / Gallon")

write.csv(f_kw,"R Code/f_kWh.csv", row.names=FALSE)
write.csv(f_gas, "R Code/f_gas.csv", row.names=FALSE)
#######################################################

Simulation

This simulation is relatively straightforward. First, arrays are initialized.

##############Initialize Arrays########################
set.seed(1)
RNGkind(kind="Mersenne-Twister",normal.kind="Box-Muller")
it=25 #number of iterations
results=rep(0,3*8*365*it) #for electric
results2=rep(0,3*8*365*it) #for hybrid
results=array(results, dim=c(3,8*365*it)) #for electric
results2=array(results2, dim=c(3,8*365*it)) #for hybrid
#######################################################

DOE Parameter Ranges

We investigate engineering ranges for miles per kilowatt hour (mpkwh) and miles per gallon (mpg).

#######DOE Parameters for mpkwh and mpg#################
#Set mpkWh
mpkwh=rep(0,3)
mpkwh[1]=3
mpkwh[2]=4
mpkwh[3]=5

#Set mpg
mpg=rep(0,3)
mpg[1]=40
mpg[2]=50
mpg[3]=60
#######################################################

Read in the Forecasted Data and Convert to Daily Costs

Retrieve the monthly forecasts for cost per kilowatt hour (cpkWh) and cost per gallon (cpg) and assign them to the days.

#######################################################
f_kw=read.csv("R Code/f_kWh.csv")
f_gas=read.csv("R Code/f_gas.csv")
f_kw=f_kw$x
f_gas=f_gas$x

#Convert Montly Cost Data to Daily Cents / kWh
cpkwh=rep(0,8*365*it)
cpkwh[1:31]=f_kw[1] #Jan
cpkwh[32:59]=f_kw[2] #Feb
cpkwh[60:90]=f_kw[3] #Mar
cpkwh[91:120]=f_kw[4] #Apr
cpkwh[121:151]=f_kw[5] #May
cpkwh[152:181]=f_kw[6] #Jun
cpkwh[182:212]=f_kw[7] #Jul
cpkwh[213:243]=f_kw[8] #Aug
cpkwh[244:273]=f_kw[9] #Sep
cpkwh[274:304]=f_kw[10] #Oct
cpkwh[305:334]=f_kw[11] #Nov
cpkwh[335:365]=f_kw[12] #Dec
cpkwh=rep(cpkwh[1:365], 8*it) #generate 8 years and 16 iterations of daily cost data

#Convert Montly Cost Data to Daily Cents / Gallon
cpgal=rep(0,8*365*it)
cpgal[1:31]=f_gas[1] #Jan
cpgal[32:59]=f_gas[2] #Feb
cpgal[60:90]=f_gas[3] #Mar
cpgal[91:120]=f_gas[4] #Apr
cpgal[121:151]=f_gas[5] #May
cpgal[152:181]=f_gas[6] #Jun
cpgal[182:212]=f_gas[7] #Jul
cpgal[213:243]=f_gas[8] #Aug
cpgal[244:273]=f_gas[9] #Sep
cpgal[274:304]=f_gas[10] #Oct
cpgal[305:334]=f_gas[11] #Nov
cpgal[335:365]=f_gas[12] #Dec
cpgal=rep(cpgal[1:365], 8*it) #generate 8 years and 16 iterations of daily cost data
#######################################################

Estimate the Daily Driving Distance

Generate driving distance from an exponential with a mean of 37 miles per day (DoT).

#######################################################
n=365*8*it  #generate the number of total days in the simulation based on iterations
miles=rexp(n,1/37)  #generate n iterates of a random exponential
#######################################################

Populate the Results

The results are generated without looping.

#######################################################
results[1,]=(1/mpkwh[1])*cpkwh[1:n]*miles[1:n]
results[2,]=(1/mpkwh[2])*cpkwh[1:n]*miles[1:n]
results[3,]=(1/mpkwh[3])*cpkwh[1:n]*miles[1:n]
rownames(results)=c("3 mpkWh","4 mpkWh","5 mpkWh")

results2[1,]=(1/mpg[1])*cpgal[1:n]*miles[1:n]
results2[2,]=(1/mpg[2])*cpgal[1:n]*miles[1:n]
results2[3,]=(1/mpg[3])*cpgal[1:n]*miles[1:n]
rownames(results2)=c("40 mpg","50 mpg","60 mpg")
#######################################################

#######################################################
#Evaluate the standard error
esd=apply(results,1,sd)/sqrt(length(results[1,]))
gsd=apply(results2,1,sd)/sqrt(length(results2[1,]))
esd
##   3 mpkWh   4 mpkWh   5 mpkWh 
## 0.6018337 0.4513753 0.3611002
gsd
##    40 mpg    50 mpg    60 mpg 
## 0.9340495 0.7472396 0.6226997
#######################################################

Validate and Write the Results

Validate output vs. input and write results.

#######################################################
#Validation
mean(cpkwh) #~13 cents
## [1] 13.09988
mean(miles) #~37 miles driven
## [1] 37.19744
mean(1/mpkwh) 
## [1] 0.2611111
mean(cpkwh)*mean(miles)*mean(1/mpkwh)  #in cents
## [1] 127.2347
mean(cpgal)  #in cents
## [1] 270.9223
mean(1/mpg)
## [1] 0.02055556
mean(cpgal)*mean(miles)*mean(1/mpg)  #in cents
## [1] 207.151
#######################################################


################Write the results#######################
#write.csv(t(results), "R Code/electric.csv")
#write.csv(t(results2), "R Code/gas.csv")
#######################################################