## -----------------------------------------------------------------------------------------------------------------
library(onlineforecast)


## -----------------------------------------------------------------------------------------------------------------
class(Dbuilding)


## -----------------------------------------------------------------------------------------------------------------
D <- Dbuilding
D$tday <- make_tday(D$t, 1:36)
names(D)


## -----------------------------------------------------------------------------------------------------------------
D$y <- D$heatload


## -----------------------------------------------------------------------------------------------------------------
Dtrain <- subset(D, c("2011-01-08","2011-01-25"))
Dtrain$scoreperiod <- in_range("2011-01-15", Dtrain$t)


## -----------------------------------------------------------------------------------------------------------------
model <- forecastmodel$new()
model$output = "y"
model$add_inputs(Ta = "Ta",     # The ambient temperature
                 mu = "one()") # one() generates a matrix of ones (i.e. an intercept is included)


## -----------------------------------------------------------------------------------------------------------------
model$kseq <- 1:6
datatr <- model$transform_data(Dtrain)


## -----------------------------------------------------------------------------------------------------------------
names(datatr)


## -----------------------------------------------------------------------------------------------------------------
class(datatr$Ta)
names(datatr$Ta)
head(datatr$Ta)
##
class(datatr$mu)
names(datatr$mu)
head(datatr$mu)


## -----------------------------------------------------------------------------------------------------------------
k <- 6


## -----------------------------------------------------------------------------------------------------------------
X <- as.data.frame(subset(datatr, kseq=k, lagforecasts = TRUE))
names(X)


## -----------------------------------------------------------------------------------------------------------------
X <- cbind(y = Dtrain[[model$output]], X)


## -----------------------------------------------------------------------------------------------------------------
frml <- pst("y ~ ", pst(names(X)[-1], collapse=" + "), " - 1")
frml


## -----------------------------------------------------------------------------------------------------------------
fit <- lm(frml, X)


## -----------------------------------------------------------------------------------------------------------------
yhat <- predict.lm(fit, X)


## -----------------------------------------------------------------------------------------------------------------
plot(X$y)
lines(yhat)


## -----------------------------------------------------------------------------------------------------------------
rmse(X$y - yhat)


## -----------------------------------------------------------------------------------------------------------------
model <- forecastmodel$new()
model$output = "y"
model$add_inputs(Ta = "lp(Ta, a1=0.85)", # The ambient temperature through a low-pass filter
                 mu = "one()")


## -----------------------------------------------------------------------------------------------------------------
class(model$inputs)
model$inputs$Ta$expr


## -----------------------------------------------------------------------------------------------------------------
tmp <- model$inputs$Ta$evaluate(D)
head(tmp[ ,1:5])


## -----------------------------------------------------------------------------------------------------------------
model$kseq <- 1:36
datatr <- model$transform_data(Dtrain)


## -----------------------------------------------------------------------------------------------------------------
plot(Dtrain$Ta$k6)
lines(datatr$Ta$k6)


## -----------------------------------------------------------------------------------------------------------------
X <- as.data.frame(subset(datatr, kseq=k, lagforecasts = TRUE))
X <- cbind(y = Dtrain[[model$output]], X)
fit <- lm(frml, X)


## -----------------------------------------------------------------------------------------------------------------
yhat <- predict.lm(fit, X)


## -----------------------------------------------------------------------------------------------------------------
plot(X$y)
lines(yhat)


## -----------------------------------------------------------------------------------------------------------------
rmse(X$y - yhat)


## -----------------------------------------------------------------------------------------------------------------
model$insert_prm(c(Ta__a1 = 0.99))
model$inputs$Ta$expr


## -----------------------------------------------------------------------------------------------------------------
datatr <- model$transform_data(Dtrain)
plot(Dtrain$Ta$k1)
lines(datatr$Ta$k1)


## -----------------------------------------------------------------------------------------------------------------
scorefun <- function(prm, model, k){
  model$reset_state()
  model$insert_prm(prm)
  datatr <- model$transform_data(Dtrain)
  X <- as.data.frame(subset(datatr, kseq=k, lagforecasts = TRUE))
  X <- cbind(y = Dtrain[[model$output]], X)
  frml <- pst("y ~ ", pst(names(X)[-1], collapse=" + "), " - 1")
  yhat <- predict(lm(frml, X), newdata = X)
  return(rmse(X$y - yhat))
}


## -----------------------------------------------------------------------------------------------------------------
scorefun(c(Ta__a1 = 0.85), model, k)


## -----------------------------------------------------------------------------------------------------------------
optim(c(Ta__a1 = 0.85), scorefun, model = model, k = k)


## -----------------------------------------------------------------------------------------------------------------
prm <- c(Ta__a1 = 0.938)
model$reset_state()
model$insert_prm(prm)
datatr <- model$transform_data(Dtrain)
X <- as.data.frame(subset(datatr, kseq=k, lagforecasts = TRUE))
X <- cbind(y = Dtrain[[model$output]], X)
fit <- lm(frml, X)
yhat <- predict.lm(fit, X)
plot(X$y)
lines(yhat)


## -----------------------------------------------------------------------------------------------------------------
kseq <- 1:36
rmsek <- sapply(kseq, function(k){
  scorefun(prm, model, k = k)
})
plot(1:36, rmsek)


## -----------------------------------------------------------------------------------------------------------------
model$kseq <- 6
lm_fit(prm, model, Dtrain, rmse)


## -----------------------------------------------------------------------------------------------------------------
model$add_prmbounds(Ta__a1 = c(0.8, 0.9, 0.999))
model$prmbounds


## -----------------------------------------------------------------------------------------------------------------
init <- model$get_prmbounds("init")
init


## -----------------------------------------------------------------------------------------------------------------
optim(par = init,
      fn = lm_fit,
      model = model,
      data = Dtrain,
      scorefun = rmse,
      printout = FALSE,
      returnanalysis = FALSE,
      lower = model$prmbounds[ ,"lower"],
      upper = model$prmbounds[ ,"upper"])


## -----------------------------------------------------------------------------------------------------------------
model$prm <- lm_optim(model, Dtrain, printout = FALSE)$par


## -----------------------------------------------------------------------------------------------------------------
model$add_inputs(I = "lp(I, a1=0.85)") # The global radiation through a low-pass filter
model$add_prmbounds(I__a1 = c(0.5, 0.8, 0.99))
lm_optim(model, Dtrain, printout = FALSE)


## -----------------------------------------------------------------------------------------------------------------
model$add_inputs(mu_tday = "fs(tday/24, 3)") # A diurnal curve using Fourier series with 3 harmonics (pairs of sine and cosines)
lm_optim(model, Dtrain, printout = FALSE)


## -----------------------------------------------------------------------------------------------------------------
model$prm <- lm_optim(model, Dtrain, printout = FALSE)$par


## -----------------------------------------------------------------------------------------------------------------
model$kseq <- 1:36
lm_fit(model$prm, model, Dtrain)
datatr <- model$transform_data(Dtrain)
Yhat <- lm_predict(model, datatr)
plot(X$y)
lines(lagvec(Yhat[ ,"k1"],1))
lines(lagvec(Yhat[ ,"k18"],18), col="red")


## -----------------------------------------------------------------------------------------------------------------
i <- 49:(96+36)
ifor <- 96
plot(Dtrain$t[i], Dtrain$y[i])
lines(Dtrain$t[ifor+model$kseq], Yhat[ifor, ], col="red")
abline(v = Dtrain$t[ifor])


## -----------------------------------------------------------------------------------------------------------------
val <- lm_fit(model$prm, model, Dtrain, returnanalysis = TRUE)
Residuals <- Dtrain$y - lagdf(val$Yhat, "+k")
plot(apply(Residuals, 2, rmse), type="b")


## -----------------------------------------------------------------------------------------------------------------
Tacoef <- sapply(getse(model$Lfits, "coefficients"), function(x){ x[grep("^Ta",names(x))] })
plot(Tacoef)