## ---- echo=1:2, purl=1:2------------------------------------------------------
# Load the package
library(onlineforecast)
#library(devtools)
#load_all(as.package("../../onlineforecast"))


## ---- output.lines=10---------------------------------------------------------
# Keep the data in D to simplify notation
D <- Dbuilding
# Set the score period 
D$scoreperiod <- in_range("2010-12-20", D$t)
# Set the training period
D$trainperiod <- in_range(D$t[1], D$t, "2011-02-01")
# Define a new model with low-pass filtering of the Ta input
model <- forecastmodel$new()
model$output = "heatload"
model$add_inputs(Ta = "lp(Ta, a1=0.9)",
                 mu = "one()")
model$add_regprm("rls_prm(lambda=0.9)")
model$add_prmbounds(Ta__a1 = c(0.5, 0.9, 0.9999),
                    lambda = c(0.9, 0.99, 0.9999))
model$kseq <- 1:36
# Optimize the parameters on the train period (i.e. until 2011-02-01)
rls_optim(model, subset(D,D$trainperiod), kseq = c(3,18))


## -----------------------------------------------------------------------------
# Keep a sequence with these points
iseq <- which(in_range("2010-12-15",D$t,"2011-01-01"))
# Fit the model
rls_fit(model$prm, model, subset(D, iseq))


## -----------------------------------------------------------------------------
# The data of a fit is saved in this list
str(model$Lfits[1:2])


## -----------------------------------------------------------------------------
# Next time step
(i <- iseq[length(iseq)] + 1)
# The new data for this time step
Dnew <- subset(D, i)


## -----------------------------------------------------------------------------
# Transform the new data
DnewTransformed <- model$transform_data(Dnew)


## -----------------------------------------------------------------------------
# The value of the coefficients for horizon 1, before the update
model$Lfits$k1$theta
# Update the coefficients
rls_update(model, DnewTransformed, Dnew[[model$output]])
# The value of the coefficients for horizon 1, after the update
model$Lfits$k1$theta


## -----------------------------------------------------------------------------
# Calculate the predictions
yhat <- rls_predict(model, DnewTransformed)


## -----------------------------------------------------------------------------
# The index for the predicted steps ahead
iseq <- i+model$kseq
# Plot the observations and predictions
plot(D$t[iseq], D$heatload[iseq], type = "b", xlab = "t", ylab = "y")
lines(D$t[iseq], yhat, type = "b", col = 2)
legend("topright", c("Observations",pst("Predictions (",min(model$kseq)," to ",max(model$kseq)," steps ahead)")), lty = 1, col = 1:2)


## -----------------------------------------------------------------------------
# Fit and predict on entire data
val <- rls_fit(model$prm, model, D)
# Keep the forecasts
D$Yhat1 <- val$Yhat


## -----------------------------------------------------------------------------
# The indexes of training period
itrain <- which(in_range("2010-12-15",D$t,"2011-01-01"))
# The indexes of the test period
itest <- which(in_range("2011-01-01",D$t,"2011-01-04"))

# Fit on the training period
rls_fit(model$prm, model, subset(D, itrain))

# Prepare for the forecasts
D$Yhat2 <- data.frame(matrix(NA, nrow(D$Yhat1), ncol(D$Yhat1)))
names(D$Yhat2) <- names(D$Yhat1)

# Step through the test period:
for(i in itest){
    Dnew <- subset(D, i)
    Dnewtr <- model$transform_data(Dnew)
    rls_update(model, Dnewtr, Dnew[[model$output]])
    D$Yhat2[i, ] <- as.numeric(rls_predict(model, Dnewtr))
}


## -----------------------------------------------------------------------------
D$Yhat1$k1[itest] - D$Yhat2$k1[itest]