# ElecLoad contains hourly loads and other data for 2005 and 2006 (=2*365*24 entries):
# 1. Electricity load in MW
# 2. day of weak: sunday=0, monday=1, etc
# 3. Hour of the day 0 -> 23
# 4. Public Holiday: 1 if Public Holiday, 0 otherwise
# 5. Scool vacation: 1 if no scool
# 6. Temperature in °F

# Create the weak matrix = dumy variables for the weakdays
weakmatrix<-model.matrix(~as.factor(ElecLoad[,2]))
#Remove intercept
weakmatrix<-weakmatrix[,-1]

#Generate FullTable
FullTable<-cbind(load=ElecLoad[,1], weakmatrix, ElecLoad[,4],
                 ElecLoad[,3],ElecLoad[,5],ElecLoad[,5]^2, ElecLoad[,6])
colnames(FullTable)<-c("Load","mon","tue","wed","thu","fri","sat",
                       "ScoolHol","PubHol","Temp","Temp2","Hour")

#Create the xreg = substed for a specific hour of the day (column 12 = Hour)
xreg<-subset(FullTable[,2:11], FullTable[,12] == 7)

#Create the Load time serie, also a subset of the full table
LoadTs<-ts(subset(FullTable[,1], FullTable[,12] == 7),start=1,frequency=1)

#Launch of auto.arima
ArimaLoad<-auto.arima(LoadTs, xreg=xreg, lambda=0)

plot(forecast(ArimaLoad,xreg=xreg), include=0)

plot(fitted(ArimaLoad))

#Create the frame.data
Load<-subset(FullTable[,1], FullTable[,12] == 7)
FullData<-cbind(LogLoad=log(Load), xreg)
FrameData<-data.frame(FullData)

# multilinear regression
mlin<-lm(LogLoad ~ mon+tue+wed+thu+fri+sat+ScoolHol+PubHol+Temp+`Temp2`, FrameData)
plot(exp(mlin$model$LogLoad), type="l",col="blue")
lines(exp(fitted(mlin)), col="red")