## -----------------------------------------------------------------------------------------------------------------
# Load the package
library(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
# Generate time of day in a forecast matrix
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)")
# Optimize the parameters, only on two horizons
kseqopt <- c(3,18)
rls_optim(model, D, kseqopt)


## -----------------------------------------------------------------------------------------------------------------
# Forecast for all horizons
model$kseq <- 1:36
# Fit with RLS
fit1 <- rls_fit(model$prm, model, D)
# See the summary of the fit
summary(fit1)


## ----output.lines=10----------------------------------------------------------------------------------------------
# Add a diurnal curve using fourier series
model$add_inputs(mu_tday = "fs(tday/24, nharmonics=4)")
# Optimize the parameters
rls_optim(model, D, kseq=kseqopt)


## -----------------------------------------------------------------------------------------------------------------
# Fit with RLS
fit2 <- rls_fit(model$prm, model, D)
# Check the fit
summary(fit2)


## -----------------------------------------------------------------------------------------------------------------
# Keep the forecasts from each model by just inserting them in the data.list
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))


## -----------------------------------------------------------------------------------------------------------------
# The simple persistence (forecast for same horizons as the model)
D$YhatP <- persistence(D$y, model$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, model$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


## -----------------------------------------------------------------------------------------------------------------
# Find all complete cases for all forecasts and horizons
ok <- complete_cases(D[nms])


## -----------------------------------------------------------------------------------------------------------------
sum(ok)
length(ok)


## -----------------------------------------------------------------------------------------------------------------
D$Yhat1[1:11, 1:10]


## -----------------------------------------------------------------------------------------------------------------
D$y[59:72]
D$YhatP[59:72, 1]


## -----------------------------------------------------------------------------------------------------------------
ok <- ok & D$scoreperiod


## -----------------------------------------------------------------------------------------------------------------
sum(ok)
sum(D$scoreperiod)


## -----------------------------------------------------------------------------------------------------------------
# The score as a function of the horizon
R <- residuals(D$Yhat1, D$y)
score(R, ok & D$scoreperiod)


## -----------------------------------------------------------------------------------------------------------------
# Only complete cases are used per default
score(R, D$scoreperiod) == score(R, ok & D$scoreperiod)


## -----------------------------------------------------------------------------------------------------------------
# The score as a function of the horizon
score(R, usecomplete=FALSE) == score(R)


## -----------------------------------------------------------------------------------------------------------------
RMSE <- score(residuals(D[nms], D$y), D$scoreperiod)


## ----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----------------------------------------------------------------------------------------
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("topleft", nms, lty=1, col=1:length(nms))


## ----plottrain----------------------------------------------------------------------------------------------------
D$trainperiod <- in_range(D$t[7*24]-1, D$t, "2011-02-01")    
plot(D$t, D$trainperiod)


## ----output.lines=10----------------------------------------------------------------------------------------------
# Optimize the parameters
rls_optim(model, subset(D,D$trainperiod), kseqopt)


## -----------------------------------------------------------------------------------------------------------------
# Fit with RLS
fit <- rls_fit(model$prm, model, D)


## ----scorefit-----------------------------------------------------------------------------------------------------
score(fit$Yhat, !D$trainperiod)


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


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


## ----plotacf, fig.height=figwidth, echo=-1------------------------------------------------------------------------
par(mfrow=c(2,2))
acf(D$Residuals$h1, na.action=na.pass)
ccf(lagvec(D$Ta$k1,1), D$Residuals$h1, na.action=na.pass)
ccf(lagvec(D$I$k1,1), D$Residuals$h1, na.action=na.pass)
ccf(lagvec(D$I$k1,1), D$Residuals$h1, na.action=na.pass)


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


## -----------------------------------------------------------------------------------------------------------------
par(mfrow=c(1,2))
hist(D$Residuals$h1)
qqnorm(D$Residuals$h1)
qqline(D$Residuals$h1)


## -----------------------------------------------------------------------------------------------------------------
boxplot(D$Residuals$h1 ~ D$tday$k0)