1
2#' State Space Model Predictions
3#'
4#' Function \code{predict.SSModel} predicts the future observations of a state
5#' space model of class \code{\link{SSModel}}.
6#'
7#' For non-Gaussian models, the results depend whether importance sampling is
8#' used (\code{nsim>0}). without simulations, the confidence intervals are based
9#' on the Gaussian approximation of \eqn{p(\alpha | y)}. Confidence intervals in
10#' response scale are computed in linear predictor scale, and then transformed
11#' to response scale. The prediction intervals are not supported. With
12#' importance sampling, the confidence intervals are computed as the empirical
13#' quantiles from the weighted sample, whereas the prediction intervals contain
14#' additional step of simulating the response variables from the sampling
15#' distribution \eqn{p(y|\theta^i)}.
16#'
17#' Predictions take account the uncertainty in state estimation
18#' (given the prior distribution for the initial states), but not the uncertainty
19#' of estimating the parameters in the system matrices (i.e. \eqn{Z}, \eqn{Q} etc.).
20#' Thus the obtained confidence/prediction intervals can underestimate the true
21#' uncertainty for short time series and/or complex models.
22#'
23#' If no simulations are used, the standard errors in response scale are
24#' computed using the Delta method.
25#'
26#' @export
27#' @importFrom stats predict end qnorm
28#' @aliases predict predict.SSModel
29#' @param object Object of class \code{SSModel}.
30#' @param newdata A compatible \code{SSModel} object to be added in the end of
31#'   the old object for which the predictions are required. If omitted,
32#'   predictions are either for the past data points, or if argument
33#'   \code{n.ahead} is given, \code{n.ahead} time steps ahead.
34#' @param n.ahead Number of steps ahead at which to predict. Only used if
35#'   \code{newdata} is omitted. Note that when using \code{n.ahead}, object
36#'   cannot contain time varying system matrices.
37#' @param interval Type of interval calculation.
38#' @param level Confidence level for intervals.
39#' @param type Scale of the prediction, \code{"response"} or \code{"link"}.
40#' @param states Which states are used in computing the predictions. Either a
41#'   numeric vector containing the indices of the corresponding states, or a
42#'   character vector defining the types of the corresponding states. Possible choices are
43#'   \code{"all"},  \code{"level"}, \code{"slope"} (which does not make sense as the corresponding Z is zero.),
44#'   \code{"trend"},  \code{"regression"}, \code{"arima"}, \code{"custom"},
45#'   \code{"cycle"} or \code{"seasonal"}, where \code{"trend"} extracts all states relating to trend.
46#'    These can be combined. Default is \code{"all"}.
47#' @param nsim Number of independent samples used in importance sampling. Used
48#'   only for non-Gaussian models.
49#' @param se.fit If TRUE, standard errors of fitted values are computed. Default is FALSE.
50#' @param prob if TRUE (default), the predictions in binomial case are
51#'   probabilities instead of counts.
52#' @param maxiter The maximum number of iterations used in approximation Default
53#'   is 50. Only used for non-Gaussian model.
54#' @param filtered If \code{TRUE}, compute predictions based on filtered
55#' (one-step-ahead) estimates. Default is FALSE i.e. predictions are based on
56#' all available observations given by user. For diffuse phase,
57#' interval bounds and standard errors of fitted values are set to \code{-Inf}/\code{Inf}
58#' (If the interest is in the first time points it might be useful to use
59#' non-exact diffuse initialization.).
60#' @param expected Logical value defining the approximation of H_t in case of Gamma
61#' and negative binomial distribution. Default is \code{FALSE} which matches the
62#' algorithm of Durbin & Koopman (1997), whereas \code{TRUE} uses the expected value
63#' of observations in the equations, leading to results which match with \code{glm} (where applicable).
64#' The latter case was the default behaviour of KFAS before version 1.3.8.
65#' Essentially this is the difference between observed and expected information in GLM context.
66#' @param \dots Ignored.
67#' @return A matrix or list of matrices containing the predictions, and
68#'   optionally standard errors.
69#' @examples
70#'
71#' set.seed(1)
72#' x <- runif(n=100,min=1,max=3)
73#' y <- rpois(n=100,lambda=exp(x-1))
74#' model <- SSModel(y~x,distribution="poisson")
75#' xnew <- seq(0.5,3.5,by=0.1)
76#' newdata <- SSModel(rep(NA,length(xnew))~xnew,distribution="poisson")
77#' pred <- predict(model,newdata=newdata,interval="prediction",level=0.9,nsim=100)
78#' plot(x=x,y=y,pch=19,ylim=c(0,25),xlim=c(0.5,3.5))
79#' matlines(x=xnew,y=pred,col=c(2,2,2),lty=c(1,2,2),type="l")
80#'
81#' model <- SSModel(Nile~SSMtrend(1,Q=1469),H=15099)
82#' pred <- predict(model,n.ahead=10,interval="prediction",level=0.9)
83#' pred
84predict.SSModel <- function(object, newdata, n.ahead,
85  interval = c("none", "confidence", "prediction"), level = 0.95,
86  type = c("response", "link"), states = NULL, se.fit = FALSE,  nsim = 0,
87  prob = TRUE, maxiter = 50, filtered = FALSE, expected = FALSE, ...) {
88
89  interval <- match.arg(interval)
90  type <- match.arg(type)
91  # Check that the model object is of proper form
92  is.SSModel(object, na.check = TRUE, return.logical = FALSE)
93  m <- attr(object, "m")
94  p <- attr(object, "p")
95  k <- attr(object, "k")
96  if (missing(states)) {
97    states <- as.integer(1:m)
98  } else {
99    if (is.numeric(states)) {
100      states <- as.integer(states)
101      if (min(states) < 1 | max(states) > m)
102        stop("Vector states should contain the indices or types of the states which are combined.")
103    } else {
104      states <- match.arg(arg = states, choices = c("all", "arima", "custom", "level","slope",
105        "cycle", "seasonal", "trend", "regression"),
106        several.ok = TRUE)
107      if ("all" %in% states) {
108        states <- as.integer(1:attr(object, "m"))
109      } else {
110        if("trend" %in% states)
111          states <- c(states, "level", "slope")
112        states <- which(attr(object, "state_types") %in% states)
113      }
114    }
115  }
116  gaussianmodel <- all(object$distribution == "gaussian")
117  if (!missing(newdata) && !is.null(newdata)) {
118    #expand model with newdata object
119    # Check that the model object is of proper form
120    is.SSModel(newdata, na.check = TRUE, return.logical = FALSE)
121    if (p != attr(newdata, "p"))
122      stop("Different number of time series for 'object' and 'newdata'.")
123    if (m != attr(newdata, "m"))
124      stop("Different number of states for 'object' and 'newdata'.")
125    if (k != attr(newdata, "k"))
126      stop("Different number of disturbance terms for 'object' and 'newdata'.")
127    if (!identical(object$distribution, newdata$distribution))
128      stop("Different distributions for 'object' and 'newdata'")
129    no <- attr(object, "n")
130    nn <- attr(newdata, "n")
131    n <- attr(object, "n") <- no + nn
132    timespan <- (no + 1):n
133    object$y <- ts(rbind(object$y, newdata$y),
134      start = start(object$y), frequency = frequency(object$y))
135    endtime <- end(object$y)
136    tvo <- attr(object, "tv")
137    tvn <- attr(newdata, "tv")
138    same <- function(x, y) isTRUE(all.equal(x, y, tolerance = 0,
139      check.attributes = FALSE))
140    if (tvo[1] || tvn[1] || !same(object$Z, newdata$Z)) {
141      object$Z <- array(data = c(array(object$Z, dim = c(m, p, no)),
142        array(newdata$Z, dim = c(m, p, nn))), dim = c(p, m, n))
143      attr(object, "tv")[1] <- 1L
144    }
145    if (gaussianmodel && (tvo[2] || tvn[2] || !same(object$H, newdata$H))) {
146      object$H <- array(data = c(array(object$H, dim = c(p, p, no)),
147        array(newdata$H, dim = c(p, p, nn))), dim = c(p, p, n))
148      attr(object, "tv")[2] <- 1L
149    } else if(!gaussianmodel) object$u <- rbind(object$u, matrix(newdata$u, nn, p))
150
151    if (tvo[3] || tvn[3] || !same(object$T, newdata$T)) {
152      object$T <- array(data = c(array(object$T, dim = c(m, m, no)),
153        array(newdata$T, dim = c(m, m, nn))), dim = c(m, m, n))
154      attr(object, "tv")[3] <- 1L
155    }
156    if (tvo[4] || tvn[4] || !same(object$R, newdata$R)) {
157      object$R <- array(data = c(array(object$R, dim = c(m, k, no)),
158        array(newdata$R, dim = c(m, k, nn))), dim = c(m, k, n))
159      attr(object, "tv")[4] <- 1L
160    }
161    if (tvo[5] || tvn[5] || !same(object$Q, newdata$Q)) {
162      object$Q <- array(data = c(array(data = object$Q, dim = c(k, k, no)),
163        array(data = newdata$Q, dim = c(k, k, nn))), dim = c(k, k, n))
164      attr(object, "tv")[5] <- 1L
165    }
166  } else {
167    if (!missing(n.ahead) && !is.null(n.ahead)) {
168      tv <- attr(object, "tv")
169      if(ifelse(gaussianmodel,any(tv), any(c(apply(object$u, 2, function(x) length(unique(x)) > 1))) || any(tv[-2])))
170        stop("Model contains time varying system matrices, cannot use argument 'n.ahead'. Use 'newdata' instead.")
171      timespan <- attr(object, "n") + 1:n.ahead
172      n <- attr(object, "n") <- attr(object, "n") + as.integer(n.ahead)
173      endtime<-end(object$y) + c(0, n.ahead)
174      object$y <- window(object$y, end = endtime, extend = TRUE)
175      if (any(object$distribution != "gaussian"))
176        object$u <- rbind(object$u, matrix(object$u[1, ], nrow = n.ahead,
177          ncol = ncol(object$u), byrow = TRUE))
178    } else {
179      timespan <- 1:attr(object, "n")
180      endtime <- end(object$y)
181    }
182  }
183  if (!gaussianmodel && interval == "prediction") {
184    if (type == "link")
185      stop("Prediction intervals can only be computed at response scale.")
186    if (nsim < 1)
187      stop("Cannot compute prediction intervals for non-gaussian models without importance sampling. Use 'nsim' argument.")
188  }
189  pred <- vector("list", length = p)
190  if (gaussianmodel) { #Gaussian case
191    if (identical(states, as.integer(1:m))) {
192      if (filtered) {
193        out <- KFS(model = object, filtering = "mean", smoothing = "none")
194        names(out)[match(c("m", "P_mu"), names(out))]  <- c("muhat", "V_mu")
195        if (out$d > 0) {
196        out$V_mu[,,1:out$d] <- Inf #diffuse phase
197        }
198      } else {
199        out <- KFS(model = object, filtering = "none", smoothing = "mean")
200      }
201    } else {
202      if (filtered) {
203        out <- KFS(model = object, filtering = "state", smoothing = "none")
204        d <- out$d
205        out <- signal(out,  states = states, filtered = TRUE)
206        names(out) <- c("muhat", "V_mu")
207        if (d > 0) {
208        out$V_mu[,,1:d] <- Inf #diffuse phase
209        }
210      } else {
211        out <- signal(KFS(model = object, filtering = "none", smoothing = "state"),
212          states = states)
213        names(out) <- c("muhat", "V_mu")
214      }
215
216    }
217    for (i in 1:p) {
218      pred[[i]] <- cbind(fit = out$muhat[timespan, i],
219        switch(interval, none = NULL,
220          confidence = out$muhat[timespan, i] +
221            qnorm((1 - level)/2) * sqrt(out$V_mu[i, i, timespan]) %o% c(1, -1),
222          prediction = out$muhat[timespan, i] + qnorm((1 - level)/2) *
223            sqrt(out$V_mu[i, i, timespan] +
224                object$H[i, i, if (dim(object$H)[3] > 1) timespan else 1]) %o% c(1, -1)),
225        se.fit = if (se.fit)
226          sqrt(out$V_mu[i, i, timespan]))
227      if (interval != "none")
228        colnames(pred[[i]])[2:3] <- c("lwr", "upr")
229    }
230  } else {
231    if (nsim < 1) { #Using approximating model
232      if (identical(states, as.integer(1:m))) {
233        if (filtered) {
234          out <- KFS(model = object, filtering = "signal", smoothing = "none",
235            maxiter = maxiter, expected = expected)
236          names(out)[match(c("t", "P_theta"), names(out))] <- c("thetahat", "V_theta")
237          if (out$d > 0) {
238          out$V_theta[,,1:out$d] <- Inf #diffuse phase
239          }
240        } else {
241          out <- KFS(model = object, smoothing = "signal", maxiter = maxiter, expected = expected)
242        }
243      } else {
244        if (filtered) {
245          out <- KFS(model = object, filtering = "state", smoothing = "none",
246            maxiter = maxiter, expected = expected)
247          d <- out$d
248          out <- signal(out, states = states, filtered = TRUE)
249          names(out) <- c("thetahat", "V_theta")
250          if (d > 0) {
251          out$V_theta[,,1:d] <- Inf #diffuse phase
252          }
253        } else {
254          out <- signal(KFS(model = object, smoothing = "state",
255            maxiter = maxiter, expected = expected), states = states)
256          names(out) <- c("thetahat", "V_theta")
257        }
258
259      }
260
261      for (i in 1:p) {
262        pred[[i]] <- cbind(fit = out$thetahat[timespan, i] +
263            (if (object$distribution[i] == "poisson")  log(object$u[timespan, i]) else 0),
264          switch(interval, none = NULL,
265            out$thetahat[timespan, i] +
266              (if (object$distribution[i] == "poisson") log(object$u[timespan, i]) else 0) +
267              qnorm((1 - level)/2) * sqrt(out$V_theta[i, i, timespan]) %o%
268              c(1, -1)), se.fit = if (se.fit)
269                sqrt(out$V_theta[i, i, timespan]))
270        if (interval == "confidence")
271          colnames(pred[[i]])[2:3] <- c("lwr", "upr")
272      }
273      if (type == "response") {
274        if (se.fit) {
275          tmp <- which(colnames(pred[[1]]) == "se.fit")
276          for (i in 1:p) {
277            pred[[i]][, "se.fit"] <- switch(object$distribution[i],
278              gaussian = pred[[i]][,"se.fit"],
279              poisson = pred[[i]][, "se.fit"] * exp(pred[[i]][, 1]),
280              binomial = pred[[i]][, "se.fit"] *
281                (if (!prob) object$u[timespan, i] else 1) *
282                exp(pred[[i]][, 1])/(1 + exp(pred[[i]][, 1]))^2,
283              gamma = pred[[i]][, "se.fit"] * exp(pred[[i]][, 1]),
284              `negative binomial` = pred[[i]][, "se.fit"] * exp(pred[[i]][, 1]))
285            pred[[i]][, -tmp] <- switch(object$distribution[i],
286              gaussian = pred[[i]][, -tmp],
287              poisson = exp(pred[[i]][, -tmp]),
288              binomial = (if (!prob) object$u[timespan, i] else 1) *
289                exp(pred[[i]][, -tmp])/(1 + exp(pred[[i]][, -tmp])),
290              gamma = exp(pred[[i]][, -tmp]),
291              `negative binomial` = exp(pred[[i]][, -tmp]))
292          }
293        } else {
294          for (i in 1:p) pred[[i]] <- switch(object$distribution[i],
295            gaussian = pred[[i]],
296            poisson = exp(pred[[i]]),
297            binomial = (if (!prob) object$u[timespan, i] else 1) *
298              exp(pred[[i]])/(1 + exp(pred[[i]])),
299            gamma = exp(pred[[i]]),
300            `negative binomial` = exp(pred[[i]]))
301        }
302      }
303
304    } else {# with importance sampling
305      if (filtered) {
306        d <- KFS(approxSSM(object, maxiter = maxiter, expected = expected), smoothing = "none")$d
307      }
308
309      if (interval == "none") {
310        imp <- importanceSSM(object,
311          ifelse(identical(states, as.integer(1:m)), "signal", "states"),
312          nsim = nsim, antithetics = TRUE, maxiter = maxiter, filtered = filtered,
313          expected = expected)
314        nsim <- as.integer(4 * nsim)
315        if (!identical(states, as.integer(1:m))) {
316          imp$samples <- .Fortran(fzalpha, NAOK = TRUE, as.integer(dim(object$Z)[3] > 1),
317            object$Z, imp$samples, signal = array(0, c(n, p, nsim)),
318            p, m, n, nsim, length(states), states)$signal
319        }
320
321        w <- imp$weights/sum(imp$weights)
322        if (type == "response") {
323          for (i in 1:p) {
324            imp$samples[timespan, i, ] <- switch(object$distribution[i],
325              gaussian = imp$samples[timespan, i, ],
326              poisson = object$u[timespan, i] * exp(imp$samples[timespan, i, ]),
327              binomial = (if (!prob) object$u[timespan, i] else 1) *
328                exp(imp$samples[timespan, i, ])/(1 + exp(imp$samples[timespan, i, ])),
329              gamma = exp(imp$samples[timespan, i, ]),
330              `negative binomial` = exp(imp$samples[timespan, i, ]))
331          }
332        } else {
333          for (i in 1:p) if (object$distribution[i] == "poisson")
334            imp$samples[timespan, i, ] <- imp$samples[timespan, i, ] + log(object$u[timespan,
335              i])
336        }
337        varmean <- .Fortran(fvarmeanw, NAOK = TRUE, imp$samples[timespan, , , drop = FALSE], w,
338          p, length(timespan),
339          nsim, mean = array(0, c(length(timespan), p)),
340          var = array(0, c(length(timespan), p)), as.integer(se.fit))
341
342        if (se.fit) {
343          if (filtered && d > 0) {
344            varmean$var[1:d, ] <- Inf #diffuse phase
345          }
346          pred <- lapply(1:p, function(j) cbind(fit = varmean$mean[, j],
347            se.fit = sqrt(varmean$var[, j])))
348
349        } else {
350          pred <- lapply(1:p, function(j) varmean$mean[, j])
351        }
352      } else {
353        pred <- interval(object, interval = interval, level = level, type = type,
354          states = states, nsim = nsim, se.fit = se.fit, timespan = timespan,
355          prob = prob, maxiter = maxiter, filtered = filtered, expected = expected)
356        if (filtered && d > 0) {
357          for (i in 1:p) {
358            pred[[i]][1:d, "lwr"] <- -Inf
359            pred[[i]][1:d, "upr"] <- Inf #diffuse phase
360            if (se.fit) {
361              pred[[i]][1:d, "se.fit"] <- Inf #diffuse phase
362            }
363          }
364
365        }
366      }
367    }
368  }
369  names(pred) <- colnames(object$y)
370  pred <- lapply(pred, ts,end = endtime, frequency = frequency(object$y))
371  if (p == 1) pred <- pred[[1]]
372  pred
373}
374