-
Notifications
You must be signed in to change notification settings - Fork 25
/
mbsts.R
330 lines (310 loc) · 14.5 KB
/
mbsts.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
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
# Copyright 2019 Steven L. Scott. All rights reserved.
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
# License as published by the Free Software Foundation; either
# version 2.1 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public
# License along with this library; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
# ===========================================================================
# This is a multivariate version of the univariate bsts function. For now it
# only supports Gaussian observations. In the future the plan is to expand to
# the same set of reward distributions supported by bsts. This function may
# actually merge with bsts in the future. It is quite similar.
#
# Expected usage 1:
# Matrix y, Matrix x
# ss <- AddSharedLocalLevel(list(), y)
# model <- mbsts(y ~ x, state.specification = ss, niter = 1000)
#
# Expected usage 2:
# Data frame with columns for time stamp and series identifier.
#
# There are two data formats to support. Wide and tall. If you want regressors
# then you'll need tall. Tall data requires an "id" and a "timestamp" variable.
# The tall format is more general, but might be less convenient sometimes.
#
# With the tall format it becomes clear that there are two kinds of regressors.
# Those specific to a time point, and those specific to a series at a time
# point. The second form can be obtained from the first by repeating the
# predictors.
#
# There should be a role for shrinkage across series in choosing the prior.
# That comes later.
#
# Args:
# formula: Either a formula as one would supply to lm(), or a matrix of data.
# If a matrix, rows represent time, and columns represent time series. If a
# matrix is supplied then set data.format to "wide".
# shared.state.specification: A list of SharedStateSpecification objects
# defining the components of state that are shared across multiple series.
# series.state.specification: A list of StateSpecification objects defining
# the series-specific components of state. There are two options for how
# this argument can be formatted.
# - It can be a single list with StateSpecification objects as elements. In
# this case each object will be used to define the series-specific state
# for each time series in the response.
# - It can be a list with length matching the number of time series to be
# modeled. Each list element is a list of StateSpecification objects
# defining the series-specific state for the corresponding time series.
# This option is less convenient, but allows greater control.
# data: An optional data frame containing the data referenced in 'formula.'
# timestamps: A vector of timestamps indicating the time of each observation.
# If "wide" data are used this argument is optional, as the timestamps will
# be inferred from the rows of the data matrix.
# series.id: A factor-like object indicating the time series to which each
# observation belongs. This is only necessary for data in "long" format.
# prior: The prior distribution over any regression coefficients or other
# parameters of the observatino distribution.
# contrasts: an optional list containing the names of contrast functions to
# use when converting factors numeric variables in a regression formula.
# This argument works exactly as it does in 'lm'. The names of the list
# elements correspond to factor variables in your model formula. The list
# elements themselves are the names of contrast functions (see
# help(contrast) and the 'contrasts.arg' argument to
# 'model.matrix.default'). This argument is only used if a model formula is
# specified. It can usually be ignored even then.
# na.action: What to do about missing values. The default is to allow
# missing responses, but no missing predictors. Set this to na.omit or
# na.exclude if you want to omit missing responses altogether, but do so
# with care, as that will remove observations from the time series.
# niter: a positive integer giving the desired number of MCMC draws
# ping: A scalar. If ping > 0 then the program will print a status message to
# the screen every 'ping' MCMC iterations.
# data.format: Indicates whether the data are in wide or long form. Wide data
# will be converted to long format using Bsts::WideToLong.
# seed: The seed for the C++ random number generator.
# ...: Extra arguments are passed to DefaultMbstsPrior().
mbsts <- function(formula,
shared.state.specification,
series.state.specification = NULL,
data = NULL,
timestamps = NULL,
series.id = NULL,
prior = NULL, # TODO
opts = NULL,
contrasts = NULL,
na.action = na.pass,
niter,
ping = niter / 10,
data.format = c("long", "wide"),
seed = NULL,
...) {
data.format <- match.arg(data.format)
check.nonnegative.scalar(niter)
check.scalar.integer(ping)
stopifnot(is.null(seed) || length(seed) == 1)
if (!is.null(seed)) {
seed <- as.integer(seed)
}
# The first step is to properly format the data. We need an object called
# data.list that contains the following:
# - predictors: a matrix of predictors
# - response: a vector of responses
# - series: a factor indicating the time series each element of 'response'
# belongs to.
# - timestamp.info: A list created by TimestampInfo
#
# The coming if/then block is longer than you'd like, but it needs to take
# place here so the 'model.matrix' black magic can work.
has.regression <- is.language(formula)
if (has.regression) {
function.call <- match.call()
my.model.frame <- match.call(expand.dots = FALSE)
frame.match <- match(c("formula", "data", "na.action"),
names(my.model.frame), 0L)
my.model.frame <- my.model.frame[c(1L, frame.match)]
my.model.frame$drop.unused.levels <- TRUE
# In an ordinary regression model the default action for NA's is to delete
# them. This makes sense in ordinary regression models, but is dangerous in
# time series, because it artificially shortens the time between two data
# points. If the user has not specified an na.action function argument then
# we should use na.pass as a default, so that NA's are passed through to the
# underlying C++ code.
if (! "na.action" %in% names(my.model.frame)) {
my.model.frame$na.action <- na.pass
}
my.model.frame[[1L]] <- as.name("model.frame")
my.model.frame <- eval(my.model.frame, parent.frame())
model.terms <- attr(my.model.frame, "terms")
predictors <- model.matrix(model.terms, my.model.frame, contrasts)
if (any(is.na(predictors))) {
stop("bsts does not allow NA's in the predictors, only the responses.")
}
response <- model.response(my.model.frame, "any")
if (data.format == "wide") {
response.frame <- WideToLong(response)
if (is.null(timestamps)) {
timestamps <- response.frame$time
} else {
timestamps <- timestamps[response.frame$time]
}
series.id <- response.frame$series.id
expanded.predictors <- predictors[series.id, ]
response <- response.frame$values
}
stopifnot(length(response) == nrow(predictors))
} else {
## If there is no regression component the data could be in either long or
## wide format. Make sure it is long.
response <- formula
if (data.format == "wide") {
# Handle a data frames like a matrix.
if (is.data.frame(response) && all(sapply(response, is.numeric))) {
response <- as.matrix(response)
## TODO: Ensure we don't lose zoo timestamps here.
}
response.frame <- WideToLong(response)
timestamps <- response.frame$time
series.id <- response.frame$series
response <- response.frame$values
}
if (is.data.frame(response) || is.matrix(response)) {
stop("Please change data.format to 'wide' if passing matrix-valued time series")
}
predictors <- matrix(1, nrow = length(response), ncol = 1)
colnames(predictors) <- "Intercept"
}
series.id <- as.factor(series.id)
nseries <- length(levels(series.id))
data.list <- list(
predictors = predictors,
response = response,
series.id = series.id,
timestamp.info = TimestampInfo(predictors, NULL, timestamps)
)
#------------------------------------------------------------------------
# Check the format of the state specification.
#------------------------------------------------------------------------
spec <- .CheckMbstsStateSpecification(shared.state.specification,
series.state.specification, nseries)
shared.state.specification <- spec[[1]]
series.state.specification <- spec[[2]]
#------------------------------------------------------------------------
# Ensure the prior has the proper format.
#------------------------------------------------------------------------
if (is.null(prior)) {
prior <- .DefaultMbstsPrior(predictors, response, series.id)
}
# The prior for the observation model is a list of spike and slab priors.
stopifnot(is.list(prior), length(prior) == nseries,
all(sapply(prior, inherits, "SpikeSlabPrior")))
#------------------------------------------------------------------------
# Check that options is either NULL or a list. A bit of faith is needed that
# the list is formatted correctly.
# ------------------------------------------------------------------------
if (!is.null(opts)) {
stopifnot(is.list(opts))
}
#------------------------------------------------------------------------
# Check that all the scalars are actually scalars.
#------------------------------------------------------------------------
if (!is.null(seed)) {
seed <- as.integer(seed)
stopifnot(length(seed) == 1)
}
check.scalar.integer(niter)
check.scalar.integer(ping)
#------------------------------------------------------------------------
# Do the work!
#------------------------------------------------------------------------
ans <- .Call("analysis_common_r_fit_multivariate_bsts_model_",
data.list,
shared.state.specification,
series.state.specification,
prior,
opts,
as.integer(niter),
as.integer(ping),
seed)
ans$has.regression <- has.regression
ans$shared.state.specification <- shared.state.specification
ans$series.state.specification <- series.state.specification
ans$prior <- prior
ans$niter <- niter
ans$timestamp.info <- data.list$timestamp.info
ans$series.id <- series.id
ans$original.series <- response
ans$predictors <- predictors
#------------------------------------------------------------------------
# Final formatting.
#------------------------------------------------------------------------
# Set dimnames.
series.names <- as.character(levels(series.id))
predictor.names <- colnames(predictors)
state.model.names <- sapply(shared.state.specification,
function(x) x$name)
ans$nseries <- length(series.names)
dimnames(ans$regression.coefficients) <- list(
NULL, series.names, predictor.names)
dimnames(ans$shared.state.contributions) <- list(
NULL, state.model.names, series.names, NULL)
class(ans) <- "mbsts"
return(ans)
}
.DefaultMbstsPrior <- function(predictors, response, series) {
## Set a default prior on each of the regression models in the mbsts
## observation equation.
ans <- list()
for (s in sort(unique(series))) {
index <- series == s
data.list <- list(predictors = predictors[index, , drop = FALSE],
response = response[index])
ans[[s]] <- .SetDefaultPrior(data.list, family = "gaussian")
}
return(ans)
}
.CheckMbstsStateSpecification <- function(shared.state.specification,
series.state.specification,
nseries) {
# Check that the shared- and series-specific state specifictions are filled
# with legal values.
#
# Args:
# shared.state.specification: A list of SharedStateSpecification objects
# defining the components of state that are shared across multiple series.
# series.state.specification: A list of StateSpecification objects defining
# the series-specific components of state. There are two options for how
# this argument can be formatted.
# - It can be a single list with StateSpecification objects as elements. In
# this case each object will be used to define the series-specific state
# for each time series in the response.
# - It can be a list with length matching the number of time series to be
# modeled. Each list element is a list of StateSpecification objects
# defining the series-specific state for the corresponding time series.
# This option is less convenient, but allows greater control.
# nseries: The number of time series.
#
# Returns:
# A 2-item list containing (1) the shared.state.specification, and (2) the
# series.state.specification, possibly after expanding it.
stopifnot(is.list(shared.state.specification),
all(sapply(shared.state.specification, inherits, "SharedStateModel")))
# The series.state.specification might be an empty list.
stopifnot(is.null(series.state.specification)
|| is.list(series.state.specification))
# If series state specification is not NULL it is either a list of state
# specificiations to be repeated for each time series, or it is a list of such
# specifications.
if (!is.null(series.state.specification) &&
all(sapply(series.state.specification, inherits, "StateModel"))) {
series.state.specification <- RepList(series.state.specification, nseries)
}
if (!is.null(series.state.specification)) {
stopifnot(is.list(series.state.specification),
length(series.state.specification) == nseries)
for (i in 1:nseries) {
stopifnot(is.null(series.state.specification[[i]]) ||
is.list(series.state.specification[[i]]),
all(sapply(series.state.specification[[i]], inherits, "StateModel")))
}
}
return(list(shared.state.specification, series.state.specification))
}