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