## ---- message=FALSE-----------------------------------------------------------
library(onlineforecast)
D <- Dbuilding
# Don't bother with test period here
D$scoreperiod <- in_range("2010-12-20", D$t)
D$tday <- make_tday(D$t, 1:36)
# The heat load forecast model
model <- forecastmodel$new()
model$output = "heatload"
model$add_inputs(Ta = "lp(Ta, a1=0.9)", 
                 I = "lp(I, a1=0.7)", 
                 mu_tday = "fs(tday/24, nharmonics=10)",
                 mu = "one()")
model$add_regprm("rls_prm(lambda=0.9)")
model$add_prmbounds(Ta__a1 = c(0.8, 0.9, 0.9999),
                    I__a1 =  c(0.4, 0.8, 0.9999),
                    lambda = c(0.9, 0.99, 0.9999))
rls_optim(model, D, kseq=c(1,18))
model$kseq <- 1
val <- rls_fit(model$prm, model, D, returnanalysis = TRUE)
D$heatloadHat <- val$Yhat


## -----------------------------------------------------------------------------
# Generate the one-step residuals
residualsk1 <- D$heatload - lagvec(D$heatloadHat$k1, 1)
# or the residual function give the same
tmp <- residuals(D$D$heatloadHat, D$heatload)
unique(tmp$k1 - residualsk1)

# Plot to see 
plot(residualsk1)
acf(residualsk1, na.action = na.pass)


## ---- message=FALSE-----------------------------------------------------------
# The residuals (they are k0, i.e. the observed one-step error at time t)
# Put them in the data
D$residuals <- residualsk1

# Setup the model
model2 <- forecastmodel$new()
model2$output = "residuals"
model2$add_inputs(AR = "AR(0)")
model2$add_regprm("rls_prm(lambda=0.9)")
model2$add_prmbounds(lambda = c(0.9, 0.99, 0.9999))

# Optimize the parameters
rls_optim(model2, D, kseq=1)


## -----------------------------------------------------------------------------
model2$kseq <- 1
val2 <- rls_fit(model2$prm, model2, D, returnanalysis = TRUE)


## ---- fig.height=3.5----------------------------------------------------------
D$residualsHat <- val2$Yhat
plot_ts(D, c("^residuals|^residualsHat"), c("2011-01-01","2011-02-01"), kseq = 1)


## -----------------------------------------------------------------------------
# The RMSE with no model
rmse(D$residuals[D$scoreperiod])
# RMSE on the residuals from the error model
rmse(residuals(D$residualsHat, D$residuals)$h1[D$scoreperiod])


## -----------------------------------------------------------------------------
D$heatloadHat2 <- D$heatloadHat + D$residualsHat

plot_ts(D, c("^heatload$|^heatloadHat"), c("2011-01-20","2011-02-01"), kseq = 1)

rmse(residuals(D$heatloadHat2, D$heatload)$h1[D$scoreperiod])