R/update.R

In robjhyndman/demography: Forecasting Mortality, Fertility, Migration and Population Data

#' Updating functional demographic models and coherent functional demographic models.
#'
#' \code{update.fmforecast()} updates \code{fdm} forecasts. The argument \code{object} is the output from \code{\link{forecast.fdm}} which has been subsequently modified with new coefficient forecasts. These new forecasts are used when re-calculating the forecast of the mortality or fertility rates, or net migration numbers.
#' \code{update.fmforecast2()} updates \code{fdmpr} forecasts. The argument \code{object} is the output from \code{\link{forecast.fdmpr}} which has been subsequently modified with new coefficient forecasts.
#'
#' @param object Output from either \code{\link{fdm}} or \code{\link{coherentfdm}}.
#' @param ... Extra arguments currently ignored.
#'
#' @return A list of the same class as \code{object}.
#' @author Rob J Hyndman.
#' @seealso \code{\link{forecast.fdm}}, \code{\link{forecast.fdmpr}}
#' @examples
#' \dontrun{
#' france.fit <- fdm(fr.mort, order = 2)
#' france.fcast <- forecast(france.fit, 50)
#' # Replace first coefficient model with ARIMA(0,1,2)+drift
#' france.fcast$coeff[[2]] <- forecast(Arima(france.fit$coeff[, 2],
#'   order = c(0, 1, 2), include.drift = TRUE
#' ), h = 50, level = 80)
#' france.fcast <- update(france.fcast)
#'
#' fr.short <- extract.years(fr.sm, 1950:2006)
#' fr.fit <- coherentfdm(fr.short)
#' fr.fcast <- forecast(fr.fit)
#' par(mfrow = c(1, 2))
#' plot(fr.fcast$male)
#' # Replace first coefficient model in product component with a damped ETS model:
#' fr.fcast$product$coeff[[2]] <- forecast(ets(fr.fit$product$coeff[, 2], damped = TRUE),
#'   h = 50, level = 80
#' )
#' fr.fcast <- update(fr.fcast)
#' plot(fr.fcast$male)
#' }
#' @keywords models
#' @name update
NULL

# Update fmforecast object
# Original object from forecast.fdm
# Assumed that the coefficient forecasts have been subsequently changed
# Object needs to be updated to reflect those changes.
#' @rdname update
#' @export
update.fmforecast <- function(object, ...) {
  if (!is.element("fmforecast", class(object))) {
    stop("object must be of class fmforecast")
  }
  h <- length(object$year)
  nb <- ncol(object$model$basis)
  adjust <- length(object$var$adj.factor) > 1
  objnames <- dimnames(object$rate[[1]])

  # Update in-sample fitted values and errors.
  fitted <- matrix(NA, length(object$model$year), nb)
  meanfcast <- varfcast <- matrix(NA, nrow = h, ncol = nb)
  qconf <- stats::qnorm(0.5 + object$coeff[[2]]$level[1] / 200)
  fitted[, 1] <- 1
  meanfcast[, 1] <- 1
  varfcast[, 1] <- 0
  for (i in 2:nb)
  {
    fitted[, i] <- fitted(object$coeff[[i]]$model)
    meanfcast[, i] <- object$coeff[[i]]$mean
    varfcast[, i] <- ((object$coeff[[i]]$upper - object$coeff[[i]]$lower) / (2 * qconf))^2
  }
  object$fitted$y <- object$model$basis %*% t(fitted)
  object$error$y <- object$model$y$y - object$fitted$y
  object$coeff.error <- object$model$coeff - fitted

  # Update point forecasts
  object$rate[[1]] <- object$model$basis %*% t(meanfcast)
  dimnames(object$rate[[1]]) <- objnames

  # Update forecast variances
  # Only model variance should have changed
  modelvar <- object$model$basis^2 %*% t(varfcast)
  totalvar <- sweep(modelvar, 1, object$var$error + object$var$mean, "+")
  if (adjust & nb > 1) {
    object$var$adj.factor <- rowMeans(object$error$y^2, na.rm = TRUE) / totalvar[, 1]
    totalvar <- sweep(totalvar, 1, object$var$adj.factor, "*")
  }
  # Add observational variance to total variance
  object$var$total <- sweep(totalvar, 1, object$var$observ, "+")

  # Update forecast intervals
  # Only parametric intervals computed here.
  tmp <- qconf * sqrt(object$var$total)
  object$rate$lower <- InvBoxCox(object$rate[[1]] - tmp, object$lambda)
  object$rate$upper <- InvBoxCox(object$rate[[1]] + tmp, object$lambda)
  object$rate[[1]] <- InvBoxCox(object$rate[[1]], object$lambda)
  if (object$type != "migration") {
    object$rate[[1]] <- pmax(object$rate[[1]], 0.000000001)
    object$rate$lower <- pmax(object$rate$lower, 0.000000001)
    object$rate$lower[is.na(object$rate$lower)] <- 0
    object$rate$upper <- pmax(object$rate$upper, 0.000000001)
  }

  # Return updated object
  return(object)
}


# Function to combine product and ratio forecasts
# object is output from forecast.fdmpr, but with modified forecasts
# This function simply recombines them again.
#' @rdname update
#' @export
update.fmforecast2 <- function(object, ...) {
  if (!is.element("fmforecast2", class(object))) {
    stop("object must be of class fmforecast2")
  }

  J <- length(object$ratio)
  ny <- length(object$ratio[[1]]$year)

  # GM model
  object$product <- update(object$product)

  # Obtain forecasts for each group
  is.mortality <- (object$product$type == "mortality")
  y <- as.numeric(is.mortality) # =1 for mortality and 0 for migration
  for (j in 1:J)
  {
    object$ratio[[j]] <- update(object$ratio[[j]])
    if (is.mortality) {
      object[[j]]$rate[[1]] <- object$product$rate$product * object$ratio[[j]]$rate[[1]]
    } else {
      object[[j]]$rate[[1]] <- object$product$rate$product + object$ratio[[j]]$rate[[1]]
    }
    if (is.mortality) {
      y <- y * object[[j]]$rate[[1]]
    } else {
      y <- y + object[[j]]$rate[[1]]
    }
  }

  # Adjust forecasts so they multiply appropriately.
  if (is.mortality) {
    y <- y^(1 / J) / object$product$rate$product
    for (j in 1:J) {
      object[[j]]$rate[[1]] <- object[[j]]$rate[[1]] / y
    }
  } else {
    y <- y / J - object$product$rate$product
    for (j in 1:J) {
      object[[j]]$rate[[1]] <- object[[j]]$rate[[1]] - y
    }
  }
  # Variance of forecasts
  qconf <- 2 * stats::qnorm(0.5 + object$product$coeff[[1]]$level / 200)
  for (j in 1:J)
  {
    vartotal <- object$product$var$total + object$ratio[[j]]$var$total
    tmp <- qconf * sqrt(vartotal)
    object[[j]]$rate$lower <- InvBoxCox(BoxCox(object[[j]]$rate[[1]], object$product$lambda) - tmp, object$product$lambda)
    object[[j]]$rate$upper <- InvBoxCox(BoxCox(object[[j]]$rate[[1]], object$product$lambda) + tmp, object$product$lambda)
  }

  return(object)
}