### Clean up the workspace rm(list = ls()) ## additional libraries library(devtools) library(splines) library(quantreg) ## Load the onlineforecast package - this is the development version load_all(".") # need to be under the package directory # Load the data D <- Dsolar names(D) # Short Data Anlysis tstart <- as.POSIXct("2009-01-01", tz = "GMT") tend <- as.POSIXct("2011-01-01", tz = "GMT") D <- subset(D, in_range(tstart,D$t,tend)) plotmat <- matrix(1:3, nrow = 3, byrow=TRUE) layout(plotmat, widths=c(1,1,1), heights=c(1,1,1)) par(mar=c(0, 3.5, 0.5, 1), lwd = 2, cex.lab = 1.5, cex.axis = 1.5) plot(D$t, D$Ps, type = "l", bty = "l", xlab = "", ylab = "", axes = F) axis(2, seq(0,ceiling(max(D$Ps)+0.2), by = 0.5), las = 1) plot(D$t, lagvec(D$I$k1,1), type = "l", col = "gold", bty = "l", xlab = "", ylab = "", axes = F) lines(D$t, D$Iobs, col = rgb(0, 0, 255, max = 255, alpha = 125)) axis(2, seq(0,ceiling(max(D$I$k1)+0.2), by = 100), las = 1) plot(D$t, lagvec(D$Ta$k1,1), type = "l", col = rgb(240, 37, 37, max = 255, alpha = 150), bty = "l", xlab = "", ylab = "", axes = F, ylim = c(min(D$Ta$k1) - 4,max(D$Ta$k1))) axis(2, seq(floor(min(D$Ta$k1)),ceiling(max(D$Ta$k1)+0.2), by = 10), las = 1) axis.POSIXct(side = 1, x = D$t, xaxt = "s", at = seq(D$t[1], D$t[length(D$t)], by = "1 month"), format = "%Y-%m-%d", mgp = c(4,-1,-2)) # Make a training set of 3 months, and a then test set with 1 month tstart <- "2009-06-01" ttest <- "2009-03-31" tend <- "2009-07-01" # Cut only the necessary period D$scoreperiod <- in_range("2009-06-10", D$t) Dtrain <- subset(D, in_range(tstart,D$t,tend)) ## Create the model, where N1 is the initial number of data points - for the cold start model <- qmodel$new(N1 = 50) ## The output is the Ps from the Dsolar data model$output = "Ps" ## The design matrix, here we are using 5 bsplines fixed at 5 and 18 hour of the day - ie when sun is up and down. model$add_inputs(I = "bspline(tday, df=5, Boundary.knots=c(5,18)) %**% I") ## Add the regularization parameter, lambda, and the bounds model$add_regprm("rls_prm(lambda=0.9)") model$add_prmbounds(lambda = c(0.99, 0.999, 0.9999)) ### Setup of the booking keeping matrix - ############################################################## ####### THIS SHOULD BE DONE IN quantile_fit.R ################ ############################################################## ## First 5 columns are the input Matrix model$IX <- 0:4 ## The number of predictors, we could just use the length of the input matrix or the IX model$K <- 5 ## The output column is the 6th one model$Iy <- 5 ## Just quick fun to see the transformation ## maybe have splines that sum to one always.. model$datatr <- model$transform_data(Dtrain) plot(Dtrain$I$k1[1:100], type = "l") lines(model$datatr$I.bs1$k1, col = "red") lines(model$datatr$I.bs2$k1, col = "red") lines(model$datatr$I.bs3$k1, col = "red") lines(model$datatr$I.bs4$k1, col = "red") lines(model$datatr$I.bs5$k1, col = "red") ## Select which quantiles we want to optimise and fit model$tau <- c(0.05,0.2, 0.5, 0.8,0.95) ## Select which prediction horizon model$kseq <- c(1,4,8,12,24) model$kseq <- 1:24 #opt_model <- quantile_optim(model = model, data = Dtrain) PAR <- c("lambda" = 0.999) quantile_fit(prm = PAR, model = model, data = Dtrain, q = model$tau[1]) PAR <- c("lambda" = 0.999) quantile_fit(prm = PAR, model = model, data = Dtrain, q = model$tau[2]) PAR <- c("lambda" = 0.999) quantile_fit(prm = PAR, model = model, data = Dtrain, q = model$tau[3]) PAR <- c("lambda" = 0.999) quantile_fit(prm = PAR, model = model, data = Dtrain, q = model$tau[4]) PAR <- c("lambda" = 0.999) quantile_fit(prm = PAR, model = model, data = Dtrain, q = model$tau[5]) Pred_model <- quantile_predict(model = model, datatr = model$datatr) plot(Dtrain$Ps, type = "l") lines(Pred_model$q0.05[,"k24"], type = "l", col = "green") lines(Pred_model$q0.5[,"k24"], type = "l", col = "blue") #lines(Pred_model$q0.5[,"k1"], type = "l", col = "green") lines(Pred_model$q0.95[,"k24"], type = "l", col = "red") str(Pred_model) k <- 24 colSeq <- grey(seq(0.9,0.1,len=3)) #colSeq <- colorJet(nPolygons) require( RColorBrewer ) colSeq <- brewer.pal(11 , "Spectral" ) for(i in k:(length(Dtrain$t)-1)){ par(mfrow = c(2,1)) par(mar = c(2,4,2,2)) # bottom left top right with(Pred_model, { plot(Dtrain$t[(i-10):(i+k)], Dtrain$Ps[(i-10):(i+k)], bty = "l", lwd = 2, col = "black", pch=19, cex=0.5, axes = FALSE, xaxt = "n", type="n", ylim = range(Dtrain$Ps[(i-10):(i+k)], q0.5[i,], q0.95[i,], q0.05[i,], na.rm = T), main = "TimeAdaptive with Weights", xlab = "Time", ylab = "Ps") axis(2) polygon(c(Dtrain$t[(i+1):(i+k)], rev(Dtrain$t[(i+1):(i+k)])), c(Pred_model$q0.05[i,],rev(Pred_model$q0.2[i,])), col= colSeq[3], border = NA) polygon(c(Dtrain$t[(i+1):(i+k)], rev(Dtrain$t[(i+1):(i+k)])), c(Pred_model$q0.8[i,],rev(Pred_model$q0.95[i,])), col= colSeq[3], border = NA) polygon(c(Dtrain$t[(i+1):(i+k)], rev(Dtrain$t[(i+1):(i+k)])), c(Pred_model$q0.2[i,],rev(Pred_model$q0.5[i,])), col= colSeq[4], border = NA) polygon(c(Dtrain$t[(i+1):(i+k)], rev(Dtrain$t[(i+1):(i+k)])), c(Pred_model$q0.5[i,],rev(Pred_model$q0.8[i,])), col= colSeq[4], border = NA) lines(Dtrain$t[(i-10):(i+k)], Dtrain$Ps[(i-10):(i+k)], col = "black", type = "b", lwd = 2) lines(Dtrain$t[(i+1):(i+k)], q0.5[i,], type = "b", col = "grey", lwd = 2) axis.POSIXct(side = 1, x = Dtrain$t[(i-10):(i+k)], at = seq(from = Dtrain$t[(i-10):(i+k)][1], to = Dtrain$t[(i-10):(i+k)][length(Dtrain$t[(i-10):(i+k)])], by = "1 hour"), format = "%Y/%m/%d %H:%M \n %a", las = 1, cex.axis = 1, srt = 45) #lines(Pred_model$q0.05[i,], type = "b", col = "blue") #lines(Pred_model$q0.95[i,], type = "b", col = "blue") plot(Dtrain$t[(i-10):(i+k)], c(rep(NA,11), as.numeric(Dtrain$I[i, 1:k])), type = "b", col = "steelblue", axes = FALSE, xlab = "Time", ylab = "Temp", lwd = 2, ylim = range(Dtrain$I[i, 1:k], Dtrain$Iobs[(i+1):(i+k)])) #lines(Dtrain$t[i], Dtrain$Iobs[i+1,1], type = "b", col = "red", axes = FALSE) lines(Dtrain$t[(i+1):(i+k)], Dtrain$Iobs[(i+1):(i+k)], type = "b", col = "red", lwd = 2) axis(2) abline(h=0,v=Dtrain$t[i],lty=2, col = "lightgrey", lwd = 2) legend("topleft", legend=c("I Obs", "I Pred"), col=c("red", "steelblue"), lty=1, cex=1, lwd = 2) i axis.POSIXct(side = 1, x = Dtrain$t[(i-10):(i+k)], at = seq(from = Dtrain$t[(i-10):(i+k)][1], to = Dtrain$t[(i-10):(i+k)][length(Dtrain$t[(i-10):(i+k)])], by = "1 hour"), format = "%Y/%m/%d %H:%M \n %a", las = 1, cex.axis = 1, srt = 45) }) } library(animation) k <- 24 colSeq <- grey(seq(0.9,0.1,len=3)) #colSeq <- colorJet(nPolygons) require( RColorBrewer ) colSeq <- brewer.pal(11 , "Spectral" ) saveHTML({ for(i in k:(60)){ par(mfrow = c(2,1)) par(mar = c(2,4,2,2)) # bottom left top right with(Pred_model, { plot(Dtrain$t[(i-10):(i+k)], Dtrain$Ps[(i-10):(i+k)], bty = "l", lwd = 2, col = "black", pch=19, cex=0.5, axes = FALSE, xaxt = "n", type="n", ylim = range(Dtrain$Ps[(i-10):(i+k)], q0.5[i,], q0.95[i,], q0.05[i,], na.rm = T), main = "TimeAdaptive with Weights", xlab = "Time", ylab = "Ps") axis(2) polygon(c(Dtrain$t[(i+1):(i+k)], rev(Dtrain$t[(i+1):(i+k)])), c(Pred_model$q0.05[i,],rev(Pred_model$q0.2[i,])), col= colSeq[3], border = NA) polygon(c(Dtrain$t[(i+1):(i+k)], rev(Dtrain$t[(i+1):(i+k)])), c(Pred_model$q0.8[i,],rev(Pred_model$q0.95[i,])), col= colSeq[3], border = NA) polygon(c(Dtrain$t[(i+1):(i+k)], rev(Dtrain$t[(i+1):(i+k)])), c(Pred_model$q0.2[i,],rev(Pred_model$q0.5[i,])), col= colSeq[4], border = NA) polygon(c(Dtrain$t[(i+1):(i+k)], rev(Dtrain$t[(i+1):(i+k)])), c(Pred_model$q0.5[i,],rev(Pred_model$q0.8[i,])), col= colSeq[4], border = NA) lines(Dtrain$t[(i-10):(i+k)], Dtrain$Ps[(i-10):(i+k)], col = "black", type = "b", lwd = 2) lines(Dtrain$t[(i+1):(i+k)], q0.5[i,], type = "b", col = "grey", lwd = 2) axis.POSIXct(side = 1, x = Dtrain$t[(i-10):(i+k)], at = seq(from = Dtrain$t[(i-10):(i+k)][1], to = Dtrain$t[(i-10):(i+k)][length(Dtrain$t[(i-10):(i+k)])], by = "1 hour"), format = "%Y/%m/%d %H:%M \n %a", las = 1, cex.axis = 1, srt = 45) #lines(Pred_model$q0.05[i,], type = "b", col = "blue") #lines(Pred_model$q0.95[i,], type = "b", col = "blue") plot(Dtrain$t[(i-10):(i+k)], c(rep(NA,11), as.numeric(Dtrain$I[i, 1:k])), type = "b", col = "steelblue", axes = FALSE, xlab = "Time", ylab = "Temp", lwd = 2, ylim = range(Dtrain$I[i, 1:k], Dtrain$Iobs[(i+1):(i+k)])) #lines(Dtrain$t[i], Dtrain$Iobs[i+1,1], type = "b", col = "red", axes = FALSE) lines(Dtrain$t[(i+1):(i+k)], Dtrain$Iobs[(i+1):(i+k)], type = "b", col = "red", lwd = 2) axis(2) abline(h=0,v=Dtrain$t[i],lty=2, col = "lightgrey", lwd = 2) legend("topleft", legend=c("Ta Obs", "Ta Pred"), col=c("red", "steelblue"), lty=1, cex=1, lwd = 2) i axis.POSIXct(side = 1, x = Dtrain$t[(i-10):(i+k)], at = seq(from = Dtrain$t[(i-10):(i+k)][1], to = Dtrain$t[(i-10):(i+k)][length(Dtrain$t[(i-10):(i+k)])], by = "1 hour"), format = "%Y/%m/%d %H:%M \n %a", las = 1, cex.axis = 1, srt = 45) }) }}, htmlfile = "misc-R/quantile/model.html", ani.height = 800, ani.width = 1600)