## ---- 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)