## ---- echo=1:2, purl=1:2------------------------------------------------------ # Load the package library(onlineforecast) #library(devtools) #load_all(as.package("../../onlineforecast")) ## ----------------------------------------------------------------------------- # Keep the data in D to simplify notation D <- Dbuilding # Keep the model output in y (just easier code later) D$y <- D$heatload # D$tday <- make_tday(D$t, 0:36) ## ----------------------------------------------------------------------------- # Set the score period D$scoreperiod <- in_range("2010-12-22", D$t) ## ---- output.lines=10--------------------------------------------------------- # Define a new model with low-pass filtering of the Ta input model <- forecastmodel$new() model$output = "y" model$add_inputs(Ta = "lp(Ta, a1=0.9)", I = "lp(I, a1=0.9)", mu = "one()") model$add_prmbounds(Ta__a1 = c(0.8, 0.9, 0.99), I__a1 = c(0.6, 0.9, 0.99), lambda = c(0.9, 0.99, 0.9999)) model$add_regprm("rls_prm(lambda=0.9)") model$kseq <- c(3,18) # Optimize the parameters model$prm <- rls_optim(model, D)$par ## ----------------------------------------------------------------------------- # Fit for all horizons model$kseq <- 1:36 # Fit with RLS fit1 <- rls_fit(model$prm, model, D) # Check the fit summary(fit1) ## ---- output.lines=10--------------------------------------------------------- # Add a diurnal curve using fourier series model$add_inputs(mu_tday = "fs(tday/24, nharmonics=4)") model$kseq <- c(3,18) # Optimize the parameters model$prm <- rls_optim(model, D)$par ## ----------------------------------------------------------------------------- # Fit for all horizons model$kseq <- 1:36 # Fit with RLS fit2 <- rls_fit(model$prm, model, D) # Check the fit summary(fit2) ## ----------------------------------------------------------------------------- # Keep the forecasts from each model D$Yhat1 <- fit1$Yhat D$Yhat2 <- fit2$Yhat ## ---- fig.height=figheight2--------------------------------------------------- # Plot to see the forecasts for the shortest and the longest horizon plot_ts(subset(D,D$scoreperiod), c("^y|^Yhat1","^y|^Yhat2"), kseq = c(1,36)) ## ---- fig.height=figheight2--------------------------------------------------- # Plot to see the forecasts for the shortest and the longest horizon plot_ts(subset(D,which(D$scoreperiod)[1:(14*24)]), c("^y|^Yhat1","^y|^Yhat2"), kseq = c(1,36)) ## ----------------------------------------------------------------------------- # Just keep the horizons kseq <- 1:36 # The simple persistence D$YhatP <- persistence(D$y, kseq) # Plot a few horizons plot_ts(D, c("^y$|YhatP$"), c("2011-01-05","2011-01-10"), kseq=c(1,24,36)) ## ----------------------------------------------------------------------------- D$YhatP[1:4, 1:8] ## ----------------------------------------------------------------------------- # Use the argument perlen to set the period length D$YhatDP <- persistence(D$y, kseq, perlen=24) # Plot a few horizons plot_ts(D, c("^y$|YhatDP$"), c("2011-01-05","2011-01-10"), kseq=c(1,24,36)) ## ----------------------------------------------------------------------------- # Find the forecasts in D nms <- grep("^Yhat", names(D), value=TRUE) nms ## ----------------------------------------------------------------------------- # The non-NA for the first forecast ok <- !is.na(D[[nms[1]]]) # Go through the remaining: all must be non-NA for a point for(nm in nms[-1]){ ok <- ok & !is.na(D[[nm]]) } ok <- as.data.frame(ok) names(ok) <- pst("k",kseq) # Lag to match resiuduals in time ok <- lagdf(ok, "+k") # Only the score period ok <- ok & D$scoreperiod # Finally, the vector with TRUE for all points with no NAs for any forecast ok <- apply(ok, 1, all) ## ----------------------------------------------------------------------------- sum(ok) length(ok) ## ----------------------------------------------------------------------------- # Use the residuals function R <- residuals(D$Yhat1, D$y) # And the score as a function of the horizon score_for_k(R, scoreperiod=ok)$scoreval ## ----------------------------------------------------------------------------- RMSE <- sapply(nms, function(nm){ score_for_k(residuals(D[[nm]],D$y), ok)$scoreval }) ## ---- include=FALSE----------------------------------------------------------- # sapply(kseq, function(k){ # rmse(y - lagdf(YhatDM[ ,pst("k",k)], k)) # # hej det er vilfred jeg er peders søn og jeg elsker min far go jeg god til matematik og jeg elsker også min mor # }) ## ---- fig.height=figheight2--------------------------------------------------- RMSE <- as.data.frame(RMSE) names(RMSE) <- nms plot(0, type="n", xlim=range(kseq), ylim=range(RMSE), xlab="Horizon k", ylab="RMSE (kW)") for(i in 1:length(RMSE)){ points(kseq, RMSE[ ,i], type="b", col=i) } ## ----plottrain---------------------------------------------------------------- D$trainperiod <- in_range(D$t[1]-1, D$t, "2011-02-01") plot(D$t, D$trainperiod) ## ---- output.lines=10--------------------------------------------------------- model$kseq <- c(3,18) # Optimize the parameters model$prm <- rls_optim(model, subset(D,D$trainperiod))$par ## ----------------------------------------------------------------------------- # Fit for all horizons model$kseq <- 1:36 # Fit with RLS fittmp <- rls_fit(model$prm, model, D) ## ----scorefit----------------------------------------------------------------- score_fit(fittmp, !D$trainperiod)$scoreval ## ----plot1, fig.height=figheight5--------------------------------------------- kseq <- c(1,18,36) plot_ts(fit1, kseq=kseq) ## ---- fig.height=figheight5, fig.keep="none"---------------------------------- plot_ts(fit2, kseq=kseq) ## ---- fig.height=figheight5, fig.keep="none"---------------------------------- xlim <- c("2011-01-01","2011-01-14") plot_ts(fit1, xlim=xlim, kseq=kseq) plot_ts(fit2, xlim=xlim, kseq=kseq) ## ----tscoef------------------------------------------------------------------- tmp <- plot_ts(fit2, kseq=kseq, plotit=FALSE) class(tmp) names(tmp) # Residuals plot_ts(tmp, c("^Residuals"), kseq=kseq) # All RLS coefficients nms <- names(fit2$Lfitval$k1) plot_ts(tmp, pattern=nms, kseq=kseq) ## ---- fig.height=figheight3--------------------------------------------------- plot_ts(tmp, pattern=c("^y$|^Yhat",nms), c("2011-02-06","2011-02-10"), kseq=kseq) ## ----rlscoefficients---------------------------------------------------------- str(fit1$Lfitval$k1) ## ----plotpairs, fig.height=figwidth------------------------------------------- kseq <- c(1,36) D$Residuals <- residuals(fit2)[ ,pst("h",kseq)] D$hour <- aslt(D$t)$hour pairs(D, subset=D$scoreperiod, pattern="Residuals|Ta|I|hour|^t$", kseq=kseq) ## ----------------------------------------------------------------------------- par(mfrow=c(1,2)) hist(D$Residuals$h1) qqnorm(D$Residuals$h1) qqline(D$Residuals$h1) ## ----------------------------------------------------------------------------- boxplot(D$Residuals$h1 ~ D$tday$k0)