IDSIA / bayesRecon / #22

#' Estimate covariance from naive or seasonal-naive residuals

#'

#' Estimates via shrinkage the covariance matrix of the residuals of the naive 

#' or seasonal naive forecasts.

#' If the frequency of the time series is > 1, the function chooses between the

#' two methods series-by-series according to a selection criterion.

#' If the frequency is 1 or is not provided, only the naive residuals are used.

#'

#' @param y_train Multivariate time series object or numeric matrix of historical 

#'   observations, with dimensions `T x n` (rows are time points, columns are series).

#' @param freq Positive integer seasonal frequency (optional).

#'   If not provided and `y_train` is a multivariate time series, the frequency 

#'   of the data is used. 

#' @param criterion Character string used when `freq > 1` to choose residuals.

#'   Supported values are `"RSS"` (default) and `"seas-test"`.

#'   `"RSS"` chooses the method with the lower residual sum of squares (RSS), 

#'   while `"seas-test"` uses a statistical test for seasonality (requires `forecast` package).

#'

#' @return A numeric `n x n` shrinkage covariance matrix estimated with

#'   [schaferStrimmer_cov()].

#'

#' @seealso [schaferStrimmer_cov()], [reconc_t()]

#'

#' @keywords internal

#' @export

.compute_naive_cov = function(y_train, freq = NULL, criterion = "RSS") {

  # Compute residuals of naive

  # If freq is not provided and the data is a multivariate time series, 

  # set freq to the frequency of the data

    }

  }

  } else {

    # compute residuals of seasonal naive

    }

    # Choose which residuals to use based on the criterion

      # check if forecast package is installed

        )

      }

      # Seasonality test for each time series

      # Compute RSS for both methods, for each series

    } else {

    }

    # No series is seasonal

      # All series are seasonal

      # Some series are seasonal, some are not (use criterion)

    } else {

    }

  }

}

#' Optimize degrees of freedom (nu) via LOO Cross-Validation

#'

#' @param res Matrix of residuals (n_obs x n_var).

#' @param prior_mean The prior mean covariance matrix (n_var x n_var).

#' @param trim Fraction of observations to trim (0 to 1). Default is 0.1 (10%).

#'

#' @details

#' TODO: check these details

#'

#' \strong{Leave-One-Out (LOO) Cross-Validation:}

#' This function estimates the optimal degrees of freedom \eqn{\nu} by maximizing the

#' out-of-sample predictive performance. This is achieved by computing the

#' log-density of each held-out observation \eqn{\mathbf{r}_i} given the remaining

#' data \eqn{\mathbf{R}_{-i}}. The total objective function is the sum of these

#' predictive log-densities:

#' \deqn{\mathcal{L}(\nu) = \sum_{i=1}^T \log f(\mathbf{r}_i | \mathbf{R}_{-i}, \nu)}

#'

#' \strong{The Log-Density Function:}

#' For each LOO step, the residuals are assumed to follow a

#' multivariate t-distribution. The density is expressed directly as a function of the

#' posterior sum-of-squares matrix \eqn{\Psi}, where \eqn{\Psi} scales implicitly with \eqn{\nu}:

#' \deqn{f(\mathbf{r}_i | \Psi, \nu) = \frac{\Gamma(\frac{\nu + T}{2})}{\Gamma(\frac{\nu + T - p}{2}) \pi^{p/2}} |\Psi|^{-1/2} \left( 1 + \mathbf{r}_i^\top \Psi^{-1} \mathbf{r}_i \right)^{-\frac{\nu + T}{2}}}

#' In the code, \eqn{\Psi} is constructed as:

#' \deqn{\Psi = (\nu - p - 1)\bar{\Sigma}_{prior} + \mathbf{R}^\top\mathbf{R}}

#' By using this formulation, the standard scaling factors \eqn{1/\nu} and \eqn{\nu^{-p/2}}

#' are absorbed into the matrix inverse and determinant, respectively.

#'

#' \strong{Efficient Computation via Sherman-Morrison:}

#' Rather than recomputing \eqn{\Psi_{-i}} and its inverse \eqn{T} times, the function uses

#' the full-sample matrix \eqn{\Psi} and adjusts it using the leverage

#' \eqn{h_i = \mathbf{r}_i^\top \Psi^{-1} \mathbf{r}_i}.

#'

#' Through the Matrix Determinant Lemma and the Sherman-Morrison formula, the

#' internal term \eqn{(1 + \mathbf{r}_i^\top \Psi_{-i}^{-1} \mathbf{r}_i)} simplifies to \eqn{(1 - h_i)^{-1}}.

#' The final log-density contribution used in the code is:

#' \deqn{\log f_i \propto \textrm{const} - \frac{1}{2}\log|\Psi| + \frac{\nu + T - 1}{2} \log(1 - h_i)}

#'

#' @return A list containing the optimization results:

#' * `optimal_nu`: The optimal degrees of freedom found.

#' * `min_neg_log_score`: The minimum negative log score achieved.

#' * `convergence`: The convergence status of the optimizer.

#' * `time_elapsed`: The time taken for the optimization.

#'

#' @import nloptr

#'

#' @export

multi_log_score_optimization <- function(res, prior_mean, trim = 0.1) {

  # Pre-compute cross-product of the residuals (to avoid recomputing inside the loop)

  # Define the Objective Function (Negative Log Score to minimize)

    # Compute posterior Psi

    # Compute Cholesky decomposition

    })

    # Invert (using Cholesky)

    # Calculate log-determinant

    # Log Gamma terms and determinant

    # Compute LOO Terms (using Sherman-Morrison)

    # Combine to get vector of LOO Log-Likelihoods per observation

    # Calculate how many observations to remove (based on 'trim')

    # Sort likelihoods (ascending) and remove the lowest 'n_trim' values

    # Return negative sum because optimizer minimizes

  }

  # Optimization Setup

  )

  # Run Optimization

  )

  ))

}

#' t-Rec: reconciliation via conditioning with uncertain covariance via multivariate t-distribution

#'

#' Reconciles base forecasts in a hierarchy by conditioning on the hierarchical

#' constraints, specified by the aggregation matrix A.

#' The base forecasts are assumed to be jointly Gaussian, conditionally on the

#' covariance matrix of the forecast errors. 

#' A Bayesian approach is adopted to account for the uncertainty of the covariance matrix. 

#' An Inverse-Wishart prior is specified on the covariance matrix, 

#' leading to a multivariate t-distribution for the base forecasts.

#' The reconciliation via conditioning is in closed-form, yielding a multivariate t 

#' reconciled distribution.

#'

#' @param A Matrix (n_upp x n_bott) defining the hierarchy (u = Ab).

#' @param base_fc_mean Vector of base forecasts (length n = n_upp + n_bott).

#' @param y_train mts (or matrix) of historical training data (T x n) used for setting prior parameters.

#' @param residuals Matrix (T x n) of base forecast residuals.

#' @param ... Additional arguments for advanced usage: see details.

#' @param return_upper Logical; if TRUE, also returns parameters for the upper-level reconciled distribution.

#' @param return_parameters Logical; if TRUE, also returns prior and posterior parameters.

#'

#' @details

#' \strong{Standard usage.}

#' 

#' The standard workflow for this function is to provide the in-sample \code{residuals}

#' and the historical training data \code{y_train}.

#' The parameters of the Inverse-Wishart prior distribution of the covariance matrix 

#' are set as follows: 

#' \itemize{

#'   \item Prior scale matrix (Psi): set as the covariance of the residuals of naive (or seasonal naive,

#'         a criterion is used to choose between the two) forecasts computed on \code{y_train}.

#'   \item Prior degrees of freedom (nu): estimated via Bayesian Leave-One-Out

#'         Cross-Validation (LOOCV) to maximize out-of-sample performance.

#' }

#' The posterior distribution of the covariance matrix is still Inverse-Wishart.

#' The parameters of the posterior are computed in closed-form using the sample 

#' covariance of the provided \code{residuals}.

#'

#' \strong{Advanced Options.}

#' 

#' Users can bypass the automated estimation by specifying:

#' \enumerate{

#'   \item \code{prior}: a list with entries 'nu' and 'Psi'.

#'         This skips the LOOCV step for \eqn{\nu_0} and the covariance estimation from \code{y_train}.

#'         It requires \code{residuals} to compute the posterior.

#'   \item \code{posterior}: a list with entries 'nu' and 'Psi'.

#'         This skips all internal estimation and updating logic.

#' }

#' Moreover, users can specify:

#' \itemize{

#'   \item \code{freq}: positive integer, used as frequency of data for the seasonal naive forecast in the specification of the prior scale matrix.

#'         By default, if \code{y_train} is a multivariate time series, the frequency of the data is used; otherwise, it is set to 1 (no seasonality).

#'   \item \code{criterion}: either 'RSS' (default) or 'seas-test', specifying which criterior is used to choose between 

#'                           the naive and seasonal naive forecasts for the specification of the prior scale matrix. 

#'                           'RSS' computes the residual sum of squares for both methods and chooses the one with lower RSS, 

#'                           while 'seas-test' uses a statistical test for seasonality 

#'                           (currently implemented using the number of seasonal differences suggested by the `forecast` package, 

#'                           which must be installed). 

#' }

#' 

#' \strong{Reconciled distribution.}

#' 

#' The reconciled distribution is a multivariate t-distribution, 

#' specified by a vector of means, a scale matrix, and a number of degrees of freedom.

#' These parameters are computed in closed-form.

#' By default, only the parameters of the reconciled distribution for the bottom-level 

#' series are returned. See examples.

#'

#'

#' @references

#' Carrara, C., Corani, G., Azzimonti, D., & Zambon, L. (2025). Modeling the uncertainty on the covariance

#' matrix for probabilistic forecast reconciliation. arXiv preprint arXiv:2506.19554.

#' \url{https://arxiv.org/abs/2506.19554}

#'

#' @return A list containing:

#' \itemize{

#'   \item \code{bottom_rec_mean}: reconciled bottom-level mean.

#'   \item \code{bottom_rec_scale_matrix}: reconciled bottom-level scale matrix.

#'   \item \code{bottom_rec_df}: reconciled degrees of freedom.

#' }

#' If \code{return_upper} is TRUE, also returns:

#' \itemize{

#'  \item \code{upper_rec_mean}: reconciled upper-level mean.

#'  \item \code{upper_rec_scale_matrix}: reconciled upper-level scale matrix.

#'  \item \code{upper_rec_df}: reconciled upper-level degrees of freedom.

#'  }

#' If \code{return_parameters} is TRUE, also returns:

#' \itemize{

#'   \item \code{prior_nu}: prior degrees of freedom.

#'   \item \code{prior_Psi}: prior scale matrix.

#'   \item \code{posterior_nu}: posterior degrees of freedom.

#'   \item \code{posterior_Psi}: posterior scale matrix.

#' }

#'

#' @examples

#' 

#' \donttest{

#' library(bayesRecon)

#'

#' if (requireNamespace("forecast", quietly = TRUE)) {

#' 

#'   set.seed(1234)

#'   n_obs <- 100

#' 

#'   # Simulate 2 bottom series from AR(1) processes

#'   y1 <- arima.sim(model = list(ar = 0.8), n = n_obs)

#'   y2 <- arima.sim(model = list(ar = 0.5), n = n_obs)

#'

#'   y_upper <- y1 + y2   # upper series is the sum of the two bottoms

#'   A <- matrix(c(1, 1), nrow = 1)  # Aggregation matrix

#'

#'   # Fit additive ETS models

#'   fit1 <- forecast::ets(y1, additive.only = TRUE)

#'   fit2 <- forecast::ets(y2, additive.only = TRUE)

#'   fit_upper <- forecast::ets(y_upper, additive.only = TRUE)

#'

#'   # Point forecasts (h = 1)

#'   fc_upper <- as.numeric(forecast::forecast(fit_upper, h = 1)$mean)

#'   fc1 <- as.numeric(forecast::forecast(fit1, h = 1)$mean)

#'   fc2 <- as.numeric(forecast::forecast(fit2, h = 1)$mean)

#'   base_fc_mean <- c(fc_upper, fc1, fc2)

#'

#'   # Residuals and training data (n_obs x n matrices, columns in same order as base_fc_mean)

#'   res <- cbind(residuals(fit_upper), residuals(fit1), residuals(fit2))

#'   y_train <- cbind(y_upper, y1, y2)

#'

#'   # --- 1) Generate joint reconciled samples ---

#'   result <- reconc_t(A, base_fc_mean, y_train = y_train, residuals = res)

#'

#'   # Sample from the reconciled bottom-level t-distribution

#'   n_samples <- 2000

#'   L_chol <- t(chol(result$bottom_rec_scale_matrix))

#'   z <- matrix(rt(ncol(A) * n_samples, df = result$bottom_rec_df), nrow = ncol(A))

#'   bottom_samples <- result$bottom_rec_mean + L_chol %*% z  # 2 x n_samples

#'

#'   # Aggregate bottom samples to get upper samples

#'   upper_samples <- A %*% bottom_samples              

#'   joint_samples <- rbind(upper_samples, bottom_samples)

#'   rownames(joint_samples) <- c("upper", "bottom_1", "bottom_2")

#'

#'   cat("Reconciled means (from samples):\n")

#'   print(round(rowMeans(joint_samples), 3))

#' 

#'   cat("Reconciled standard deviations (from samples):\n")

#'   print(round(apply(joint_samples, 1, sd), 3))

#'

#'   # --- 2) 95% prediction intervals via t-distribution quantiles ---

#'   result2 <- reconc_t(A, base_fc_mean, y_train = y_train,

#'                       residuals = res, return_upper = TRUE)

#'

#'   alpha <- 0.05

#'   # Bottom series intervals

#'   for (i in seq_len(ncol(A))) {

#'     s_i <- sqrt(result2$bottom_rec_scale_matrix[i, i])

#'     lo  <- result2$bottom_rec_mean[i] + s_i * qt(alpha / 2,     df = result2$bottom_rec_df)

#'     hi  <- result2$bottom_rec_mean[i] + s_i * qt(1 - alpha / 2, df = result2$bottom_rec_df)

#'     cat(sprintf("Bottom %d: 95%% PI = [%.3f, %.3f]\n", i, lo, hi))

#'   }

#'   # Upper series interval

#'   s_u <- sqrt(result2$upper_rec_scale_matrix[1, 1])

#'   lo  <- result2$upper_rec_mean[1] + s_u * qt(alpha / 2,     df = result2$upper_rec_df)

#'   hi  <- result2$upper_rec_mean[1] + s_u * qt(1 - alpha / 2, df = result2$upper_rec_df)

#'   cat(sprintf("Upper:    95%% PI = [%.3f, %.3f]\n", lo, hi))

#' }

#' }

#' 

#' @seealso 

#' [reconc_gaussian()]

#'

#' @export

reconc_t <- function(A,

                     base_fc_mean,

                     y_train = NULL,

                     residuals = NULL,

                     ...,

                     return_upper = FALSE,

                     return_parameters = FALSE) {

  }

  ##############################################################################

  ### CASE 1 ###

  # If posterior is provided, check if is a list with entries nu and Psi and extract values

    ### CASE 2 ###

    # If posterior not provided, first check that residuals are provided

  } else {

    # Compute sample covariance of the residuals

    # Shrink to diagonal for numerical reasons

    ### CASE 2a ###

    # If prior is provided, check if is a list with entries nu and Psi and extract values

      ### CASE 2b ###

      # If prior not provided:

      # - compute Psi using the (shrinked) covariance matrix of the residuals of the naive

      #   or seasonal naive forecasts

      # - set nu using LOOCV

    } else {

      # Compute the covariance residuals of the naive or seasonal naive forecasts on training data

      # Use it to set the prior

    }

    # Compute posterior parameters

  }

  ##############################################################################

  # Reconcile via conditioning the t-distribution in closed form

    } else {

    }

  }

}

#' Core reconciliation via multivariate t-distribution.

#'

#' Internal function that performs the core reconciliation logic for the t-distribution based

#' reconciliation method. This function assumes an uncertain covariance matrix with an Inverse-Wishart prior.

#'

#' @param A Matrix (n_upper x n_bottom) defining the hierarchy where upper = A %*% bottom.

#' @param base_fc_mean Vector of length (n_upper + n_bottom) containing the base forecast means

#'   for both upper and bottom levels (upper first, then bottom).

#' @param Psi_post Scale matrix (n_upper + n_bottom x n_upper + n_bottom) of the posterior

#'   multivariate t-distribution.

#' @param nu_post Degrees of freedom of the posterior multivariate t-distribution.

#' @param return_upper Logical, whether to return the reconciled parameters for the upper variables (default is FALSE).

#'

#' @return A list containing:

#' \itemize{

#'   \item \code{bottom_rec_mean}: reconciled bottom-level mean.

#'   \item \code{bottom_rec_scale_matrix}: reconciled bottom-level scale matrix.

#'   \item \code{bottom_rec_df}: reconciled degrees of freedom.

#' }

#' If \code{return_upper} is TRUE, also returns:

#' \itemize{

#'  \item \code{upper_rec_mean}: reconciled upper-level mean.

#'  \item \code{upper_rec_scale_matrix}: reconciled upper-level scale matrix.

#'  \item \code{upper_rec_df}: reconciled upper-level degrees of freedom.

#'  }

#'

#' @keywords internal

#' @export

.core_reconc_t <- function(A, base_fc_mean, Psi_post, nu_post, return_upper = FALSE,

                           return_parameters = FALSE) {

  # Indices for Upper and Bottom

  # Extract Psi blocks

  # Extract means

  # Compute Q = Psi_U - (Psi_UB %*% t(A)) - (A %*% t(Psi_UB)) + (A %*% Psi_B %*% t(A))

  # Invert Q using Cholesky (should be p.d.)

  # If Cholesky fails, stop

  }

  # Incoherence

  # Lambda = Psi_UB^T - Psi_B A^T

  # Compute b_tilde = b_hat + Lambda * Q^{-1} * delta

  # Compute nu_tilde = nu' - n_b + 1

  # Compute C = (1 + delta^T Q^{-1} delta) / nu_tilde

  # Compute Sigma_B_tilde = C * [ Psi_B - Lambda * Q^{-1} * Lambda^T ]

  # Prepare output

  )

    # Compute the parameters of the uppers using closure property

  }

}