R/methods.R

In dfms: Dynamic Factor Models

#' @srrstats {TS2.5} *Incorporate a system to ensure that both row and column orders follow the same ordering as the underlying time series data. This may, for example, be done by including the `index` attribute of the time series data as an attribute of the auto-covariance matrix.*
#' @srrstats {TS2.6} *Where applicable, auto-covariance matrices should also include specification of appropriate units.*
#' @srrstats {TS3.0} *Provide tests to demonstrate at least one case in which errors widen appropriately with forecast horizon.*
#' @srrstats {TS3.1} *If possible, provide at least one test which violates TS3.0*
#' -> currently I don't forecast the covariance matrices. This could be implemented in the future.
#' @srrstats {TS3.2} *Document the general drivers of forecast errors or horizons, as demonstrated via the particular cases of TS3.0 and TS3.1*
#' @srrstats {TS3.3} *Either:*
#' @srrstats {TS3.3a} *Document, preferable via an example, how to trim forecast values based on a specified error margin or equivalent; or*
#' @srrstats {TS3.3b} *Provide an explicit mechanism to trim forecast values to a specified error margin, either via an explicit post-processing function, or via an input parameter to a primary analytic function.*
#' @srrstats {TS4.0} *Return values should either:*
#' @srrstats {TS4.0a} *Be in same class as input data, for example by using the [`tsbox` package](https://www.tsbox.help/) to re-convert from standard internal format (see 1.4, above); or*
#' @srrstats {TS4.0b} *Be in a unique, preferably class-defined, format.*
#' @srrstats {TS4.1} *Any units included as attributes of input data should also be included within return values.*
#' @srrstats {TS4.2} *The type and class of all return values should be explicitly documented.*
#' @srrstats {TS4.3} *Return values should explicitly include all appropriate units and/or time scales*
#' @srrstats {TS4.4} *Document the effect of any such transformations on forecast data, including potential effects on both first- and second-order estimates.*
#' @srrstats {TS4.5} *In decreasing order of preference, either:*
#' @srrstats {TS4.5a} *Provide explicit routines or options to back-transform data commensurate with original, non-stationary input data*
#' @srrstats {TS4.5b} *Demonstrate how data may be back-transformed to a form commensurate with original, non-stationary input data.*
#' @srrstats {TS4.5c} *Document associated limitations on forecast values*
#' @srrstats {TS4.6} *Time Series Software which implements or otherwise enables forecasting should return either:*
#' @srrstats {TS4.6a} *A distribution object, for example via one of the many packages described in the CRAN Task View on [Probability Distributions](https://cran.r-project.org/web/views/Distributions.html) (or the new [`distributional` package](https://pkg.mitchelloharawild.com/distributional/) as used in the [`fable` package](https://fable.tidyverts.org) for time-series forecasting).*
#' @srrstats {TS4.6b} *For each variable to be forecast, predicted values equivalent to first- and second-order moments (for example, mean and standard error values).*
#' @srrstats {TS4.6c} *Some more general indication of error associated with forecast estimates.*
#' @srrstats {TS4.7} *Ensure that forecast (modelled) values are clearly distinguished from observed (model or input) values, either (in this case in no order of preference) by*
#' @srrstats {TS4.7a} *Returning forecast values alone*
#' @srrstats {TS4.7b} *Returning distinct list items for model and forecast values*
#' @srrstats {TS4.7c} *Combining model and forecast values into a single return object with an appropriate additional column clearly distinguishing the two kinds of data.*
#' @srrstats {TS5.0} *Implement default `plot` methods for any implemented class system.*
#' @srrstats {TS5.1} *When representing results in temporal domain(s), ensure that one axis is clearly labelled "time" (or equivalent), with continuous units.*
#' @srrstats {TS5.2} *Default to placing the "time" (or equivalent) variable on the horizontal axis.*
#' @srrstats {TS5.3} *Ensure that units of the time, frequency, or index variable are printed by default on the axis.*
#' @srrstats {TS5.5} *Provide options to determine whether plots of data with missing values should generate continuous or broken lines.*
#' @srrstats {TS5.6} *By default indicate distributional limits of forecast on plot*
#' @srrstats {TS5.7} *By default include model (input) values in plot, as well as forecast (output) values*
#' @srrstats {TS5.8} *By default provide clear visual distinction between model (input) values and forecast (output) values.*


#' @name summary.dfm
#' @aliases print.dfm
#' @aliases summary.dfm
#' @aliases print.dfm_summary
#'
#' @title DFM Summary Methods
#'
#' @description Summary and print methods for class 'dfm'. \code{print.dfm} just prints basic model information and the factor transition matrix \eqn{\textbf{A}}{A}, \code{coef.dfm} returns \eqn{\textbf{A}}{A} and \eqn{\textbf{C}}{C} in a plain list, whereas
#' \code{summary.dfm} returns all system matrices and additional residual and goodness of fit statistics---with a print method allowing full or compact printout.
#'
#' @param x,object an object class 'dfm'.
#' @param digits integer. The number of digits to print out.
#' @param \dots not used.
#'
#' @seealso \link{dfms-package}
#' @importFrom collapse qsu frange
#' @export
print.dfm <- function(x, digits = 4L, ...) {

  X <- x$X_imp
  A <- x$A
  r <- dim(A)[1L]
  p <- dim(A)[2L]/r
  if(length(qv <- x$quarterly.vars)) {
    cat("Mixed Frequency Dynamic Factor Model\nn = ", dim(X)[2L], ", nm = ", dim(X)[2L] - length(qv), ", nq = ", length(qv), ", T = ", dim(X)[1L], ", r = ", r, ", p = ", p,
        "\n%NA = ", round(sum(attr(X, "missing"))/prod(dim(X))*100, digits), ", %NAm = ", round(sum(attr(X, "missing")[, -ckmatch(qv, colnames(X))])/(nrow(X)*(ncol(X)-length(qv)))*100, digits), "\n", sep = "")
  } else {
  cat("Dynamic Factor Model: n = ", dim(X)[2L], ", T = ", dim(X)[1L], ", r = ", r, ", p = ", p, ", %NA = ",
      if(x$anyNA) round(sum(attr(X, "missing"))/prod(dim(X))*100, digits) else 0,"\n", sep = "")
  }
  if(length(x$rho)) cat("   with AR(1) errors: mean(abs(rho)) =", round(mean(abs(x$rho)), 3), "\n")
  fnam <- paste0("f", seq_len(r))
  cat("\nFactor Transition Matrix [A]\n")
  print(round(A, digits))
  return(invisible(x))
}

#' @rdname summary.dfm
#' @export
coef.dfm <- function(object, ...) list(A = object$A, C = object$C)

#' @rdname summary.dfm
#' @export
logLik.dfm <- function(object, ...) object$loglik[length(object$loglik)]

#' @rdname summary.dfm
#' @param method character. The factor estimates to use: one of \code{"qml"}, \code{"2s"} or \code{"pca"}.
#' @param \dots not used.
#' @return Summary information following a dynamic factor model estimation. \code{coef()} returns \eqn{\textbf{A}}{A} and \eqn{\textbf{C}}{C}.
#' @importFrom stats cov
#' @importFrom collapse pwcov
#' @export
summary.dfm <- function(object, method = switch(object$em.method, none = "2s", "qml"), ...) {

  X <- object$X_imp
  Fa <- switch(tolower(method), pca = object$F_pca, `2s` = object$F_2s, qml = object$F_qml, stop("Unkown method", method))
  A <- object$A
  r <- dim(A)[1L]
  p <- dim(A)[2L] / r
  C <- object$C
  idio_ar1 <- !is.null(object$rho)
  res <- if(idio_ar1) object[["e"]] else  X - tcrossprod(Fa, C)
  anymissing <- object$anyNA
  if(!idio_ar1 && anymissing) res[attr(X, "missing")] <- NA
  rescov <- pwcov(res, use = if(!idio_ar1 && anymissing) "pairwise.complete.obs" else "everything", P = TRUE)
  ACF <- if(idio_ar1) object$rho else AC1(res, anymissing)
  R2 <- 1 - diag(rescov[,, 1L])
  summ <- list(info = c(n = dim(X)[2L], T = dim(X)[1L], r = r, p = p, nq = length(object$quarterly.vars),
                        `%NA` = if(anymissing) sum(attr(X, "missing")) / prod(dim(X)) * 100 else 0,
                        `%NAm` = if(length(object$quarterly.vars)) sum(attr(X, "missing")[, -ckmatch(object$quarterly.vars, colnames(X))])/(nrow(X)*(ncol(X)-length(object$quarterly.vars)))*100 else NA),
               call = object$call,
               idio_ar1 = idio_ar1,
               F_stats = msum(Fa),
               A = A,
               F_cov = pwcov(Fa, P = TRUE),
               Q = object$Q,
               C = C,
               R_diag = diag(object$R),
               res_cov = rescov,
               res_ACF = ACF,
               res_ACF_stats = msum(ACF),
               R2 = R2,
               R2_stats = msum(R2))
  class(summ) <- "dfm_summary"
  return(summ)
}

#' @rdname summary.dfm
#' @param compact integer. Display a more compact printout: \code{0} prints everything, \code{1} omits the observation matrix \eqn{\textbf{C}}{C} and residual covariance matrix \code{cov(resid(model))}, and \code{2} omits all disaggregated information on the input data. Sensible default are chosen for different sizes of the input dataset so as to limit large printouts.
#' @param \dots not used.
#'
#' @examples
#' mod <- DFM(diff(BM14_Q), 2, 3)
#' print(mod)
#' summary(mod)
#'
#' @export
print.dfm_summary <- function(x,
                              digits = 4L,
                              compact = sum(x$info["n"] > 15, x$info["n"] > 40), ...) {

  inf <- as.integer(x$info[1:4])
  if(x$info["nq"] > 0) {
    cat("Mixed Frequency Dynamic Factor Model\nn = ", inf[1L], ", nm = ", inf[1L] - x$info["nq"], ", nq = ", x$info["nq"], ", T = ", inf[2L], ", r = ", inf[3L], ", p = ", inf[4L],
        "\n%NA = ", round(x$info["%NA"], digits), ", %NAm = ", round(x$info["%NAm"], digits), "\n", sep = "")
  } else {
    cat("Dynamic Factor Model: n = ", inf[1L], ", T = ", inf[2L], ", r = ", inf[3L], ", p = ", inf[4L],
        ", %NA = ", round(x$info["%NA"], digits), "\n", sep = "")
  }
  if(x$idio_ar1) cat("   with AR(1) errors: mean(abs(rho)) =", round(mean(abs(x$res_ACF)), 3), "\n")
  cat("\nCall: ", deparse(x$call))
  # cat("\nModel: ", ))
  cat("\n\nSummary Statistics of Factors [F]\n")
  print(x$F_stats, digits)
  cat("\nFactor Transition Matrix [A]\n")
  print(x$A, digits = digits)
  cat("\nFactor Covariance Matrix [cov(F)]\n")
  print(x$F_cov, digits)
  cat("\nFactor Transition Error Covariance Matrix [Q]\n")
  print(round(x$Q, digits))
  if(compact == 0L) {
  cat("\nObservation Matrix [C]\n")
  print(round(x$C, digits))
  }
  if(compact < 2L) {
  cat("\nObservation Error Covariance Matrix [diag(R) - Restricted]\n")
  # cat("\n Estimated Diagonal (DFM Assumes R is Diagonal)\n")
  print(round(x$R_diag, digits))
  }
  if(compact == 0L) {
  cat("\nObservation Residual Covariance Matrix [cov(resid(DFM))]\n")
  print(x$res_cov, digits)
  }
  if(compact < 2L) {
  cat("\nResidual AR(1) Serial Correlation\n")
  print(x$res_ACF, digits) # TODO: Add P-Value
  }
  cat("\nSummary of Residual AR(1) Serial Correlations\n")
  print(x$res_ACF_stats, digits)
  if(compact < 2L) {
  cat("\nGoodness of Fit: R-Squared\n")
  print(x$R2, digits)
  }
  cat("\nSummary of Individual R-Squared's\n")
  print(x$R2_stats, digits)
  return(invisible(x))
}


#' Plot DFM
#' @param x an object class 'dfm'.
#' @param method character. The factor estimates to use: one of \code{"qml"}, \code{"2s"}, \code{"pca"} or \code{"all"} to plot all estimates.
#' @param type character. The type of plot: \code{"joint"}, \code{"individual"} or \code{"residual"}.
#' @param scale.factors logical. Standardize factor estimates, this usually improves the plot since the factor estimates corresponding to the greatest PCA eigenvalues tend to have a greater variance than the data.
#' @param \dots for \code{plot.dfm}: further arguments to \code{\link{plot}}, \code{\link{ts.plot}}, or \code{\link{boxplot}}, depending on the \code{type} of plot. For \code{screeplot.dfm}: further arguments to \code{\link{screeplot.ICr}}.
#' @returns Nothing.
#' @seealso \link{dfms-package}
#' @examples \donttest{
#' # Fit DFM with 3 factors and 3 lags in the transition equation
#' mod <- DFM(diff(BM14_M), r = 3, p = 3)
#' plot(mod)
#' plot(mod, type = "individual", method = "all")
#' plot(mod, type = "residual")
#' }
#'
#' @importFrom graphics boxplot axis box mtext plot.default
#' @importFrom collapse unlist2d ckmatch na_rm seq_row
#' @export
plot.dfm <- function(x,
                     method = switch(x$em.method, none = "2s", "qml"),
                     type = c("joint", "individual", "residual"),
                     scale.factors = TRUE, ...) {
  Fa <- switch(tolower(method[1L]),
              all = cbind(x$F_pca, setColnames(x$F_2s, paste("TwoStep", colnames(x$F_2s))),
                          if(length(x$F_qml)) setColnames(x$F_qml, paste("QML", colnames(x$F_qml))) else NULL),
              pca = x$F_pca, `2s` = x$F_2s, qml = x$F_qml, stop("Unknown method:", method[1L]))

  nf <- dim(Fa)[2L]
  allests <- tolower(method[1L]) == "all"
  dots <- list(...)

  switch(tolower(type[1L]),
    joint = {
      Xr <- frange(x$X_imp)
      if(scale.factors) Fa <- fscale(Fa)
      Fr <- frange(Fa)
      ts.plot(x$X_imp, col = "grey85", ylim = c(min(Xr[1L], Fr[1L]), max(Xr[2L], Fr[2L])),
              ylab = if(is.null(dots$ylab)) "Value" else dots$ylab,
              main = if(is.null(dots$main)) "Standardized Series and Factor Estimates" else dots$main, ...)
      cols <- rainbow(nf)
      for (i in seq_len(nf)) lines(Fa[, i], col = cols[i])
      legend("topleft", colnames(Fa), col = cols, lty = 1, bty = "n", ncol = if(allests) 3L else 1L)
    },
    individual = {
      # if(allests) {
      if(scale.factors) Fa <- fscale(Fa)
      qml <- !is.null(x$F_qml)
      if(allests) nf <- nf / (2L + qml)

      # Extracted from plot.ts()...
      cex.lab = par("cex.lab")
      col.lab = par("col.lab")
      font.lab = par("font.lab")
      oldpar <- par(mar = c(0, 5.1, 0, 2.1), oma = c(6, 0, 5, 0), mfrow = c(nf, 1L))
      on.exit(par(oldpar))

      for(i in seq_len(nf)) {
        plot.default(Fa[, i], axes = FALSE, xlab = "", ylab = "", type = "n")
        lines(Fa[, i], type = 'l', col = if(allests) "red" else "black", ...)
        if(allests) {
          lines(Fa[, i + nf], type = 'l', col = "orange", ...)
          if(qml) lines(Fa[, i + 2L * nf], type = 'l', col = "blue", ...)
          if(i == 1L) legend("topleft", c("PCA", "TwoStep", if(qml) "QML"),
                             col = c("red", "orange", "blue"), lty = 1, bty = "n")
        }
        box(...)
        axis(2, xpd = NA, ...)
        if(i == nf) axis(1, xpd = NA, ...)
        mtext(paste("f", i), 2, line = 3, cex = cex.lab, col = col.lab, font = font.lab, ...)
        if(i == nf) mtext(if(is.null(dots$xlab)) "Time" else dots$xlab, side = 1, line = 3, cex = cex.lab, col = col.lab, font = font.lab, ...)
      }
      par(mfrow = c(1, 1)) # on.exit above takes care of changes to parameters.
      mtext(if(is.null(dots$main)) paste(if(scale.factors) "Standardized", "Factor Estimates") else dots$main,
            side = 3, line = 3, cex = par("cex.main"), font = par("font.main"), col = par("col.main"), ...)
      # } else {
      #   oldpar <- par(mfrow = c(nf, 1L))
      #   on.exit(par(oldpar))
      #   cnF <- colnames(Fa)
      #   for (i in seq_len(nf)) plot(Fa[, i], type = 'l', main = cnF[i], ylab = "Value",
      #                               xlab = if(i == nf) "Time" else "" , ...)
      # }
    },
    residual = {
      if(allests) stop("Need to choose a specific method for residual plots")
      oldpar <- par(mar = c(11.5, 4.1, 4.1, 2.1))
      on.exit(par(oldpar))
      # Use residuals method which handles MQ and idio.ar1 properly
      res <- residuals.dfm(x, method = method[1L], standardized = TRUE, na.keep = TRUE)
      boxplot(res, main = if(is.null(dots$main)) "Residuals by input variable" else dots$main, las = 2, ...)
    },
    stop("Unknown plot type: ", type[1L])
  )
}


#' Extract Factor Estimates in a Data Frame
#' @param x an object class 'dfm'.
#' @param method character. The factor estimates to use: any of \code{"qml"}, \code{"2s"}, \code{"pca"} (multiple can be supplied) or \code{"all"} for all estimates.
#' @param pivot character. The orientation of the frame: \code{"long"}, \code{"wide.factor"} or \code{"wide.method"}, \code{"wide"} or \code{"t.wide"}.
#' @param time a vector identifying the time dimension, or \code{NULL} to omit a time variable.
#' @param stringsAsFactors make factors from method and factor identifiers. Same as option to \code{\link{as.data.frame.table}}.
#' @param \dots not used.
#'
#' @return A data frame of factor estimates.
#' @seealso \link{dfms-package}
#'
#' @examples \donttest{
#' library(xts)
#' # Fit DFM with 3 factors and 3 lags in the transition equation
#' mod <- DFM(diff(BM14_M), r = 3, p = 3)
#'
#' # Taking a single estimate:
#' print(head(as.data.frame(mod, method = "qml")))
#' print(head(as.data.frame(mod, method = "qml", pivot = "wide")))
#'
#' # Adding a proper time variable
#' time <- index(BM14_M)[-1L]
#' print(head(as.data.frame(mod, method = "qml", time = time)))
#'
#' # All estimates: different pivoting methods
#' for (pv in c("long", "wide.factor", "wide.method", "wide", "t.wide")) {
#'    cat("\npivot = ", pv, "\n")
#'    print(head(as.data.frame(mod, pivot = pv, time = time), 3))
#' }
#' }
#'
#' @importFrom collapse ckmatch na_rm seq_row t_list unattrib
#' @importFrom stats setNames
#' @export
as.data.frame.dfm <- function(x, ...,
                              method = "all",
                              pivot = c("long", "wide.factor", "wide.method", "wide", "t.wide"),
                              time = seq_row(x$F_pca),
                              stringsAsFactors = TRUE) {

  estm <- c(PCA = "pca", TwoStep = "2s", QML = "qml")
  if(length(method) > 1L || method != "all")
     estm <- estm[ckmatch(method, estm, e = "Unknown method:")]
  estlist <- x[paste0("F_", estm)]
  names(estlist) <- names(estm)
  estlist <- na_rm(estlist) # Also removes NULL elements

  nam <- names(estlist)
  m <- length(estlist)
  TT <- nrow(estlist[[1L]])
  r <- ncol(estlist[[1L]])

  if(!is.null(time) && length(time) != TT) {
    if(length(x$rm.rows)) time <- time[-x$rm.rows]
    if(length(time) != TT) stop(sprintf("time must be a length %s vector or NULL", TT))
  }

  res <- switch(tolower(pivot[1L]),
    long = list(Method = if(stringsAsFactors) setAttrib(rep(1:m, each = TT*r), list(levels = nam, class = "factor")) else rep(nam, each = TT*r),
                Factor = if(stringsAsFactors) setAttrib(rep(1:r, times = m, each = TT), list(levels = paste0("f", 1:r), class = "factor")) else rep(paste0("f", 1:r), times = m, each = TT),
                Time = if(length(time)) rep(time, times = m*r) else NULL,
                Value = unlist(estlist, use.names = FALSE)),
    wide.factor = c(list(Method = if(stringsAsFactors) setAttrib(rep(1:m, each = TT), list(levels = nam, class = "factor")) else rep(nam, each = TT),
                         Time = if(length(time)) rep(time, times = m) else NULL),
                    setNames(lapply(t_list(unattrib(lapply(estlist, mctl))), unlist, FALSE, FALSE), paste0("f", 1:r))),
    wide.method = c(list(Factor = if(stringsAsFactors) setAttrib(rep(1:r, each = TT), list(levels = paste0("f", 1:r), class = "factor")) else rep(paste0("f", 1:r), each = TT),
                         Time = if(length(time)) rep(time, times = r) else NULL),
                    lapply(estlist, unattrib)),
    # If only one method, do not do combine names e.g. "QML_f1"? -> most of the time people just want a simple frame like this...
    wide = c(list(Time = time), setNames(unlist(lapply(estlist, mctl), FALSE, FALSE), if(length(nam) == 1L) paste0("f", 1:r) else outer(paste0("f", 1:r), nam, paste, sep = "_"))),
    t.wide = c(list(Time = time), setNames(unlist(t_list(lapply(estlist, mctl)), FALSE, FALSE), if(length(nam) == 1L) paste0("f", 1:r) else t(outer(paste0("f", 1:r), nam, paste, sep = "_")))),
    stop("Unknown pivot option:", pivot[1L])
  )

  if(is.null(time)) res <- na_rm(res)
  attr(res, "methods") <- estm
  attr(res, "row.names") <- .set_row_names(length(res[[1L]]))
  class(res) <- "data.frame"
  return(res)
}

predict_dfm_core <- function(object, method, use.full.state = TRUE) {
  method <- tolower(method)
  if(isTRUE(use.full.state) && method != "pca" && length(object$ss_full)) {
    ss_full <- object$ss_full
    if(!is.null(ss_full) && !is.null(ss_full$C) && !is.null(ss_full$F_smooth)) {
      if(nrow(ss_full$F_smooth) == nrow(object$X_imp) && nrow(ss_full$C) == ncol(object$X_imp)) {
        res <- tcrossprod(ss_full$F_smooth, ss_full$C)
        dimnames(res) <- dimnames(object$X_imp)
        return(res)
      }
    }
  }
  Fa <- switch(method,
               pca = object$F_pca, `2s` = object$F_2s, qml = object$F_qml,
               stop("Unkown method", method))
  if(is.null(object$quarterly.vars)) return(tcrossprod(Fa, object$C))
  qind <- ckmatch(object$quarterly.vars, dimnames(object$C)[[1L]])
  res_m <- tcrossprod(Fa, object$C[-qind,, drop = FALSE])
  Fa_lags <- flag(Fa, 0:4, fill = 0, stubs = FALSE)
  Cq_lags <- object$C[qind, rep(1:ncol(Fa), each = 5), drop = FALSE] %r*% rep(c(1, 2, 3, 2, 1), ncol(Fa))
  res_q <- tcrossprod(Fa_lags, Cq_lags)
  cbind(res_m, res_q)
}

#' @name residuals.dfm
#' @aliases residuals.dfm
#' @aliases resid.dfm
#' @aliases fitted.dfm
#'
#' @title DFM Residuals and Fitted Values
#' @description The residuals \eqn{\textbf{e}_t = \textbf{x}_t - \textbf{C} \textbf{F}_t}{e(t) = x(t) - C F(t)} or fitted values \eqn{\textbf{C} \textbf{F}_t}{C F(t)} of the DFM observation equation.
#'
#' @param object an object of class 'dfm'.
#' @param method character. The factor estimates to use: one of \code{"qml"}, \code{"2s"} or \code{"pca"}.
#' @param orig.format logical. \code{TRUE} returns residuals/fitted values in a data format similar to \code{X}.
#' @param standardized logical. \code{FALSE} will put residuals/fitted values on the original data scale.
#' @param na.keep logical. \code{TRUE} inserts missing values where \code{X} is missing (default \code{TRUE} as residuals/fitted values are only defined for observed data). \code{FALSE} returns the raw prediction, which can be used to interpolate data based on the DFM. For residuals, \code{FALSE} returns the difference between the prediction and the initial imputed version of \code{X} use for PCA to initialize the Kalman Filter.
#' @param use.full.state logical. Use the full state-space (if available) for fitted values and residuals. This includes idiosyncratic state components when \code{idio.ar1 = TRUE}, so fitted values reflect the full observation equation and residuals measure what is left after both factor and idiosyncratic components. Set to \code{FALSE} to obtain factor-only fitted values and residuals. Falls back to the compact form if unavailable or if \code{method = "pca"}.
#' @param \dots not used.
#'
#' @return A matrix of DFM residuals or fitted values. If \code{orig.format = TRUE} the format may be different, e.g. a data frame.
#'
#' @seealso \link{dfms-package}
#'
#' @examples \donttest{
#' library(xts)
#' # Fit DFM with 3 factors and 3 lags in the transition equation
#' mod <- DFM(diff(BM14_M), r = 3, p = 3)
#'
#' # Residuals
#' head(resid(mod))
#' plot(resid(mod, orig.format = TRUE)) # this is an xts object
#'
#' # Fitted values
#' head(fitted(mod))
#' head(fitted(mod, orig.format = TRUE)) # this is an xts object
#' }
#'
#' @importFrom collapse TRA.matrix mctl setAttrib pad
#' @export
residuals.dfm <- function(object,
                          method = switch(object$em.method, none = "2s", "qml"),
                          orig.format = FALSE,
                          standardized = FALSE,
                          na.keep = TRUE,
                          use.full.state = TRUE, ...) {
  X <- object$X_imp
  if(!(standardized && length(object[["e"]]))) {
    X_pred <- predict_dfm_core(object, method, use.full.state = use.full.state)
    if(!standardized) {  # TODO: What if AR(1) resid available?
      stats <- attr(X, "stats")
      X_pred <- unscale(X_pred, stats)
      res <- unscale(X, stats) - X_pred
    } else res <- X - X_pred
  } else res <- object[["e"]]
  if(na.keep && object$anyNA) res[attr(X, "missing")] <- NA
  if(orig.format) {
    if(length(object$rm.rows)) res <- pad(res, object$rm.rows, method = "vpos")
    if(attr(X, "is.list")) res <- mctl(res)
    return(setAttrib(res, attr(X, "attributes")))
  }
  return(qM(res))
}

#' @rdname residuals.dfm
#' @export
fitted.dfm <- function(object,
                       method = switch(object$em.method, none = "2s", "qml"),
                       orig.format = FALSE,
                       standardized = FALSE,
                       na.keep = TRUE,
                       use.full.state = TRUE, ...) {
  X <- object$X_imp
  res <- predict_dfm_core(object, method, use.full.state = use.full.state)
  if(!standardized) res <- unscale(res, attr(X, "stats"))
  if(na.keep && object$anyNA) res[attr(X, "missing")] <- NA
  if(orig.format) {
    if(length(object$rm.rows)) res <- pad(res, object$rm.rows, method = "vpos")
    if(attr(X, "is.list")) res <- mctl(res)
    return(setAttrib(res, attr(X, "attributes")))
  }
  return(qM(res))
}

#' News Decomposition
#'
#' Compute the Banbura and Modugno (2014) news decomposition of forecast updates.
#' Given an old vintage and an updated vintage, the function decomposes the
#' forecast revision at \code{t.fcst} into contributions from new releases.
#'
#' Let \eqn{y_t^{old}} and \eqn{y_t^{new}} be the old and new forecasts of a target
#' series at \eqn{t = t_{fcst}}. For each new release \eqn{i} (a previously missing
#' observation that becomes observed), the innovation is
#' \deqn{\nu_i = x_i^{new} - \hat{x}_i^{old},}
#' where \eqn{\hat{x}_i^{old}} is the smoothed estimate from the old vintage.
#' The revision is decomposed as
#' \deqn{y_t^{new} - y_t^{old} = \sum_i g_i \nu_i,}
#' with gain weights computed from Kalman smoother covariances:
#' \deqn{g = \sigma_y C_y P_1 P_2^{-1}.}
#' Here \eqn{\sigma_y} is the target series standard deviation, \eqn{C_y} is the
#' loading row for the target series, \eqn{P_1} collects cross-covariances between
#' the target and each news item, and \eqn{P_2} is the covariance matrix of the
#' news items (including measurement error where appropriate). See Section 2.3 and
#' Appendix D in Banbura and Modugno (2014).
#'
#' The function uses the system matrices and scaling from the new vintage. The old
#' data are re-standardized to the new-vintage scale before smoothing so that
#' innovations and gains are computed on a consistent scale. Set
#' \code{standardized = FALSE} to report results on the original data scale.
#'
#' @note This implementation is translated from the original MATLAB codes and is
#' consistent with the BM2014 news decomposition formulas.
#' If the model was estimated with \code{max.missing < 1} and
#' \code{na.rm.method = "LE"} in \code{\link{tsnarmimp}} (called by \code{DFM()}), leading or trailing rows with many missing values
#' may be removed by \code{DFM()}. If old and new vintages are both dfm objects, and they drop different rows,
#' then \code{t.fcst} can become out of bounds. When \code{comparison} is provided
#' as raw data, \code{news()} drops \code{object$rm.rows} from the new dataset (if present) and
#' forces \code{max.missing = 1} for the re-estimation call to keep row alignment.
#' To avoid issues, estimate both vintages with \code{max.missing = 1}.
#' For mixed-frequency or idiosyncratic AR(1) models, \code{news()} relies on the full
#' state-space matrices stored in \code{dfm$ss_full}.
#'
#' @param object a \code{dfm} object for the old vintage.
#' @param comparison a \code{dfm} object or a new dataset for the updated vintage.
#' @param t.fcst integer. Forecast target time index.
#' @param target.vars Integer or character identifying target variables. Defaults to all variables.
#' @param series optional character vector for naming variables.
#' @param standardized logical. Return results on standardized scale?
#' @param \dots not used.
#' @return For a single target, a \code{dfm_news} object with elements:
#' \itemize{
#' \item \code{y_old}: old forecast for the target variable at \code{t.fcst}.
#' \item \code{y_new}: new forecast for the target variable at \code{t.fcst}.
#' \item \code{news_df}: data frame with one row per series and columns:
#' \itemize{
#' \item \code{series}: series name.
#' \item \code{actual}: actual release (if any).
#' \item \code{forecast}: old-vintage forecast of the release.
#' \item \code{news}: total innovation for the series on the output scale. If there is a
#' single release, \code{news} equals \code{actual - forecast}. With multiple releases,
#' \code{news} aggregates those innovations for the series.
#' \item \code{gain}: effective weight on \code{news} such that \code{impact = news * gain}
#' (on the output scale).
#' \item \code{gain_std}: effective weight on the standardized innovations.
#' \item \code{impact}: contribution of the series to the target revision.
#' }
#' }
#' If \code{target.vars} selects multiple targets, a \code{dfm_news_list} object is returned,
#' where each element is a \code{dfm_news} object and list names correspond to targets.
#'
#' @references
#' Banbura, M., & Modugno, M. (2014). Maximum likelihood estimation of factor
#' models on datasets with arbitrary pattern of missing data.
#' *Journal of Applied Econometrics, 29*(1), 133-160.
#'
#' @seealso \link{dfms-package}
#'
#' @examples \donttest{
#' # (1) Monthly DFM example
#' X <- collapse::qM(BM14_M)[, BM14_Models$medium[BM14_Models$freq == "M"]]
#' X_old <- X
#' # Creating earlier vintage
#' X_old[nrow(X) - 1, sample(which(is.finite(X[nrow(X) - 1, ]) & is.na(X[nrow(X), ])), 5)] <- NA
#' X_old[nrow(X), sample(which(is.finite(X[nrow(X), ])), 5)] <- NA
#' # Estimating DFM
#' dfm <- DFM(X_old, r = 2, p = 2, em.method = "none")
#' # News computation (second DFM fit internally with same settings and rows)
#' res <- news(dfm, X, target.vars = c("ip_tot_cstr", "orders", "urx"))
#' # See results
#' print(res)
#' head(res$news_df)
#'
#' # (2) MQ nowcast of GDP (idio.ar1 = FALSE for speed)
#' library(magrittr)
#' library(xts)
#' # Creating MQ dataset
#' BM14 <- merge(BM14_M, BM14_Q)
#' BM14[, BM14_Models$log_trans] %<>% log()
#' BM14[, BM14_Models$freq == "M"] %<>% diff()
#' BM14[, BM14_Models$freq == "Q"] %<>% diff(3)
#' X <- BM14[-1, BM14_Models$small]
#' quarterly.vars <- BM14_Models$series[BM14_Models$small & BM14_Models$freq == "Q"]
#' # Creating earlier vintage
#' X_old <- X
#' X_old[355, c("ip_tot_cstr", "new_cars")] <- NA
#' X_old[356, c("new_cars", "pms_pmi", "euro325", "capacity")] <- NA
#' # Estimating DFM
#' dfm <- DFM(X_old, r = 2, p = 2, quarterly.vars = quarterly.vars, max.missing = 1)
#' # News computation (second DFM fit internally with same settings and rows)
#' res_mq <- news(dfm, X, t.fcst = 356, target.vars = "gdp")
#' # See results
#' print(res_mq)
#' head(res_mq$news_df)
#' }
#'
#' @export
news <- function(object, ...) UseMethod("news")

#' @rdname news
#'
#' @export
news.dfm <- function(object,
                     comparison,
                     t.fcst = nrow(object$X_imp),
                     target.vars = NULL,
                     series = NULL,
                     standardized = FALSE, ...) {
  if(inherits(comparison, "dfm")) {
    dfm_new <- comparison
  } else {
    cl <- object$call
    if(is.null(cl)) stop("dfm object has no call to re-estimate comparison")
    if(length(object$rm.rows)) comparison <- comparison[-object$rm.rows, , drop = FALSE]
    cl$max.missing <- 1
    if(is.null(cl$X)) cl[[2L]] <- comparison else cl$X <- comparison
    dfm_new <- eval(cl, envir = parent.frame())
    if(!inherits(dfm_new, "dfm")) stop("comparison data did not produce a 'dfm' object")
  }
  X_old <- dfm_news_restore_missing(object$X_imp)
  X_new <- dfm_news_restore_missing(dfm_new$X_imp)

  stats_old <- attr(object$X_imp, "stats")
  stats_new <- attr(dfm_new$X_imp, "stats")
  if(is.null(stats_new)) stop("dfm_new lacks stats to standardize data")
  if(is.null(stats_old)) stop("dfm object lacks stats to standardize data")

  X_old <- dfm_news_scale(unscale(X_old, stats_old), stats_new)

  if(nrow(X_old) != nrow(X_new) || ncol(X_old) != ncol(X_new)) {
    stop("dfm objects have incompatible data dimensions")
  }
  if(!is.null(colnames(X_old)) && !is.null(colnames(X_new)) && any(colnames(X_old) != colnames(X_new))) {
    stop("dfm objects have different variable ordering")
  }

  if(!is.numeric(t.fcst) || length(t.fcst) != 1L) stop("t.fcst must be a single integer index")
  t_fcst <- as.integer(t.fcst)
  if(t_fcst < 1L || t_fcst > nrow(X_old)) stop("t.fcst is out of bounds")

  vars_idx <- resolve_vars(target.vars, ncol(X_old), colnames(X_old))

  r_old <- nrow(object$A)
  if(ncol(object$A) %% r_old != 0L) stop("Invalid transition matrix dimensions in dfm object")
  p_old <- ncol(object$A) / r_old
  r <- nrow(dfm_new$A)
  if(ncol(dfm_new$A) %% r != 0L) stop("Invalid transition matrix dimensions in dfm_new")
  p <- ncol(dfm_new$A) / r
  if(r_old != r || p_old != p) stop("dfm objects have incompatible r or p")

  n <- ncol(X_old)
  if(is.null(series)) series <- colnames(X_old)
  if(is.null(series)) series <- paste0("Series", seq_len(n))
  if(length(series) != n) stop("series must have length equal to the number of variables")

  stats <- dfm_news_stats(dfm_new$X_imp)
  Mx <- stats$Mx
  Wx <- stats$Wx
  scale_vec <- if(standardized) rep(1, n) else Wx

  # Detect model type: MQ and/or AR1
  quarterly.vars <- dfm_new$quarterly.vars
  idio.ar1 <- !is.null(dfm_new$rho)
  state <- dfm_news_state(dfm_new, require_full = length(quarterly.vars) || idio.ar1)

  # Compute releases ONCE (shared across all targets)
  rel <- which(is.na(X_old) & !is.na(X_new), arr.ind = TRUE)
  n_news <- nrow(rel)

  # Pre-compute KFS results and shared matrices
  if(n_news > 0L) {
    t_miss <- rel[, 1L]
    v_miss <- rel[, 2L]
    lag <- t_fcst - t_miss
    k <- as.integer(max(c(abs(lag), max(lag) - min(lag))))

    Res_old <- dfm_news_kfs(X_old, state, k)
    Res_new <- dfm_news_kfs(X_new, state, 0L)

    # Pre-compute innovations
    innov <- numeric(n_news)
    for(i in seq_len(n_news)) {
      innov[i] <- X_new[t_miss[i], v_miss[i]] - Res_old$X_sm[t_miss[i], v_miss[i]]
    }

    # Pre-compute P1 (r x n_news)
    P <- Res_old$P
    state_dim <- Res_old$state_dim
    C_use <- Res_old$C
    R_use <- Res_old$R
    P1 <- matrix(0, state_dim, n_news)
    for(i in seq_len(n_news)) {
      h <- abs(t_fcst - t_miss[i])
      m <- max(t_miss[i], t_fcst)
      idx <- (h * state_dim + 1L):(h * state_dim + state_dim)
      Pp <- matrix(P[1:state_dim, idx, m], state_dim, state_dim)
      if(t_miss[i] > t_fcst) Pp <- t(Pp)
      P1[, i] <- Pp %*% t(C_use[v_miss[i], , drop = FALSE])
    }

    # Pre-compute P2 with symmetry optimization (upper triangle only)
    P2 <- matrix(0, n_news, n_news)
    for(i in seq_len(n_news)) {
      for(j in i:n_news) {
        h <- abs(lag[i] - lag[j])
        m <- max(t_miss[i], t_miss[j])
        idx <- (h * state_dim + 1L):(h * state_dim + state_dim)
        Pp <- matrix(P[1:state_dim, idx, m], state_dim, state_dim)
        if(t_miss[j] > t_miss[i]) Pp <- t(Pp)
        WW <- if(v_miss[i] == v_miss[j] && t_miss[i] != t_miss[j]) 0 else R_use[v_miss[i], v_miss[j]]
        P2[i, j] <- drop(C_use[v_miss[i], , drop = FALSE] %*% Pp %*% t(C_use[v_miss[j], , drop = FALSE])) + WW
      }
    }
    # Mirror to lower triangle (P2 is symmetric)
    P2[lower.tri(P2)] <- t(P2)[lower.tri(P2)]

    if(qr(P2)$rank < n_news) stop("P2 is singular; cannot compute news weights")

    # Pre-compute P1 %*% solve(P2) for efficiency
    P1_P2inv <- P1 %*% solve(P2)

    # Pre-compute actual and forecasts (shared across targets)
    actual <- forecast <- rep(NA_real_, n)
    for(i in seq_len(n_news)) {
      actual[v_miss[i]] <- X_new[t_miss[i], v_miss[i]]
      forecast[v_miss[i]] <- Res_old$X_sm[t_miss[i], v_miss[i]]
    }

  } else {
    # No new releases
    Res_old <- dfm_news_kfs(X_old, state, 0L)
    Res_new <- dfm_news_kfs(X_new, state, 0L)
  }

  # Inner function: now only computes target-specific outputs
  compute_news <- function(v_news) {
    # Case 1: Target is directly observed in new data
    if(!is.na(X_new[t_fcst, v_news])) {
      y_old <- Res_old$X_sm[t_fcst, v_news]
      y_new <- X_new[t_fcst, v_news]
      if(!standardized) {
        y_old <- y_old * Wx[v_news] + Mx[v_news]
        y_new <- y_new * Wx[v_news] + Mx[v_news]
      }

      na_vec <- rep(NA_real_, n)
      gain <- gain_std <- impact <- news <- numeric(n)
      impact[v_news] <- news[v_news] <- y_new - y_old
      gain[v_news] <- 1
      gain_std[v_news] <- 1

      return(list(y_old = y_old, y_new = y_new,
                  news_df = data.frame(
                    series = series,
                    actual = na_vec,
                    forecast = na_vec,
                    news = news,
                    gain = gain,
                    gain_std = gain_std,
                    impact = impact)))
    }

    # Case 2: No new releases
    if(n_news == 0L) {
      y_old <- Res_old$X_sm[t_fcst, v_news]
      y_new <- Res_new$X_sm[t_fcst, v_news]
      if(!standardized) {
        y_old <- y_old * Wx[v_news] + Mx[v_news]
        y_new <- y_new * Wx[v_news] + Mx[v_news]
      }
      na_vec <- rep(NA_real_, n)
      gain <- gain_std <- impact <- news <- numeric(n)

      return(list(y_old = y_old, y_new = y_new,
                  news_df = data.frame(
                    series = series,
                    actual = na_vec,
                    forecast = na_vec,
                    news = news,
                    gain = gain,
                    gain_std = gain_std,
                    impact = impact)))
    }

    # Case 3: Main news decomposition (uses pre-computed P1_P2inv, innov)
    y_old <- Res_old$X_sm[t_fcst, v_news]
    y_new <- Res_new$X_sm[t_fcst, v_news]
    if(!standardized) {
      y_old <- y_old * Wx[v_news] + Mx[v_news]
      y_new <- y_new * Wx[v_news] + Mx[v_news]
    }

    # Compute gain using pre-computed P1_P2inv (standardized innovations)
    gain_std_release <- drop(scale_vec[v_news] * (C_use[v_news, , drop = FALSE] %*% P1_P2inv))
    temp <- gain_std_release * innov

    actual <- forecast <- rep(NA_real_, n)
    impact <- news <- news_std <- numeric(n)
    X_old_sm <- Res_old$X_sm
    for(i in seq_len(n_news)) {
      v <- v_miss[i]
      actual_i <- X_new[t_miss[i], v]
      forecast_i <- X_old_sm[t_miss[i], v]
      if(!standardized) {
        actual_i <- actual_i * Wx[v] + Mx[v]
        forecast_i <- forecast_i * Wx[v] + Mx[v]
      }
      actual[v] <- actual_i
      forecast[v] <- forecast_i
      news_std[v] <- news_std[v] + innov[i]
      news[v] <- news[v] + (actual_i - forecast_i)
      impact[v] <- impact[v] + temp[i]
    }

    gain_out <- numeric(n)
    nz <- news != 0
    gain_out[nz] <- impact[nz] / news[nz]

    gain_std <- numeric(n)
    nz_std <- news_std != 0
    gain_std[nz_std] <- impact[nz_std] / news_std[nz_std]

    list(y_old = y_old,
         y_new = y_new,
         news_df = data.frame(
           series = series,
           actual = actual,
           forecast = forecast,
           news = news,
           gain = gain_out,
           gain_std = gain_std,
           impact = impact))
  }

  res <- lapply(vars_idx, compute_news)
  if(length(vars_idx) == 1L) {
    res <- res[[1L]]
    res$target.var <- setNames(vars_idx, series[vars_idx])
    res$t.fcst <- t_fcst
    res$standardized <- standardized
    class(res) <- "dfm_news"
    return(res)
  }
  if(!is.null(colnames(X_old))) names(res) <- colnames(X_old)[vars_idx]
  attr(res, "target.vars") <- setNames(vars_idx, series[vars_idx])
  attr(res, "t.fcst") <- t_fcst
  attr(res, "standardized") <- standardized
  class(res) <- "dfm_news_list"
  res
}

#' @rdname news
#' @param x an object of class 'dfm_news' or 'dfm_news_list'.
#' @param digits integer. Number of digits to print.
#' @param \dots not used.
#' @export
print.dfm_news <- function(x, digits = 4L, ...) {
  cat("DFM News\n")
  cat("Target variable:", names(x$target.var), "\n")
  cat("Target time:", x$t.fcst, "\n")
  cat("Old forecast:", round(x$y_old, digits), "\n")
  cat("New forecast:", round(x$y_new, digits), "\n")
  cat("Revision:", round(x$y_new - x$y_old, digits), "\n")
  cat("Standardized:", isTRUE(x$standardized), "\n")
  return(invisible(x))
}

#' @rdname news
#' @export
print.dfm_news_list <- function(x, digits = 4L, ...) {
  t_fcst <- attr(x, "t.fcst")
  standardized <- attr(x, "standardized")
  cat("DFM News (Multiple Targets)\n")
  cat("Target time:", t_fcst, "\n")
  cat("Targets:", length(x), "\n")
  cat("Standardized:", isTRUE(standardized), "\n")
  tbl <- sapply(x, function(res) c(y_old = unattrib(res$y_old), y_new = unattrib(res$y_new), revision = unattrib(res$y_new - res$y_old)))
  if(is.null(dim(tbl))) tbl <- matrix(tbl, nrow = 3L)
  print(round(t(tbl), digits))
  return(invisible(x))
}

#' @rdname news
#' @param name character. Element name.
#' @param i index. Element position or name.
#' @export
`$.dfm_news_list` <- function(x, name) {
  i <- match(name, names(x))
  if(is.na(i)) return(NULL)
  res <- unclass(x)[[i]]
  res$target.var <- attr(x, "target.vars")[i]
  res$t.fcst <- attr(x, "t.fcst")
  res$standardized <- attr(x, "standardized")
  class(res) <- "dfm_news"
  res
}

#' @rdname news
#' @export
`[[.dfm_news_list` <- function(x, i) {
  if(is.character(i)) {
    i <- match(i, names(x))
    if(is.na(i)) return(NULL)
  }
  res <- unclass(x)[[i]]
  res$target.var <- attr(x, "target.vars")[i]
  res$t.fcst <- attr(x, "t.fcst")
  res$standardized <- attr(x, "standardized")
  class(res) <- "dfm_news"
  res
}

#' @rdname news
#' @export
`[.dfm_news_list` <- function(x, i) {
  res <- unclass(x)[i]
  attr(res, "target.vars") <- attr(x, "target.vars")[i]
  attr(res, "t.fcst") <- attr(x, "t.fcst")
  attr(res, "standardized") <- attr(x, "standardized")
  class(res) <- "dfm_news_list"
  res
}

#' @rdname news
#' @importFrom collapse rowbind
#' @export
as.data.frame.dfm_news_list <- function(x, ...) {
  res <- lapply(x, .subset2, "news_df") |>
   rowbind(idcol = "target")
  attr(res, "target.vars") <- attr(x, "target.vars")
  attr(res, "t.fcst") <- attr(x, "t.fcst")
  attr(res, "standardized") <- attr(x, "standardized")
  res
}

#% @aliases forecast.dfm
#' @name predict.dfm
#' @aliases print.dfm_forecast
#' @aliases plot.dfm_forecast
#'
#' @title DFM Forecasts
#'
#' @description This function produces h-step ahead forecasts of both the factors and the data,
#' with an option to also forecast autocorrelated residuals with a univariate method and produce a combined forecast.
#'
#' @param object an object of class 'dfm'.
#' @param h integer. The forecast horizon.
#' @param method character. The factor estimates to use: one of \code{"qml"}, \code{"2s"} or \code{"pca"}.
#' @param standardized logical. \code{FALSE} will return data forecasts on the original scale.
#' @param use.full.state logical. Use the full state-space (if available) when computing residuals for optional residual forecasting. When \code{idio.ar1 = TRUE}, this yields residuals after both factor and idiosyncratic components; set to \code{FALSE} to use factor-only residuals. Falls back to the compact form if unavailable or if \code{method = "pca"}.
#' @param resFUN an (optional) function to compute a univariate forecast of the residuals.
#' The function needs to have a second argument providing the forecast horizon (\code{h}) and return a vector of forecasts. See Examples.
#' @param resAC numeric. Threshold for residual autocorrelation to apply \code{resFUN}: only residual series where AC1 > resAC will be forecasted.
#' @param \dots further arguments to \code{resFUN}.
#'
#' @returns A list-like object of class 'dfm_forecast' with the following elements:
#'  \item{\code{X_fcst}}{\eqn{h \times n}{h x n} matrix with the forecasts of the variables. }
#'  \item{\code{F_fcst}}{\eqn{h \times r}{h x r} matrix with the factor forecasts. }
#'  \item{\code{X}}{\eqn{T \times n}{T x n} matrix with the standardized (scaled and centered) data - with attributes attached allowing reconstruction of the original data:
#'       \tabular{llll}{
#'       \code{"stats"} \tab\tab is a \eqn{n \times 5}{n x 5} matrix of summary statistics of class \code{"qsu"} (see \code{\link[collapse]{qsu}}). Only attached if \code{standardized = TRUE}. \cr\cr
#'       \code{"attributes"} \tab\tab contains the \code{\link{attributes}} of the original data input.\cr\cr
#'       \code{"is.list"} \tab\tab is a logical value indicating whether the original data input was a list / data frame. \cr\cr
#'       }
#'      }
#'  \item{\code{F}}{\eqn{T \times r}{T x r} matrix of factor estimates. }
#'  \item{\code{method}}{the factor estimation method used.}
#'  \item{\code{anyNA}}{logical indicating whether \code{X} contains any missing values.}
#'  \item{\code{h}}{the forecast horizon.}
#'  \item{\code{resid.fc}}{logical indicating whether a univariate forecasting function was applied to the residuals.}
#'  \item{\code{resid.fc.ind}}{indices indicating for which variables (columns of \code{X}) the residuals were forecasted using the univariate function.}
#'  \item{\code{call}}{call object obtained from \code{match.call()}.}
#'
#' @seealso \link{dfms-package}
#'
#' @examples \donttest{
#' library(xts)
#' library(collapse)
#'
#' # Fit DFM with 3 factors and 3 lags in the transition equation
#' mod <- DFM(diff(BM14_M), r = 3, p = 3)
#'
#' # 15 period ahead forecast
#' fc <- predict(mod, h = 15)
#' print(fc)
#' plot(fc, xlim = c(300, 370))
#'
#' # Also forecasting autocorrelated residuals with an AR(1)
#' fcfun <- function(x, h) predict(ar(na_rm(x)), n.ahead = h)$pred
#' fcar <- predict(mod, resFUN = fcfun, h = 15)
#' plot(fcar, xlim = c(300, 370))
#'
#' # Retrieving a data frame of the forecasts
#' head(as.data.frame(fcar, pivot = "wide")) # Factors
#' head(as.data.frame(fcar, use = "data"))   # Data
#' head(as.data.frame(fcar, use = "both"))   # Both
#' }
#' @export
# TODO: Option for prediction in original format??
predict.dfm <- function(object,
                        h = 10L,
                        method = switch(object$em.method, none = "2s", "qml"),
                        standardized = TRUE,
                        use.full.state = TRUE,
                        resFUN = NULL,
                        resAC = 0.1, ...) {
  method <- tolower(method)
  X <- object$X_imp
  Fa <- switch(method,
              pca = object$F_pca, `2s` = object$F_2s, qml = object$F_qml,
              stop("Unkown method", method))

  ss_full <- if(isTRUE(use.full.state) && method != "pca") object$ss_full else NULL
  use_state_forecast <- !is.null(ss_full) &&
    !is.null(ss_full$A) && !is.null(ss_full$C) && !is.null(ss_full$F_smooth)

  if(use_state_forecast) {
    A_full <- ss_full$A
    C_full <- ss_full$C
    F_last <- drop(ftail(ss_full$F_smooth, 1L))
    state_dim <- ncol(ss_full$F_smooth)
    r <- nrow(object$A)
    ny <- nrow(C_full)

    F_fc <- matrix(NA_real_, nrow = h, ncol = r)
    X_fc <- matrix(NA_real_, nrow = h, ncol = ny)
    for(i in seq_len(h)) {
      F_last <- drop(A_full %*% F_last)
      F_fc[i, ] <- F_last[seq_len(r)]
      X_fc[i, ] <- C_full %*% F_last
    }
  } else {
    nf <- dim(Fa)[2L]
    C <- object$C
    ny <- dim(C)[1L]
    A <- object$A
    r <- dim(A)[1L]
    p <- dim(A)[2L] / r

    F_fc <- matrix(NA_real_, nrow = h, ncol = nf)
    X_fc <- matrix(NA_real_, nrow = h, ncol = ny)

    # DFM forecasting loop
    if(is.null(object$quarterly.vars)) {
      F_last <- ftail(Fa, p)   # dimnames(F_last) <- list(c("L2", "L1"), c("f1", "f2"))
      for (i in seq_len(h)) {
        F_reg <- ftail(F_last, p)
        F_fc[i, ] <- tmp <- A %*% vec(t(F_reg)[, p:1, drop = FALSE])
        dim(tmp) <- NULL
        X_fc[i, ] <- C %*% tmp
        F_last <- rbind(F_reg, tmp)
      }
    } else {
      F_last <- ftail(Fa, max(p, 5))
      qind <- ckmatch(object$quarterly.vars, dimnames(X)[[2L]])
      Cq_lags <- object$C[qind, rep(seq_len(ncol(Fa)), each = 5), drop = FALSE] %r*% rep(c(1, 2, 3, 2, 1), ncol(Fa))
      # Mixed frequency forecasting loop
      for (i in seq_len(h)) {
        F_reg <- ftailrev(F_last, p)
        F_fc[i, ] <- tmp <- A %*% vec(t(F_reg))
        dim(tmp) <- NULL
        X_fc[i, -qind] <- C[-qind,, drop = FALSE] %*% tmp
        F_last <- rbind(F_last, tmp)
        X_fc[i, qind] <- Cq_lags %*% vec(ftailrev(F_last, 5))
      }
    }
  }

  # Additional univariate forecasting of the residuals?
  fcr <- NULL
  if(!is.null(resFUN)) {
    if(!is.function(resFUN)) stop("resFUN needs to be a forecasting function with second argument h that produces a numeric h-step ahead forecast of a univariate time series")
    # If X is a multivariate time series object for which the univariate forecasting function could have methods.
    ofl <- !attr(X, "is.list") && length(attr(X, "attributes")[["class"]])
    rsid <- residuals(object, method, orig.format = ofl, standardized = TRUE, na.keep = FALSE,
                      use.full.state = use.full.state) # TODO: What about missing values??
    if(ofl && length(object$rm.rows)) rsid <- rsid[-object$rm.rows, , drop = FALSE]
    ACF <- AC1(rsid, object$anyNA)
    fcr <- which(abs(ACF) >= abs(resAC)) # TODO: Check length of forecast??
    for (i in fcr) X_fc[, i] <- X_fc[, i] + as.numeric(resFUN(rsid[, i], h, ...))
  } else if(!is.null(res <- object[["e"]]) && !use_state_forecast) {
    rho <- object$rho
    last_res <- res[nrow(res), ]
    for (i in seq_len(h)) {
      last_res <- last_res * rho
      X_fc[i, ] <- X_fc[i, ] + last_res
    }
  }

  if(!standardized) { # Unstandardize factors with the average mean and SD??
    stats <- attr(X, "stats")
    attr(X, "stats") <- NULL
    X_fc <- unscale(X_fc, stats)
    X <- unscale(X, stats)
  }

  dimnames(X_fc) <- dimnames(X)
  dimnames(F_fc) <- dimnames(Fa)

  if(object$anyNA) {
    X[attr(X, "missing")] <- NA
    attr(X, "missing") <- NULL
  }

  # model = object, # Better only save essential objects...
  res <- list(X_fcst = X_fc,
              F_fcst = F_fc,
              X = X,
              F = Fa,
              method = method,
              anyNA = object$anyNA,
              h = h,
              resid.fc = !is.null(resFUN), # TODO: Rename list elements??
              resid.fc.ind = fcr,
              call = match.call())

  class(res) <- "dfm_forecast"
  return(res)
}
# forecast.dfm <- predict.dfm

#' @rdname predict.dfm
#' @param x object of type 'dfm_forecast', returned from \code{predict.dfm}.
#' @param digits integer. The number of digits to print out.
#' @param \dots not used.
#' @export
print.dfm_forecast <- function(x,
                               digits = 4L, ...) {
  h <- x$h
  cat(h, "Step Ahead Forecast from Dynamic Factor Model\n\n")
  cat("Factor Forecasts\n")
  F_fcst <- x$F_fcst
  dimnames(F_fcst)[[1L]] <- seq_len(h)
  print(round(F_fcst, digits))
  cat("\nSeries Forecasts\n")
  X_fcst <- x$X_fcst
  dimnames(X_fcst)[[1L]] <- seq_len(h)
  print(round(X_fcst, digits))
  return(invisible(x))
}

#' @rdname predict.dfm
#' @param main,xlab,ylab character. Graphical parameters passed to \code{\link{ts.plot}}.
#' @param factors integers indicating which factors to display. Setting this to \code{NA}, \code{NULL} or \code{0} will omit factor plots.
#' @param scale.factors logical. Standardize factor estimates, this usually improves the plot since the factor estimates corresponding to the greatest PCA eigenvalues tend to have a greater variance than the data.
#' @param factor.col,factor.lwd graphical parameters affecting the colour and line width of factor estimates plots. See \code{\link{par}}.
#' @param fcst.lty integer or character giving the line type of the forecasts of factors and data. See \code{\link{par}}.
#' @param data.col character vector of length 2 indicating the colours of historical data and forecasts of that data. Setting this to \code{NA}, \code{NULL} or \code{""} will not plot data and data forecasts.
#' @param legend logical. \code{TRUE} draws a legend in the top-left of the chart.
#' @param legend.items character names of factors for the legend.
#' @param grid logical. \code{TRUE} draws a grid on the background of the plot.
#' @param vline logical. \code{TRUE} draws a vertical line deliminating historical data and forecasts.
#' @param vline.lty,vline.col graphical parameters affecting the appearance of the vertical line. See \code{\link{par}}.
#' @param \dots further arguments passed to \code{\link{ts.plot}}. Sensible choices are \code{xlim} and \code{ylim} to restrict the plot range.
#' @importFrom collapse setop
#' @export
# TODO: multiple plot types...# , type = c("joint", "individual")
# also arguments show = c("both", "factors", "data"), and
# Also put plot on original timescale if ts object
# TODO: Option to unstandardize factors.
plot.dfm_forecast <- function(x,
                              main = paste(x$h, "Period Ahead DFM Forecast"),
                              xlab = "Time", ylab = "Standardized Data",
                              factors = seq_len(ncol(x$F)),
                              scale.factors = TRUE,
                              factor.col = rainbow(length(factors)),
                              factor.lwd = 1.5,
                              fcst.lty = "dashed",
                              data.col = c("grey85", "grey65"),
                              legend = TRUE,
                              legend.items = paste0("f", factors),
                              grid = FALSE, vline = TRUE,
                              vline.lty = "dotted", vline.col = "black", ...) {

  dcl <- is.character(data.col[1L]) && nzchar(data.col[1L])
  ffl <- length(factors) && !is.na(factors[1L]) && factors[1L] > 0L
  nyliml <- !(...length() && any(names(list(...)) == "ylim")) # ...names() -> Added after R 3.3.0
  if(!ffl) factors <- 1L
  Fa <- x$F[, factors, drop = FALSE]
  r <- ncol(Fa)
  TT <- nrow(Fa)
  if(ffl) {
    F_fcst <- x$F_fcst[, factors, drop = FALSE]
    if(scale.factors) {
      fcstat <- qsu(Fa)
      F_fcst <- setop(TRA.matrix(F_fcst, fcstat[, "Mean"], "-"), "/", fcstat[, "SD"], rowwise = TRUE) # Unscale ??
      Fa <- fscale(Fa)
    }
    if(nyliml) Fr <- frange(Fa)
    Fa <- rbind(Fa, matrix(NA_real_, x$h, r))
  } else Fr <- NULL
  if(dcl) {
    X <- x$X
    n <- ncol(X)
    if(nyliml) {
      Xr <- frange(X, na.rm = TRUE)
      Pr <- frange(if(ffl) c(F_fcst, x$X_fcst) else x$X_fcst)
    }
    X_fcst <- rbind(matrix(NA_real_, TT-1L, n), X[TT, , drop = FALSE], x$X_fcst)
    X <- rbind(X, matrix(NA_real_, x$h, n))
  } else {
    data.col <- Xr <- NULL
    X <- Fa[, 1L]
    if(nyliml) Pr <- frange(F_fcst)
  }
  if(ffl) F_fcst <- rbind(matrix(NA_real_, TT-1L, r), Fa[TT, , drop = FALSE], F_fcst)
  if(nyliml) {
    ts.plot(X, col = data.col[1L],
            ylim = c(min(Xr[1L], Fr[1L], Pr[1L]), max(Xr[2L], Fr[2L], Pr[1L])),
            main = main, xlab = xlab, ylab = ylab, ...)
  } else ts.plot(X, col = data.col[1L], main = main, xlab = xlab, ylab = ylab, ...)
  if(grid) grid()
  if(dcl) for (i in seq_len(n)) lines(X_fcst[, i], col = data.col[2L], lty = fcst.lty)
  if(ffl) for (i in seq_len(r)) {
    lines(Fa[, i], col = factor.col[i], lwd = factor.lwd)
    lines(F_fcst[, i], col = factor.col[i], lwd = factor.lwd, lty = fcst.lty)
  }
  if(ffl && legend) legend("topleft", legend.items, col = factor.col,
                           lwd = factor.lwd, lty = 1L, bty = "n")
  if(vline) abline(v = TT, col = vline.col, lwd = 1L, lty = vline.lty)
}


#' @rdname predict.dfm
#' @param x an object class 'dfm_forecast'.
#' @param use character. Which forecasts to use \code{"factors"}, \code{"data"} or \code{"both"}.
#' @param pivot character. The orientation of the frame: \code{"long"} or \code{"wide"}.
#' @param time a vector identifying the time dimension, must be of length T + h, or \code{NULL} to omit a time variable.
#' @param stringsAsFactors logical. If \code{TRUE} and \code{pivot = "long"} the 'Variable' column is created as a factor. Same as option to \code{\link{as.data.frame.table}}.
#' @param \dots not used.
#'
#' @export
as.data.frame.dfm_forecast <- function(x, ...,
                              use = c("factors", "data", "both"),
                              pivot = c("long", "wide"),
                              time = seq_len(nrow(x$F) + x$h),
                              stringsAsFactors = TRUE) {

  mat <- switch(tolower(use[1L]),
                factors = rbind(x$F, x$F_fcst),
                data = rbind(x$X, x$X_fcst),
                both = cbind(rbind(x$F, x$F_fcst), rbind(x$X, x$X_fcst)),
                stop("Unknown use option:", use[1L]))

  fcvec <- c(rep(FALSE, nrow(x$F)), rep(TRUE, x$h))
  TT <- nrow(mat)
  r <- ncol(mat)
  if(!is.null(time) && length(time) != TT) stop(sprintf("time must be a length %s vector or NULL", TT))

  res <- switch(tolower(pivot[1L]),
      long = list(Variable = if(stringsAsFactors) setAttrib(rep(1:r, each = TT), list(levels = dimnames(mat)[[2L]], class = "factor")) else
                                                rep(dimnames(mat)[[2L]], each = TT),
                  Time = if(length(time)) rep(time, r) else NULL,
                  Forecast = rep(fcvec, r),
                  Value = unattrib(mat)),
      wide = c(list(Time = time, Forecast = fcvec), mctl(mat, TRUE)),
      stop("Unknown pivot option:", pivot[1L])
  )

  if(is.null(time)) res <- na_rm(res)
  attr(res, "row.names") <- .set_row_names(length(res[[1L]]))
  class(res) <- "data.frame"
  return(res)
}

# interpolate.dfm <- function(x, method = "qml", interpolate = TRUE) {
#   W <- is.na(data)
#   stats <- qsu(data)
#   STDdata <- fscale(data)
#   if(nrow(x$C) != ncol(data)) stop("dimension mismatch")
#   Fcst <- tcrossprod(x[[method]], x$C)
#   # TODO: Make this work for data.table...
#   STDdata[W] <- Fcst[W]
#   STDdata <- ((STDdata %r*% stats[, "SD"]) %r+% stats[, "Mean"])
#   data[W] <- STDdata[W]
#   data
# }
#
# nowcast.dfm <- function(x, method = "qml", ...) {
# }
#
# backcast.dfm <- function(x, method = "qml", ...) {
# }

# Adapted from: https://github.com/nmecsys/nowcasting/blob/master/R/ICfactors.R
#' @title Information Criteria to Determine the Number of Factors (r)
#' @description Minimizes 3 information criteria proposed by Bai and Ng (2002) to determine the optimal number of factors r* to be used in an approximate factor model.
#' A Screeplot can also be computed to eyeball the number of factors in the spirit of Onatski (2010).
#' @param X a \code{T x n} numeric data matrix or frame of stationary time series.
#' @param max.r integer. The maximum number of factors for which IC should be computed (or eigenvalues to be displayed in the screeplot).
#' @param x an object of type 'ICr'.
#' @param \dots further arguments to \code{\link{ts.plot}} or \code{\link{plot}}.
#'
#' @return A list of 4 elements:
#' \item{F_pca}{\code{T x n} matrix of principle component factor estimates.}
#' \item{eigenvalues}{the eigenvalues of the covariance matrix of \code{X}.}
#' \item{IC}{\code{r.max x 3} 'table' containing the 3 information criteria of Bai and Ng (2002), computed for all values of \code{r} from \code{1:r.max}.}
#' \item{r.star}{vector of length 3 containing the number of factors (\code{r}) minimizing each information criterion.}
#'
#' @details Following Bai and Ng (2002) and De Valk et al. (2019), let \eqn{NSSR(r)}{NSSR(r)} be the normalized sum of squared residuals \eqn{SSR(r) / (n \times T)}{SSR(r) / (n x T)} when r factors are estimated using principal components.
#' Then the information criteria can be written as follows:
#'
#' \deqn{IC_{r1} = \ln(NSSR(r)) + r\left(\frac{n + T}{nT}\right) + \ln\left(\frac{nT}{n + T}\right)}{ICr1 = ln(NSSR(r)) + r * (n + T)/(n * T) + ln((n * T)/(n + T))}
#' \deqn{IC_{r2} = \ln(NSSR(r)) + r\left(\frac{n + T}{nT}\right) + \ln(\min(n, T))}{ICr2 = ln(NSSR(r)) + r * (n + T)/(n * T) + ln(min(n, T))}
#' \deqn{IC_{r3} = \ln(NSSR(r)) + r\left(\frac{\ln(\min(n, T))}{\min(n, T)}\right)}{ICr3 = ln(NSSR(r)) + r * ln(min(n, T))/min(n, T)}
#'
#' The optimal number of factors r* corresponds to the minimum IC. The three criteria are are asymptotically equivalent, but may give significantly
#' different results for finite samples. The penalty in \eqn{IC_{r2}}{ICr2} is highest in finite samples.
#'
#' In the Screeplot a horizontal dashed line is shown signifying an eigenvalue of 1, or a share of variance corresponding to 1 divided by the number of eigenvalues.
#'
#' @note To determine the number of lags (\code{p}) in the factor transition equation, use the function \code{vars::VARselect} with r* principle components (also returned by \code{ICr}).
#'
#' @seealso \link{dfms-package}
#'
#' @examples
#' library(xts)
#' library(vars)
#'
#' ics <- ICr(diff(BM14_M))
#' print(ics)
#' plot(ics)
#' screeplot(ics)
#'
#' # Optimal lag-order with 6 factors chosen
#' VARselect(ics$F_pca[, 1:6])
#'
#' @references
#' Bai, J., Ng, S. (2002). Determining the Number of Factors in Approximate Factor Models. \emph{Econometrica, 70}(1), 191-221. https://doi.org/10.1111/1468-0262.00273.
#'
#' Onatski, A. (2010). Determining the number of factors from empirical distribution of eigenvalues. \emph{The Review of Economics and Statistics, 92}(4), 1004-1016.
#'
#' De Valk, S., de Mattos, D., & Ferreira, P. (2019). Nowcasting: An R package for predicting economic variables using dynamic factor models. \emph{The R Journal, 11}(1), 230-244.
#' @export
ICr <- function(X, max.r = min(20, ncol(X)-1)) {

  # Converting to matrix and standardizing
  X <- fscale(qM(X), na.rm = TRUE)
  dimnames(X) <- NULL

  if(anyNA(X)) {
    message("Missing values detected: imputing data with tsnarmimp() with default settings")
    X <- tsnarmimp(X)
    attributes(X) <- list(dim = dim(X))
  }

  n <- ncol(X)
  TT <- nrow(X)

  # defining rmax and checking if it is a positive integer
  if(!is.numeric(max.r) || max.r < 1) stop("max.r needs to be a positive integer")
  max.r <- if(max.r > n) n else as.integer(max.r)

  # Eigen decomposition
  eigen_decomp <- eigen(cov(X), symmetric = TRUE)
  evs <- eigen_decomp$vectors
  F_pca <- X %*% evs

  # Various constant terms, according to the 3 criteria of Bai and Ng (2002)
  Tn <- TT * n
  npTdTn <- (n + TT) / Tn
  minnT <- min(n, TT)
  c1 <- npTdTn * log(1/npTdTn)
  c2 <- npTdTn * log(minnT)
  c3 <- log(minnT) / minnT
  cvec <- c(c1, c2, c3)
  result <- matrix(0, max.r, 3)

  # Calculating the IC
  for (r in 1:max.r) {
    # Residuals from r PC's
    res <- X - tcrossprod(F_pca[, 1:r], evs[, 1:r])
    # Log normalized sum of squared errors
    logV <- log(sum(colSums(res^2)/Tn))
    # Computing criteria
    result[r, ] <- logV + r * cvec
  }

  dimnames(result) <- list(r = 1:max.r, IC = paste0("IC", 1:3))
  class(result) <- "table"
  colnames(F_pca) <- paste0("PC", 1:n)

  res_obj <- list(F_pca = F_pca, eigenvalues = eigen_decomp$values, IC = result, r.star = apply(result, 2, which.min))
  class(res_obj) <- "ICr"
  return(res_obj)
}

#' @rdname ICr
#' @export
print.ICr <- function(x, ...) {
  cat("Optimal Number of Factors (r) from Bai and Ng (2002) Criteria\n\n")
  print(x$r.star)
  return(invisible(x))
}

#' @rdname ICr
#' @importFrom collapse fmin.matrix
#' @importFrom graphics grid points
#' @export
plot.ICr <- function(x, ...) {

  ts.plot(x$IC, gpars = list(xlab = "Number of Factors (r)", ylab = "IC Value", lty = c(2L, 1L, 3L),
                             main = "Optimal Number of Factors (r) from Bai and Ng (2002) Criteria"), ...)
  # grid()
  legend("topleft", paste0(names(x$r.star), ", r* = ", x$r.star), lty = c(2L, 1L, 3L)) # , bty = "n"
  points(x = x$r.star, y = fmin.matrix(x$IC), pch = 19, col ="red")

}

#' @rdname ICr
#' @param type character. Either \code{"ev"} (eigenvalues), \code{"pve"} (percent variance explained), or \code{"cum.pve"} (cumulative PVE). Multiple plots can be requested.
#' @param show.grid logical. \code{TRUE} shows gridlines in each plot.
#' @importFrom stats screeplot
#' @export
screeplot.ICr <- function(x, type = "pve", show.grid = TRUE, max.r = 30, ...) {
  ev <- x$eigenvalues
  n <- length(ev)
  pve <- (ev / sum(ev)) * 100
  cs_pve <- cumsum(pve)

  if(length(ev) > max.r) {
    ev <- ev[1:max.r]
    pve <- pve[1:max.r]
    cs_pve <- cs_pve[1:max.r]
  }

  ## This is smarter, but less flexible...
  # if(length(ev) > 20) {
  #   if(cs_pve[21] > 90) {
  #     pve = pve[1:20]
  #     cs_pve = cs_pve[1:20]
  #   } else {
  #     pve = pve[cs_pve < 90]
  #     cs_pve = cs_pve[cs_pve < 90]
  #   }
  # }

  if(length(type) > 1L) {
    oldpar <- par(mfrow = c(1, length(type)))
    on.exit(par(oldpar))
  }
  if(any(type == "ev")) {
    plot(ev, type = "o", ylab = "Eigenvalue", xlab = "Principal Component", col = "dodgerblue4", ...)
    if(show.grid) grid()
    abline(h = 1, lty = 2)
  }
  if(any(type == "pve")) {
    plot(pve, type = "o", ylab = "% Variance Explained", xlab = "Principal Component", col = "dodgerblue4", ...)
    if(show.grid) grid()
    abline(h = 100 / n, lty = 2)
  }
  if(any(type == "cum.pve")) {
    plot(cs_pve, type = "o", ylab = "Cumulative % Variance Explained", xlab = "Number of Principal Components", col = "brown3", ...)
    if(show.grid) grid()
  }
}

#' @rdname plot.dfm
#' @export
screeplot.dfm <- function(x, ...) {
  xl <- list(eigenvalues = x$eigen$values)
  screeplot.ICr(xl, ...)
}