Intro

This vignette presents an example of using the onlineforecasting package for fitting a model and calculating forecasts as carried out by Bacher et al. (2013).

Data

Data for the forecasting examples are taken from a data set collected in Sønderborg, Denmark. It comprises heat load measurements for sixteen houses, together with local climate observations and weather forecasts (NWPs). The houses are generally built in the sixties and seventies, with a floor plan in the range of 85 to 170 \(\mathrm{m^2}\), and constructed in bricks. For each house only the total heat load, including both space heating and hot tap water heating, is available. The climate observations are measured at the local district heating plant within 10 kilometers from the houses. The NWPs are from the HIRLAM-S05 model and provided by the Danish Meteorological Institute. All times are in UTC and the time stamp for average values are set to the end of the time interval.

Load the package:

library(onlineforecast)

The Dbuilding data is included in the package. Its a data.list (see the vignette onlineforecasting.pdf):

class(Dbuilding)
##     [1] "data.list" "list"

Keep it in D and see the content:

D <- Dbuilding
names(D)
##     [1] "t"             "heatload"      "heatloadtotal" "Taobs"        
##     [5] "Iobs"          "Ta"            "I"

The time:

head(D$t)
##     [1] "2010-12-15 01:00:00 UTC" "2010-12-15 02:00:00 UTC"
##     [3] "2010-12-15 03:00:00 UTC" "2010-12-15 04:00:00 UTC"
##     [5] "2010-12-15 05:00:00 UTC" "2010-12-15 06:00:00 UTC"

The observed heat load (in kW) of a house is kept in a data.frame:

head(D$heatload)
##     [1] 5.92 5.85 5.85 5.88 5.85 5.83

The Numerical Weather Predictions (NWPs) of ambient temperature steps 1 to 9 hours ahead are:

head(D$Ta[ ,1:9])
##          k1    k2    k3    k4    k5    k6    k7    k8    k9
##     1 -2.82 -3.20 -3.12 -3.09 -3.13 -3.16 -3.17 -3.09 -2.77
##     2 -2.90 -3.12 -3.09 -3.13 -3.16 -3.17 -3.09 -2.77 -2.32
##     3 -2.94 -3.09 -3.13 -3.16 -3.17 -3.09 -2.77 -2.32 -1.95
##     4 -2.89 -3.11 -3.05 -3.11 -3.12 -2.81 -2.37 -2.01 -1.81
##     5 -2.81 -3.05 -3.11 -3.12 -2.81 -2.37 -2.01 -1.81 -1.86
##     6 -2.77 -3.11 -3.12 -2.81 -2.37 -2.01 -1.81 -1.86 -2.28

So at “2010-12-15 01:00:00 GMT” the latest available forecasts is the first row of Ta.

We will add a y value into the data.list for simplicity and convention:

D$y <- D$heatload

Create a time of the day matrix for the forecast

D$tday <- make_tday(D$t, 1:36)
head(D$tday)
##       k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 k13 k14 k15 k16 k17 k18 k19 k20 k21
##     1  2  3  4  5  6  7  8  9 10  11  12  13  14  15  16  17  18  19  20  21  22
##     2  3  4  5  6  7  8  9 10 11  12  13  14  15  16  17  18  19  20  21  22  23
##     3  4  5  6  7  8  9 10 11 12  13  14  15  16  17  18  19  20  21  22  23   0
##     4  5  6  7  8  9 10 11 12 13  14  15  16  17  18  19  20  21  22  23   0   1
##     5  6  7  8  9 10 11 12 13 14  15  16  17  18  19  20  21  22  23   0   1   2
##     6  7  8  9 10 11 12 13 14 15  16  17  18  19  20  21  22  23   0   1   2   3
##       k22 k23 k24 k25 k26 k27 k28 k29 k30 k31 k32 k33 k34 k35 k36
##     1  23   0   1   2   3   4   5   6   7   8   9  10  11  12  13
##     2   0   1   2   3   4   5   6   7   8   9  10  11  12  13  14
##     3   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15
##     4   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16
##     5   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17
##     6   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18

A time series plot, see “?plot_ts.data.list” (Note how the forecasts Ta are lagged to be synced with observations (i.e. then also with each other)):

plot_ts(D, c("^y","Ta"), kseq=c(1,12))