# R/iarima.R In StMoMo: Stochastic Mortality Modelling

#### Documented in forecast.iarimaiarimasimulate.iarima

#' Fit independent arima series to a multivariate time series
#'
#' Fits independent arima series to \code{x}, a multivariate
#' time series.
#'
#'
#' @param x numeric matrix with a multivariate time series.
#' Series are arranged in rows with columns representing time.
#' @param order an optional matrix with one row per time series
#' specifying the ARIMA models: for the ith row the three
#' components \eqn{(p, d, q)} are the AR order, the degree of
#' differencing, and the MA order. If absent the arima models
#' @param include.constant an optional vector of logical values
#' indicating if the ARIMA model for the ith series should include a
#' constant value. The default is \code{TRUE}. This parameter is ignored
#' if \code{order} is \code{NULL}.
#'
#' @return an object of class \code{"iarima"} with components:
#' \item{models}{a list with the arima models fitted to each time
#' series.}
#' \item{x}{the original time series.}
#'
#' @details
#' The fitting of the ARIMA models for each time series is done with
#' \pkg{forecast}. See the latter function for further details on
#' input arguments \code{kt.order} and \code{kt.include.constant}.
#'
#' @export
iarima <- function(x, order = NULL, include.constant = TRUE, ...) {
x <- as.matrix(x)
if (ncol(x) == 1L)
x <- t(x)
N <- nrow(x)

if(!is.null(order)){
order <- as.matrix(order)
if (ncol(order) == 1L)
order <- t(order)
if (ncol(order) != 3)
stop("order must have three columns representing (p,d,q)")
if (nrow(order) == 1)
order <- matrix(order, nrow = N, ncol = 3, byrow = TRUE)
if (nrow(order) != N)
stop("number of rows in order does not match number of series.")
if (length(include.constant) == 1)
include.constant <- rep(include.constant, N)
if (length(include.constant) != N)
stop("Length of include.constant does not match number of series")
}

models <- list()
if(is.null(order)){
for (i in 1:N)  {
models[[i]] <- forecast::auto.arima(x[i, ], ...)
}
} else {
for (i in 1:N)  {
models[[i]] <- forecast::Arima(x[i, ], order = order[i, ],
include.constant = include.constant[i])
}
}
structure(list(models = models, x = x), class = "iarima")
}

#' Forecast independent arima series
#'
#' Returns forecasts and other information for a group of
#' independent arima series.
#'
#' @param object an object of class \code{"iarima"}.
#' @param h Number of periods for forecasting.
#' @param level confidence level for prediction intervals.
#' @param fan if \code{TRUE}, level is set to \code{seq(50, 99, by = 1)}.
#' This is suitable for fan plots.
#' @param ... other arguments.
#'
#' @return An object of class \code{"iarimaForecast"} with components:
#' \item{model}{a list containing information about the fitted arima models.}
#' \item{mean}{ array with the central forecast.}
#' \item{lower}{ three dimensional array with lower limits for prediction
#'  intervals.}
#' \item{upper}{ three dimensional array with upper limits for prediction
#'  intervals.}
#'  \item{level}{ the confidence values associated with the prediction
#'  intervals.}
#'  @export
forecast.iarima <- function(object, h = 10, level = c(80,95), fan = FALSE, ...) {

x <- object$x nYear <- ncol(x) N <- nrow(x) yearsFor <- (as.numeric(colnames(x)[nYear]) + 1):(as.numeric(colnames(x)[nYear]) + h) mean <- array(NA, dim = c(N, h), dimnames = list(rownames(x), yearsFor)) if (fan) level <- seq(51, 99, by = 3) else { if (min(level) > 0 & max(level) < 1) level <- 100 * level else if (min(level) < 0 | max(level) > 99.99) stop("Confidence limit out of range") } nconf <- length(level) lower <- upper <- array(NA, c(N, h, nconf), dimnames = list(rownames(x), yearsFor, paste(level, "%", sep = ""))) for(i in 1:N){ kt.for <- forecast(object$models[[i]], h = h, level = level)
mean[i, ] <- kt.for$mean lower[i, ,] <- kt.for$lower
upper[i, ,] <- kt.for$upper } structure(list(model = object, level = level, mean = mean, lower = lower, upper = upper), class = "iarimaForecast") } #' Simulate independent arima series #' #' Returns one simulated path of the group of independent #' arima series in \code{object}. #' #' @param object An object of class \code{"iarima"}. #' @param nsim number of periods for the simulated series. #' @param seed either \code{NULL} or an integer that will be used in a #' call to \code{\link{set.seed}} before simulating the time series. #' The default, \code{NULL} will not change the random generator state. #' @param ... other arguments. #' #' @export simulate.iarima <- function(object, nsim = 10, seed = NULL, ...) { #Hack to remove notes in CRAN check x <- NULL if (!exists(".Random.seed", envir = .GlobalEnv)) runif(1) if (is.null(seed)) RNGstate <- .Random.seed else { R.seed <- .Random.seed set.seed(seed) RNGstate <- structure(seed, kind = as.list(RNGkind())) on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv)) } sim <- t(simplify2array(lapply(object$models, FUN = simulate, nsim = nsim)))
rownames(sim) <- rownames(object$x) x <- object$x
nYear <- ncol(x)
yearsSim <- (as.numeric(colnames(x)[nYear]) + 1):(as.numeric(colnames(x)[nYear]) + nsim)