R/fts.criterion.R

In ottosven/dffm: Approximate Functional Factor Model

#' Information Criteria to Estimate the Number of Factors for Functional Factor Models
#'
#' `fts.criterion` estimates the optimal number of factors (K) and lags (p) for functional factor models based on information criteria. It applies the methodology proposed by Otto and Salish (2024)
#' 
#' Estimates the number of factors K and lags p for the approximate functional factor model using the
#' information criteria proposed in Otto and Salish (2024).
#'
#' @param fdaobj An object of class 'fdaobj', typically the output from `fda.preprocess` or `fts.cumAC`.
#' @param K.max An optional integer specifying the maximum number of factors to consider. If NULL (default), `K.max` will be set to the highest number supported by the given 'fdaobj'.
#' @param p.max An optional integer specifying the maximum number of lags to consider. If NULL (default), `p.max` will be determined based on Schwert's rule of thumb.
#'
#' @return
#' Returns a list containing:
#' \item{IC.min}{A data.frame with the estimated number of factors (K) and lags (p) based on the Bayesian Information Criterion (BIC), Hannan-Quinn Criterion (HQC), and the final prediction error (fFPE) criterion of Aue et al. (2015).}
#' \item{MSE.matrix}{A matrix of Mean Squared Error (MSE) values for all combinations of considered factors and lags.}
#'
#' @export
#' @references
#' * Otto, S., & Salish, N. (2024). Approximate Factor Models For Functional Time Series. arXiv:2201.02532.
#' * Aue, A., Norinho, D. D., & Hörmann, S. (2015). On the prediction of stationary functional time series. Journal of the American Statistical Association, 110(509), 378-392.
#'
#' @examples
#' fed = load.fed()
#' fpcaobj = fda.preprocess(data = fed)
#' criterion_output = fts.criterion(fpcaobj)
#' print(criterion_output$IC.min)
fts.criterion = function(fdaobj, K.max = NULL, p.max = NULL){
  ## check input
  if(!(class(fdaobj) == "fdaobj")){
    stop("Please insert an object of class fdaobj, typically generated by the function fda.preprocess or fts.cumAC.")
  }
  if((fdaobj$operator == "cumAC_operator")){
    FPCA = fda.preprocess(fdaobj$raw.data, observationgrid = fdaobj$observationgrid, workinggrid = fdaobj$workinggrid)
  } else {
    FPCA = fdaobj
    cumAC = fts.cumAC(fdaobj, q0=1)
  }
  ## 
  n = dim(cumAC$scores)[1]
  if(is.null(p.max)) (p.max = max(floor(8*(n/100)^(1/4)),1)) #Schwert rule of thumb
  if(is.null(K.max)) (K.max = min(dim(cumAC$scores)[2],floor(n/(2*min(4,p.max)))))
  ## Evaluate IC
  eval.criterion = function(J,m){
    if(J*m >= (n-J*m)) (return(rep(Inf,4)))
    mse = fts.VARforecast(cumAC, K=J, p=m, h=0)$MSE
    ## Aue et al. (2015) criterion
    factors = FPCA$scores[,1:J, drop=FALSE]
    VARfit = lm(embed(factors, m+1)[,1:J] ~ -1 + embed(factors, m+1)[,-(1:J)])
    if(VARfit$rank != J*m) (return(rep(Inf,4)))
    tr.sigmahat = sum((residuals(VARfit))^2)/(n-(J*m))
    out = c(
      mse,
      log(mse) + J*m*log(n)/n,
      log(mse) + 2*J*m*log(log(n))/n,
      (n+m*J)/n * tr.sigmahat + sum(FPCA$eigenvalues[-c(1:J)])
    )
    return(out)
  }
  IC=list()
  IC.min=matrix(nrow = 3, ncol=2, dimnames = list(c("BIC", "HQC", "fFPE"), c("K.opt", "p.opt")))
  for(j in 1:4) IC[[j]] = matrix(nrow = K.max, ncol = p.max, dimnames = list(paste("K.", 1:K.max, sep = ""), paste("p.", 1:p.max, sep = "")))
  for (i in 1:K.max){
    aux = sapply(1:p.max, eval.criterion, J = i)
    for(j in 1:4) IC[[j]][i,] = aux[j,]
  }
  for(j in 1:3) IC.min[j,] = which(IC[[j+1]] == min(na.omit(IC[[j+1]])), arr.ind = TRUE)[1,]
  MSE.matrix = IC[[1]]
  return(list("IC.min" = IC.min,"MSE.matrix" = MSE.matrix))
}