setwd('E:\\c\\gep') train = read.table("https://s3-us-west-2.amazonaws.com/geptestdata/data/output_train.txt", header=TRUE, sep=" ") valid = read.table("https://s3-us-west-2.amazonaws.com/geptestdata/data/output_valid.txt", header=TRUE, sep=" ") mtrain = lm(target~predict, data=train) summary(mtrain) mvalid = lm(target~predict, data=valid) summary(mvalid) jpeg('E:\\c\\gep\\TrainFit_elecdemand.jpg', quality = 200, bg = "white", pointsize = 12, width = 480, height = 480, units = "px") plot(train$predict, train$target, xlab="Fitted", ylab="Actual", main="Training (14 days)", col="blue", cex=1.5, cex.main=1.8, cex.lab=1.5,cex.axis=1.5) abline(mtrain,lwd=2) text(3300,4500,"R-Square = 0.948", cex=1.5) dev.off() jpeg('E:\\c\\gep\\ValidFit_elecdemand.jpg', quality = 200, bg = "white", pointsize = 12, width = 480, height =480, units = "px") plot(valid$predict, valid$target, xlab="Predicted", ylab="Actual", main="Validation (7 days)", col="red", cex=1.5, cex.main=1.8, cex.lab=1.5,cex.axis=1.5) abline(mtrain,lwd=2) text(3500,5000,"R-Square = 0.930", cex=1.5) dev.off() MAPE = mean(abs(valid$target-valid$predict)/valid$target) MAE = mean((valid$target-valid$predict)/valid$target) jpeg('E:\\c\\gep\\ValidTS_elecdemand.jpg', quality = 200, bg = "white", pointsize = 12, width = 600, height =480, units = "px") plot(valid$target, col="red", lty=2, pch=16, xlab="Hour", ylab="Load (MW)", main="Prediction (7 days)", cex=1, cex.main=1.8, cex.lab=1.5,cex.axis=1.5) lines(valid$target, lty=2, col="red") lines(valid$predict, lty=1, col="blue") legend(120,5300, c("Actual", "Predicted"), col = c("red", "blue"), lty = c(2, 1), pch=c(16,NA),merge = TRUE,cex=1.5) text(70,5200,paste("MAPE = ",round(MAPE*100,digits=2),"%"),cex=1.5) text(70,5050,paste("MPE = ",round(MPE*100,digits=2),"%"),cex=1.5) dev.off() # Arima model library(forecast) all = rbind(train,valid) fc_log = data.frame(fc=double(0)) ts_all = ts(all$target,start=c(1,1),frequency=24) for(day in seq(0,6)){ ts_train = window(ts_all,start=c(1+day,1),end=c(14+day,24)) fit = arima(ts_train,order=c(3,0,4),seasonal=list(order=c(0,0,2))) fc = predict(fit,n.ahead=24) df = data.frame(fc=as.numeric(fc$pred)) fc_log = rbind(fc_log,df) } tsMAPE = mean(abs(valid$target-fc_log$fc)/valid$target) tsMPE = mean((valid$target-fc_log$fc)/valid$target) jpeg('E:\\c\\gep\\ValidArima_elecdemand.jpg', quality = 200, bg = "white", pointsize = 12, width = 600, height =480, units = "px") plot(valid$target, col="red", lty=2, pch=16, xlab="Hour", ylab="Load (MW)", main="ARIMA Model Day-ahead Forecast", cex=1, cex.main=1.8, cex.lab=1.5,cex.axis=1.5) lines(valid$target, lty=2, col="red") lines(fc_log$fc, lty=1, col="blue") legend(120,5300, c("Actual", "Predicted"), col = c("red", "blue"), lty = c(2, 1), pch=c(16,NA),merge = TRUE,cex=1.5) text(70,5200,paste("MAPE = ",round(tsMAPE*100,digits=2),"%"),cex=1.5) text(70,5050,paste("MPE = ",round(tsMPE*100,digits=2),"%"),cex=1.5) dev.off()