## ---- message=FALSE----------------------------------------------------------- # For storing data and cache files cachedir <- tempdir() # Load the package and data library(onlineforecast) path <- paste0(cachedir,"/DdhwMultiFamily.RDS") download.file("https://onlineforecasting.org/examples/data/DdhwMultiFamily.RDS", path) D <- as.data.list(readRDS(path)) # See the variables str(D) ## ----------------------------------------------------------------------------- plot_ts(D) ## ----------------------------------------------------------------------------- xlim <- c("2019-01-14","2019-01-21") plot_ts(D, xlim=xlim) ## ----------------------------------------------------------------------------- D$tday <- make_tday(D$t, 0:36) head(D$tday) ## ----------------------------------------------------------------------------- boxplot(D$dhw ~ D$tday$k0) ## ----------------------------------------------------------------------------- D$scoreperiod <- in_range(D$t[1]+14*24*3600, D$t) ## ----------------------------------------------------------------------------- # New model model <- forecastmodel$new() # Clip the output to be inside this range model$output <- "dhw" model$outputrange <- c(0,Inf) # The forgetting factor model$add_regprm("rls_prm(lambda=0.9)") model$add_prmbounds(lambda = c(0.9, 0.999, 0.999999)) model$kseq <- 1:36 model$kseqopt <- c(1,18) ## ----------------------------------------------------------------------------- model1 <- model$clone_deep() model1$add_inputs(mu_tday = "fs(tday/24, nharmonics=8)", mu = "one()") ## ----------------------------------------------------------------------------- prm <- list(mu_tday__nharmonics = c(min=2, init=5, max=12)) Lstep1 <- step_optim(model1, D, prm, method="Brent", keepinputs = TRUE, cachedir=cachedir) # Keep it model1 <- Lstep1$final$model # The trace to see how the level of improvement getse(Lstep1,"score") ## ---- fig.height=figheight2--------------------------------------------------- fit1 <- rls_fit(model1$prm, model1, D) D$Yhat1 <- fit1$Yhat plot_ts(D, c("^dhw","^Yhat1"), kseq = c(1,18)) ## ----------------------------------------------------------------------------- plot_ts(D, c("^dhw|^Yhat1"), xlim, kseq=model$kseqopt) ## ----------------------------------------------------------------------------- # wday == 0 is Sunday and 6 is Saturday workday <- as.numeric(aslt(D$t)$wday %in% 1:5) D$workday <- lagdf(make_input(workday,0:36), 0:-36) head(D$workday) # One for weekend D$weekend <- 1 - D$workday ## ----------------------------------------------------------------------------- plot_ts(D, c("workday","weekend")) ## ---- fig.height=figheight3--------------------------------------------------- Residuals <- D$scoreperiod * residuals(fit1) par(mfrow=c(3,1)) boxplot(Residuals$h1 ~ D$tday$k0) boxplot(Residuals$h1 ~ D$tday$k0, subset=D$workday$k0==1) boxplot(Residuals$h1 ~ D$tday$k0, subset=D$weekend$k0==1) ## ----------------------------------------------------------------------------- model2 <- model$clone_deep() model2$add_inputs(mu_tday_workday = "fs(tday/24, nharmonics=10) %**% workday", mu_tday_weekend = "fs(tday/24, nharmonics=10) %**% weekend", mu_workday = "workday", mu_weekend = "weekend") ## ----------------------------------------------------------------------------- # Now optimize prm <- list(mu_tday_workday__nharmonics = c(min=4, init=12, max=12), mu_tday_weekend__nharmonics = c(min=4, init=8, max=12)) Lstep2 <- step_optim(model2, D, prm, method="Brent", keepinputs = TRUE, cachedir=cachedir) model2 <- Lstep2$final$model # Fit to get the forecasts fit2 <- rls_fit(model2$prm, model2, D) D$Yhat2 <- fit2$Yhat ## ----------------------------------------------------------------------------- plot_ts(D, c("dhw|Yhat2"), xlim, kseq=model$kseqopt) ## ---- fig.height=figheight2--------------------------------------------------- plot_ts(D, c("dhw","Yhat2"), xlim=c("2019-02-01","2019-03-01"), kseq=model$kseqopt) ## ----------------------------------------------------------------------------- RMSE <- score(residuals(D[c("Yhat1","Yhat2")], D$dhw), D$scoreperiod) plot(0, type="n", xlim=range(model$kseq), ylim=range(RMSE), xlab="Horizon k", ylab="RMSE (kW)") for(i in 1:ncol(RMSE)){ points(model$kseq, RMSE[ ,i], type="b", col=i) } legend("right", c("Single daily","Workdays and weekend"), lty=1, col=1:2) ## ---- fig.height=figheight3--------------------------------------------------- Residuals <- D$scoreperiod * residuals(fit2) par(mfrow=c(3,1)) boxplot(Residuals$h1 ~ D$tday$k0) boxplot(Residuals$h1 ~ D$tday$k0, subset=D$workday$k0==1) boxplot(Residuals$h1 ~ D$tday$k0, subset=D$weekend$k0==1) ## ---- fig.height=figheight2--------------------------------------------------- nms <- names(fit2$Lfitval$k1) D$coefwork <- fit2$Lfitval$k1[grep("workday",nms)] D$coefweekend <- fit2$Lfitval$k1[grep("weekend",nms)] plot_ts(D, c("coefwork","coefweekend"), xlim=c("2019-02-01","2019-03-01"), legendcex=0.7, legendspace=15) ## ----------------------------------------------------------------------------- model1$output <- "y" D$y <- D$dhw D$y[D$workday$k0 != 1] <- NA prm <- list(mu_tday__nharmonics = c(min=2, init=5, max=12)) m <- step_optim(model1, D, prm, method="Brent", keepinputs = TRUE, cachedir=cachedir)$final$model fitworkday <- rls_fit(m$prm, m, D) ## ----------------------------------------------------------------------------- D$y <- D$dhw D$y[D$weekend$k0 != 1] <- NA m <- step_optim(model1, D, prm, method="Brent", keepinputs = TRUE, cachedir=cachedir)$final$model fit <- rls_fit(m$prm, m, D) ## ---- fig.height=figheight3--------------------------------------------------- # Combine them i <- D$workday[ ,pst("k",fit$model$kseq)] == 1 i[is.na(i)] <- FALSE fit$Yhat[i] <- fitworkday$Yhat[i] # Plot them D$Yhat <- fit$Yhat plot_ts(D, c("dhw","^Yhat$","Yhat2"), xlim=c("2019-02-01","2019-03-01"), kseq=model$kseqopt) ## ----------------------------------------------------------------------------- RMSE <- score(residuals(D[c("Yhat","Yhat2")], D$dhw), D$scoreperiod) plot(0, type="n", xlim=range(model$kseq), ylim=range(RMSE), xlab="Horizon k", ylab="RMSE (kW)") for(i in 1:ncol(RMSE)){ points(model$kseq, RMSE[ ,i], type="b", col=i) } legend("right", c("Two models","One model"), lty=1, col=1:2) ## ---- message=FALSE----------------------------------------------------------- path <- paste0(cachedir,"/Ddhw11Houses.RDS") download.file("https://onlineforecasting.org/examples/data/Ddhw11Houses.RDS", path) D <- as.data.list(readRDS(path)) ## ---- fig.height=figheight3--------------------------------------------------- D <- subset(D, c("2010-02-01","2010-07-01")) plot_ts(D, c(pst("Pdhw",1:11,"$"),"Pdhw$")) ## plot_ts(D, c(pst("Pdhw",1:11,"$"),"Pdhw$"), xlim=c("2010-04-01","2010-04-14")) ## ---- fig.height=figheight3--------------------------------------------------- D$tday <- make_tday(D$t, 0:36) head(D$tday) D$scoreperiod <- in_range(D$t[1]+14*24*3600, D$t) # Take the this one, which is made as the mean of all output <- "Pdhw" ## ---- fig.height=figheight3--------------------------------------------------- # New model model <- forecastmodel$new() model$output <- output # Clip the output to be inside this range model$outputrange <- c(0,Inf) # The forgetting factor model$add_regprm("rls_prm(lambda=0.9)") model$add_prmbounds(lambda = c(0.9, 0.999, 0.999999)) # The horizons model$kseq <- 1:36 model$kseqopt <- c(1,18) # Add inputs model1 <- model$clone_deep() model1$add_inputs(mu_tday = "fs(tday/24, nharmonics=8)", mu = "one()") # Step to optimize prm <- list(mu_tday__nharmonics = c(min=2, init=8, max=12)) model1 <- step_optim(model1, D, prm, method="Brent", keepinputs = TRUE, cachedir=cachedir)$final$model # fit1 <- rls_fit(model1$prm, model1, D) D$Yhat1 <- fit1$Yhat plot_ts(D, c(pst("^",output),"^Yhat1"), kseq=model$kseqopt) # plot_ts(D, c(pst("^",output),"^Yhat1"), xlim=c("2010-04-01","2010-04-14"), kseq=model$kseqopt) # wday == 0 is Sunday and 6 is Saturday workday <- as.numeric(aslt(D$t)$wday %in% 1:5) D$workday <- lagdf(make_input(workday,0:36), 0:-36) head(D$workday) # One for weekend D$weekend <- 1 - D$workday ## ----------------------------------------------------------------------------- Residuals <- D$scoreperiod * residuals(fit1) par(mfrow=c(3,1)) boxplot(Residuals$h1 ~ D$tday$k0) boxplot(Residuals$h1 ~ D$tday$k0, subset=D$workday$k0==1) boxplot(Residuals$h1 ~ D$tday$k0, subset=D$weekend$k0==1) ## ---- fig.height=figheight2--------------------------------------------------- D$y <- D[[output]] model1$output <- "y" D$y[D$workday$k0 != 1] <- NA m <- step_optim(model1, D, prm, method="Brent", keepinputs = TRUE, cachedir=cachedir)$final$model fitworkday <- rls_fit(m$prm, m, D) # Similarly, for the weekends: D$y <- D[[output]] D$y[D$weekend$k0 != 1] <- NA m <- step_optim(model1, D, prm, method="Brent", keepinputs = TRUE, cachedir=cachedir)$final$model fit <- rls_fit(m$prm, m, D) D$Yhat <- fit$Yhat # Combine them i <- D$workday[ ,pst("k",fit$model$kseq)] == 1 i[is.na(i)] <- FALSE D$Yhat[i] <- fitworkday$Yhat[i] ## ----------------------------------------------------------------------------- Residuals <- D$scoreperiod * residuals(fit$Yhat, D[[output]]) plot(Residuals$h24) plot(Residuals$h1) ## ----------------------------------------------------------------------------- RMSE <- score(residuals(D[c("Yhat","Yhat1")], D[[output]]), D$scoreperiod) plot(0, type="n", xlim=range(model$kseq), ylim=range(RMSE), xlab="Horizon k", ylab="RMSE (kW)") for(i in 1:ncol(RMSE)){ points(model$kseq, RMSE[ ,i], type="b", col=i) } legend("right", c("Two models","One model"), lty=1, col=1:2)