-
Notifications
You must be signed in to change notification settings - Fork 25
/
add.regression.holiday.R
301 lines (285 loc) · 12.5 KB
/
add.regression.holiday.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
.ValidateHolidayList <- function(holiday.list) {
## Args:
## holiday.list: a list of objects of class Holiday.
## Returns:
## If holiday.list is a list of objects of class Holiday then it is returned
## unchanged. Otherwise an error is raised.
stopifnot(is.list(holiday.list),
all(sapply(holiday.list, inherits, "Holiday")))
}
.DefaultRegressionHolidayModelCoefficientPrior <- function(sdy) {
## Args:
## sdy: The standard deviation of the target series.
## Returns:
## A default prior of class NormalPrior, describing variation (across days
## and across holidays) among daily holiday effects.
return(NormalPrior(0, sdy))
}
AddRegressionHoliday <- function(state.specification = NULL,
y,
holiday.list,
time0 = NULL,
prior = NULL,
sdy = sd(as.numeric(y), na.rm = TRUE)) {
## Add a regression-based state model describing the effects of one or more
## holidays. Each day of each holiday exerts a constant effect, specific to
## that day, on the observed time series.
##
## Args:
## state.specification: A list of state components. If omitted, an empty
## list is assumed.
## y: A numeric vector. The time series being modeled. See 'sdy' below.
## holiday.list: A list of objects of type 'Holiday'. See ?Holiday. The
## width of the influence window should be the same number of days for all
## the holidays in this list. See below if this case does not apply.
## time0: Either NULL, or an object convertible to class Date giving the
## date of the first day in the training data y. If NULL and y is of type
## 'zoo' then an attempt will be made to infer time0 from the index of y.
## prior: An object of type NormalPrior describing the a priori
## expected variation among holiday effects.
## sdy: The standard deviation of the 'y' argument, used to speicfy default
## priors in case prior is missing.
##
## Returns:
## A list containing the information needed to specify this state model to
## the underlying C++ code, in the expected format.
.ValidateHolidayList(holiday.list)
if (missing(state.specification)) {
state.specification <- list()
}
if (is.null(prior)) {
prior <- .DefaultRegressionHolidayModelCoefficientPrior(sdy)
}
stopifnot(inherits(prior, "NormalPrior"))
spec <- list(name = "RegressionHolidays",
holidays = holiday.list,
time0 = as.Date(.SetTimeZero(time0, y)),
size = 1,
prior = prior)
class(spec) <- c("RegressionHolidayStateModel", "HolidayStateModel", "StateModel")
state.specification[[length(state.specification) + 1]] <- spec
return(state.specification)
}
.PlotHolidayRegressionCoefficients <- function(coefficients, ylim, ...) {
## A driver function used to implement common plotting needs for
## plot.RegressionHolidayStateModel and
## plot.HierarchicalRegressionHolidayStateModel.
##
## Args:
## coefficients: A matrix of holiday coefficients. Each row is a Monte Carlo
## draw, and each column is a specific holiday coefficient. The column
## names for the matrix are expected to give the name of the holiday
## corresponding to each coefficient. All coefficients for a given
## holiday window are expected to be grouped together in ascending order.
## ylim: Limits for the vertical axis.
## ...: Extra arguments passed to boxplot.
##
## Side Effects:
## A plot is added to the current graphics device. Side-by-side boxplots
## show the posterior distribution of the effect of each holiday.
##
## Returns:
## invisible(NULL)
if(is.null(ylim)) {
ylim = range(coefficients)
}
variable.names <- colnames(coefficients)
old.mai <- par("mai")
old.mai[1] <- max(strheight(variable.names, units = "inches"))
old.mai[2] <- max(strwidth(variable.names, units = "inches")) +
strwidth(" -- ", units = "inches")
oldpar <- par(mai = old.mai)
on.exit(par(oldpar))
## Reversing the order of the coefficients means the plot can be read from top
## to bottom, which is the natural order.
new.holiday.index <- (1:ncol(coefficients))[colnames(coefficients) != ""][-1]
coefficients <- coefficients[, ncol(coefficients):1]
new.holiday.index <- 1 + ncol(coefficients) - new.holiday.index
boxplot(coefficients, horizontal = TRUE, ylim = ylim, las = 1, cex = .6, ...)
abline(v = 0, lty = 3)
abline(h = new.holiday.index + .5, lty = 3)
return(invisible(NULL))
}
plot.RegressionHolidayStateModel <- function(x,
bsts.object,
burn = NULL,
time = NULL,
style = NULL,
ylim = NULL,
...) {
## S3 method for plotting a RegressionHolidayStateModel
##
## Args:
## x: An object inheriting from RegressionHolidayStateModel.
## bsts.object: A bsts model that includes state.specification in its state
## specification.
## burn: The number of MCMC iterations to discard as burn-in.
## time: Not used. Here to match the signature of plot.StateModel.
## style: Not used. Here to match the signature of plot.StateModel.
## ylim: Limits on the vertical axis.
## ...: Extra arguments passed boxplot
##
## Side Effects:
## A plot is added to the current graphics device. Side-by-side boxplots
## show the posterior distribution of the effect of each holiday.
##
## Returns:
## invisible(NULL)
state.specification <- x
stopifnot(inherits(state.specification, "RegressionHolidayStateModel"))
stopifnot(inherits(bsts.object, "bsts"))
if (is.null(.FindStateSpecification(state.specification, bsts.object))) {
stop("The state specification is not part of the bsts object.")
}
holidays <- state.specification$holidays
coefficients <- NULL
for (i in 1:length(holidays)) {
new.coefficients <- bsts.object[[holidays[[i]]$name]]
colnames(new.coefficients) <-
c(holidays[[i]]$name, rep("", ncol(new.coefficients) - 1))
coefficients <- cbind(coefficients, new.coefficients)
}
if (is.null(burn)) {
burn <- 0
}
stopifnot(is.numeric(burn))
if (burn > 0) {
coefficients <- coefficients[-(1:burn), , drop = FALSE]
}
.PlotHolidayRegressionCoefficients(coefficients, ylim, ...)
title("Regression Holiday Effects")
return(invisible(NULL))
}
plot.HierarchicalRegressionHolidayStateModel <- function(x,
bsts.object,
burn = NULL,
time = NULL,
style = NULL,
ylim = NULL,
...) {
## S3 method for plotting a HierarchicalRegressionHolidayStateModel.
##
## Args:
## x: An object inheriting from HierarchicalRegressionHolidayStateModel.
## bsts.object: A bsts model that includes state.specification in its state
## specification.
## burn: The number of MCMC iterations to discard as burn-in.
## time: Not used. Here to match the signature of plot.StateModel.
## style: Not used. Here to match the signature of plot.StateModel.
## ylim: Limits on the vertical axis.
## ...: Extra arguments passed boxplot
##
## Side Effects:
## A plot is added to the current graphics device. Side-by-side boxplots
## show the posterior distribution of the effect of each holiday.
##
## Returns:
## invisible(NULL)
state.specification <- x
stopifnot(inherits(
state.specification, "HierarchicalRegressionHolidayStateModel"))
stopifnot(inherits(bsts.object, "bsts"))
if (is.null(.FindStateSpecification(state.specification, bsts.object))) {
stop("The state specification is not part of the bsts object.")
}
coefficient.array <- bsts.object$holiday.coefficients
if (is.null(burn)) {
burn <- 0
}
if (burn > 0) {
coefficient.array <- coefficient.array[-(1:burn), , , drop = FALSE]
}
holiday.window.size <- dim(coefficient.array)[3]
number.of.holidays <- dim(coefficient.array)[2]
coefficients <- matrix(nrow = dim(coefficient.array)[1],
ncol = holiday.window.size * number.of.holidays)
colnames(coefficients) <- rep("", ncol(coefficients))
colnames(coefficients)[seq(from = 1, by = holiday.window.size,
len = number.of.holidays)] <- dimnames(coefficient.array)[[2]]
start <- 0
for (i in 1:number.of.holidays) {
coefficients[, start + (1:holiday.window.size)] <- coefficient.array[, i, ]
start <- start + holiday.window.size
}
.PlotHolidayRegressionCoefficients(coefficients, ylim, ...)
title("Hierarchical Regression\nHoliday Effects")
return(invisible(NULL))
}
AddHierarchicalRegressionHoliday <- function(
state.specification = NULL,
y,
holiday.list,
coefficient.mean.prior = NULL,
coefficient.variance.prior = NULL,
time0 = NULL,
sdy = sd(as.numeric(y), na.rm = TRUE)) {
## Args:
## state.specification: A list of state components. If omitted, an empty
## list is assumed.
## y: A numeric vector. The time series being modeled. See 'sdy' below.
## holiday.list: A list of objects of type 'Holiday'. See ?Holiday. The
## width of the influence window should be the same number of days for all
## the holidays in this list. See below if this case does not apply.
## coefficient.mean.prior: An object of type MvnPrior giving the hyperprior
## for the average effect of a holiday in each day of the influence
## window.
## coefficient.variance.prior: An object of type InverseWishartPrior
## describing the prior belief about the variation in holiday effects from
## one holiday to the next.
## time0: Either NULL, or an object convertible to class Date giving the
## date of the first day in the training data y. If NULL and y is of type
## 'zoo' then an attempt will be made to infer time0 from the index of y.
## sdy: The standard deviation of the 'y' argument, used to speicfy default
## priors in case either coefficient.mean.prior or
## coefficient.variance.prior are missing. If both coefficient.mean.prior
## and coefficient.variance.prior are specified then neither 'y' nor 'sdy'
## is needed. If 'sdy' is specified then 'y' is not needed.
##
## Returns:
## A list containing the information needed to specify this state model to
## the C++ code, in the expected format.
## Make sure the holiday list is a list of Holiday objects, and that it
## contains enough holidays to be useful.
.ValidateHolidayList(holiday.list)
if (length(holiday.list) < 3) {
stop("You need 3 or more holidays to fit the hierarchical model in ",
"AddHierarchicalRegressionHolidayModel.")
}
if (missing(state.specification)) {
state.specification <- list()
}
if (is.null(coefficient.mean.prior) || is.null(coefficient.variance.prior)) {
if (missing(sdy)) {
stopifnot(is.numeric(y))
sdy <- sd(as.numeric(y), na.rm = TRUE)
}
max.window.width <- unique(sapply(holiday.list, MaxWindowWidth))
if (length(max.window.width) != 1) {
stop("All holidays must have the same window width.")
}
}
if (is.null(coefficient.mean.prior)) {
coefficient.mean.prior <- MvnPrior(
mean = rep(0, max.window.width),
variance = diag(rep(sdy, max.window.width))
)
}
stopifnot(inherits(coefficient.mean.prior, "MvnPrior"))
if (is.null(coefficient.variance.prior)) {
coefficient.variance.prior <- InverseWishartPrior(
variance.guess = diag(rep(sdy / 10, max.window.width)),
variance.guess.weight = max.window.width + 1)
}
stopifnot(inherits(coefficient.variance.prior,
"InverseWishartPrior"))
spec <- list(name = "HierarchicalRegressionHolidays",
holidays = holiday.list,
time0 = as.Date(.SetTimeZero(time0, y)),
size = 1,
coefficient.mean.prior = coefficient.mean.prior,
coefficient.variance.prior = coefficient.variance.prior)
class(spec) <- c("HierarchicalRegressionHolidayStateModel",
"HolidayStateModel", "StateModel")
state.specification[[length(state.specification) + 1]] <- spec
return(state.specification)
}