##########################################################
# Forecasting electricity consumption in a functional way: 
# R/S-plus commandlines
##########################################################
#
# Entering and organizing dataset 
#################################
CONSELDAT <- as.matrix(read.table("npfda-electricity.dat"))
attributes(CONSELDAT)$dimnames <- NULL
learning <- 1:26
testing <- 27
elec.past.learn<-CONSELDAT[learning,]  # sample of explanatory curves
elec.past.testing<-CONSELDAT[testing,] # The 27th year
s <- 1                                 # forecasting horizon 1
elec.futur.s <- CONSELDAT[2:27,s]      # sample of real responses
##
# Functional nonparametric forecasting
##########################################
# Kernel functional regression forecasting
#-----------------------------------------
result.pred.reg.step.s <- funopare.knn.lcv(elec.futur.s,elec.past.learn,elec.past.testing,5,kind.of.kernel="quadratic",semimetric="pca")
# Kernel functional median forecasting
#-------------------------------------
result.pred.median.step.s <- funopare.quantile.lcv(elec.futur.s,elec.past.learn,elec.past.testing,2,alpha=0.5, Knearest=NULL, kind.of.kernel="quadratic",semimetric="pca")
# Kernel functional mode forecasting
#-----------------------------------
result.pred.mode.step.s <- funopare.mode.lcv(elec.futur.s,elec.past.learn,elec.past.testing,2,Knearest=NULL, kind.of.kernel="quadratic",semimetric="pca")
# Median estimation and 90% prediction band
#------------------------------------------
result.pred.quantiles.s<-funopare.quantile.lcv(elec.futur.s,elec.past.learn,elec.past.testing,2,alpha=c(0.05,0.5,0.95), Knearest=NULL,kind.of.kernel="quadratic",semimetric="pca")
##
# Collecting the forecasting results
####################################
# Forecasted value with regression method
#----------------------------------------
result.pred.reg.step.s$Predicted.values
# Forecasted value with median method
#------------------------------------
result.pred.median.step.s$Predicted.values
# Forecasted value with mode method
#----------------------------------
result.pred.mode.step.s$Predicted.values
#.05, .5 and .95 estimated quantiles
#-----------------------------------
result.pred.quantiles.s$Predicted.values
##
# Forecasting the 28th year
###########################
pred.reg <- 0
pred.median <- 0
pred.mode <- 0
for(s in 1:12){
        elec.futur.s <- CONSELDAT[2:27,s] # forecasting horizon "s"
        result.pred.step.s <- funopare.knn.lcv(elec.futur.s,elec.past.learn,elec.past.testing,5,kind.of.kernel="quadratic",semimetric="pca")
        pred.reg[s] <- result.pred.step.s$Predicted.values
        result.pred.median.step.s <- funopare.quantile.lcv(elec.futur.s,elec.past.learn,elec.past.testing,2,alpha=0.5,Knearest=NULL, kind.of.kernel="quadratic", semimetric="pca")
        pred.median[s] <- result.pred.median.step.s$Predicted.values
        result.pred.mode.step.s <- funopare.mode.lcv(elec.futur.s,elec.past.learn,elec.past.testing,2,Knearest=NULL,kind.of.kernel="quadratic",semimetric="pca")
        pred.mode[s] <- result.pred.mode.step.s$Predicted.values
}
##
# Plotting forecasted values
############################
year28 <- CONSELDAT[28,]
mse.reg <- round(sum((pred.reg-year28)^2)/12,4)
mse.median <- round(sum((pred.median-year28)^2)/12,4)
mse.mode <- round(sum((pred.mode-year28)^2)/12,4)
par(mfrow=c(2,2))
plot(1:12,pred.reg,xlab='28th year', ylab='', main=paste('Regression: MSE=',mse.reg,sep=''), type='l',lty=2,ylim=range(c(pred.reg,year28)))
par(new=T)
plot(1:12,year28,type='l',lty=1,axes=F,xlab='', ylab='')
plot(1:12,pred.median,xlab='28th year', ylab='', main=paste('Median: MSE=',mse.median,sep=''), type='l',lty=2,ylim=range(c(pred.median,year28)))
par(new=T)
plot(1:12,year28,type='l',lty=1,axes=F,xlab='', ylab='', ylim=range(c(pred.median,year28)))
plot(1:12,pred.mode,xlab='28th year', ylab='', main=paste('Mode: MSE=',mse.mode,sep=''), type='l',lty=2,ylim=range(c(pred.mode,year28)))
par(new=T)
plot(1:12,year28,type='l',lty=1,axes=F,xlab='', ylab='', ylim=range(c(pred.mode,year28)))