Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

A bug fix in the forecast.gts function #53

Open
wants to merge 12 commits into
base: master
Choose a base branch
from
Open

A bug fix in the forecast.gts function #53

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
9bfeb3b
allow.changes argument is added.
ShanikaLW Mar 26, 2020
07a80bf
Block principal pivoting algorithm for OLS and WLS
ShanikaLW Mar 26, 2020
733b725
Block principal pivoting algorithm for MinT
ShanikaLW Mar 26, 2020
e1c40dd
A helper function for bpv
ShanikaLW Mar 26, 2020
10de90f
Implements non-negative forecast reconciliation
ShanikaLW Mar 26, 2020
4b3c31e
Implements non-negative forecast reconciliation
ShanikaLW Mar 26, 2020
dd1b3cf
Implements non-negative forecast reconciliation
ShanikaLW Mar 26, 2020
6cfeeec
Encoding UTF-8 is added
ShanikaLW Mar 26, 2020
c9081e3
Reference style is updated
ShanikaLW Mar 26, 2020
c9d9d40
A minor change
ShanikaLW Mar 26, 2020
81e563d
updated .Rd files
ShanikaLW Mar 26, 2020
5f9845f
A bug fixed
ShanikaLW Apr 18, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion DESCRIPTION
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,5 +39,6 @@ URL: http://pkg.earo.me/hts
BugReports: https://github.com/earowang/hts/issues
License: GPL (>= 2)
VignetteBuilder: knitr
RoxygenNote: 6.0.1
RoxygenNote: 7.1.0
Roxygen: list(markdown = TRUE, roclets=c('rd', 'collate', 'namespace'))
Encoding: UTF-8
278 changes: 170 additions & 108 deletions R/MinT.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,30 +35,16 @@ shrink.estim <- function(x, tar)
round(lambda, digits = 4))))
}

## Arguments

#fcasts: Matrix of forecasts for all levels of the hierarchical time series.
# Each row represents one forecast horizon and each column represents one time series from the hierarchy

# nodes: If the object class is hts, a list contains the number of child nodes referring to hts.
# groups: If the object is gts, a gmatrix is required, which is the same as groups in the function gts.
# residuals: Matrix of insample residuals for all time series in the hierarchy. Each column referring to one time series.
# covariance: Type of the covariance matrix to be used. Sample covariance matrix ("sam") or shrinking towards a diagonal unequal variances ("shr").
# algorithms: Algorithm used to compute inverse of the matrices.
# keep: Return a gts object or the reconciled forecasts at the bottom level.


# MinT - Trace minimization approach


#' Trace minimization for hierarchical or grouped time series
#'
#' Using the method of Wickramasuriya et al. (2015), this function combines the
#' Using the method of Wickramasuriya et al. (2019), this function combines the
#' forecasts at all levels of a hierarchical or grouped time series. The
#' \code{\link{forecast.gts}} calls this function when the \code{MinT} method
#' is selected.
#'
#'
#' @param fcasts Matrix of forecasts for all levels of a hierarchical or
#' grouped time series. Each row represents one forecast horizon and each
#' column represents one time series of aggregated or disaggregated forecasts.
Expand All @@ -70,23 +56,36 @@ shrink.estim <- function(x, tar)
#' disaggregated time series. The columns must be in the same order as
#' \code{fcasts}.
#' @param covariance Type of the covariance matrix to be used. Shrinking
#' towards a diagonal unequal variances ("shr") or sample covariance matrix
#' ("sam").
#' towards a diagonal unequal variances (\code{"shr"}) or sample covariance matrix
#' (\code{"sam"}).
#' @param nonnegative Logical. Should the reconciled forecasts be non-negative?
#' @param algorithms Algorithm used to compute inverse of the matrices.
#' @param keep Return a gts object or the reconciled forecasts at the bottom
#' level.
#' @param parallel Logical. Import parallel package to allow parallel processing.
#' @param num.cores Numeric. Specify how many cores are going to be used.
#' @param control.nn A list of control parameters to be passed on to the
#' block principal pivoting algorithm. See 'Details'.
#' @return Return the reconciled \code{gts} object or forecasts at the bottom
#' level.
#' @details
#' The \code{control.nn} argument is a list that can supply any of the following components:
#' \describe{
#' \item{\code{ptype}}{Permutation method to be used: \code{"fixed"} or \code{"random"}. Defaults to \code{"fixed"}.}
#' \item{\code{par}}{The number of full exchange rules that may be tried. Defaults to 10.}
#' \item{\code{gtol}}{The tolerance of the convergence criteria. Defaults to \code{sqrt(.Machine$double.eps)}.}
#' }
#' @author Shanika L Wickramasuriya
#' @seealso \code{\link[hts]{hts}}, \code{\link[hts]{gts}},
#' \code{\link[hts]{forecast.gts}}, \code{\link[hts]{combinef}}
#' @references Wickramasuriya, S. L., Athanasopoulos, G., & Hyndman, R. J.
#' (2015). Forecasting hierarchical and grouped time series through trace
#' minimization. \emph{Working paper 15/15, Department of Econometrics &
#' Business Statistics, Monash University.}
#' \url{http://robjhyndman.com/working-papers/mint/}
#' @references Wickramasuriya, S. L., Athanasopoulos, G., & Hyndman, R. J. (2019).
#' Optimal forecast reconciliation for hierarchical and grouped time series through trace minimization.
#' \emph{Journal of the American Statistical Association}, \bold{114}(526), 804--819. \url{http://robjhyndman.com/working-papers/mint/}
#'
#' Hyndman, R. J., Lee, A., & Wang, E. (2015). Fast computation of reconciled
#' Wickramasuriya, S. L., Turlach, B. A., & Hyndman, R. J. (to appear). Optimal non-negative forecast reconciliation.
#' \emph{Statistics and Computing}. \url{https://robjhyndman.com/publications/nnmint/}
#'
#' Hyndman, R. J., Lee, A., & Wang, E. (2016). Fast computation of reconciled
#' forecasts for hierarchical and grouped time series. \emph{Computational
#' Statistics and Data Analysis}, \bold{97}, 16--32.
#' \url{http://robjhyndman.com/papers/hgts/}
Expand All @@ -113,6 +112,20 @@ shrink.estim <- function(x, tar)
#' y.f_cg <- MinT(allf, get_nodes(htseg1), residual = res, covariance = "shr",
#' keep = "all", algorithms = "cg")
#' }
#'
#' \dontrun{h <- 12
#' ally <- abs(aggts(htseg2))
#' allf <- matrix(NA, nrow = h, ncol = ncol(ally))
#' res <- matrix(NA, nrow = nrow(ally), ncol = ncol(ally))
#' for(i in 1:ncol(ally)) {
#' fit <- auto.arima(ally[, i], lambda = 0, biasadj = TRUE)
#' allf[,i] <- forecast(fit, h = h)$mean
#' res[,i] <- na.omit(ally[, i] - fitted(fit))
#' }
#' b.f <- MinT(allf, get_nodes(htseg2), residual = res, covariance = "shr", keep = "bottom", algorithms = "lu")
#' b.nnf <- MinT(allf, get_nodes(htseg2), residual = res, covariance = "shr", keep = "bottom", algorithms = "lu",
#' nonnegative = TRUE, parallel = TRUE)
#' }
#'
#' # gts example
#' \dontrun{abc <- ts(5 + matrix(sort(rnorm(200)), ncol = 4, nrow = 50))
Expand All @@ -136,126 +149,175 @@ shrink.estim <- function(x, tar)
#' plot(y.f)}
#'
#' @export MinT
MinT <- function (fcasts, nodes, groups, residual, covariance = c("shr", "sam"),
algorithms = c("lu", "cg", "chol"), keep = c("gts", "all", "bottom"))
MinT <- function (fcasts, nodes = NULL, groups = NULL, residual, covariance = c("shr", "sam"),
nonnegative = FALSE, algorithms = c("lu", "cg", "chol"),
keep = c("gts", "all", "bottom"), parallel = FALSE, num.cores = 2, control.nn = list())
{
if (is.null(nodes) && is.null(groups)) {
stop("Please specify the hierarchical or the grouping structure.", call. = FALSE)
}

if (!xor(is.null(nodes), is.null(groups))) {
stop("Please specify either nodes or groups argument, not both.", call. = FALSE)
}

alg <- match.arg(algorithms)
keep <- match.arg(keep)
covar <- match.arg(covariance)
res <- residual
fcasts <- stats::as.ts(fcasts)
tspx <- stats::tsp(fcasts)
cnames <- colnames(fcasts)

if(missing(residual))
{
stop("MinT needs insample residuals.", call. = FALSE)
}
if(covar=="sam")
{
n <- nrow(res)
w.1 <- crossprod(res) / n
if(is.positive.definite(w.1)==FALSE)

if (!nonnegative) {
if (missing(residual))
{
stop("MinT needs covariance matrix to be positive definite.", call. = FALSE)
stop("MinT needs insample residuals.", call. = FALSE)
}
}else{
tar <- lowerD(res)
shrink <- shrink.estim(res, tar)
w.1 <- shrink[[1]]
lambda <- shrink[[2]]
if(is.positive.definite(w.1)==FALSE)
if (covar=="sam")
{
stop("MinT needs covariance matrix to be positive definite.", call. = FALSE)
}
}

if (missing(groups)) { # hts class
totalts <- sum(Mnodes(nodes))
if (!is.matrix(fcasts)) {
fcasts <- t(fcasts)
}
h <- nrow(fcasts)
if (ncol(fcasts) != totalts) {
stop("Argument fcasts requires all the forecasts.", call. = FALSE)
n <- nrow(res)
w.1 <- crossprod(res) / n
if(is.positive.definite(w.1)==FALSE)
{
stop("MinT needs covariance matrix to be positive definite.", call. = FALSE)
}
} else {
tar <- lowerD(res)
shrink <- shrink.estim(res, tar)
w.1 <- shrink[[1]]
lambda <- shrink[[2]]
if (is.positive.definite(w.1)==FALSE)
{
stop("MinT needs covariance matrix to be positive definite.", call. = FALSE)
}
}
gmat <- GmatrixH(nodes)
fcasts <- t(fcasts)
if (alg == "chol") {
smat <- Smatrix(gmat)
if (!is.null(w.1)) {
w.1 <- as.matrix.csr(w.1)

if (is.null(groups)) { # hts class
totalts <- sum(Mnodes(nodes))
if (!is.matrix(fcasts)) {
fcasts <- t(fcasts)
}
allf <- CHOL(fcasts = fcasts, S = smat, weights = w.1)
h <- nrow(fcasts)
if (ncol(fcasts) != totalts) {
stop("Argument fcasts requires all the forecasts.", call. = FALSE)
}
else {
smat <- SmatrixM(gmat)
gmat <- GmatrixH(nodes)
fcasts <- t(fcasts)
if (alg == "chol") {
smat <- Smatrix(gmat)
if (!is.null(w.1)) {
w.1 <- as.matrix.csr(w.1)
}
allf <- CHOL(fcasts = fcasts, S = smat, weights = w.1, allow.changes = FALSE)
}
else {
smat <- SmatrixM(gmat)
if (!is.null(w.1)) {
weights <- methods::as(w.1, "sparseMatrix")
}
if (alg == "lu") {
allf <- LU(fcasts = fcasts, S = smat, weights = weights)
allf <- LU(fcasts = fcasts, S = smat, weights = weights, allow.changes = FALSE)
}
else if (alg == "cg") {
allf <- CG(fcasts = fcasts, S = smat, weights = weights)
allf <- CG(fcasts = fcasts, S = smat, weights = weights, allow.changes = FALSE)
}
}

if (keep == "all") {
if (keep == "all") {
out <- t(allf)
}
else {
bottom <- totalts - (ncol(smat):1L) + 1L
bf <- t(allf[bottom, ])
if (keep == "gts") {
bf <- ts(bf, start = tspx[1L], frequency = tspx[3L])
out <- suppressMessages(hts(bf, nodes = nodes))
}
else {
out <- bf
bottom <- totalts - (ncol(smat):1L) + 1L
bf <- t(allf[bottom, ])
if (keep == "gts") {
bf <- ts(bf, start = tspx[1L], frequency = tspx[3L])
out <- suppressMessages(hts(bf, nodes = nodes))
}
else {
out <- bf
}
}
}
}
else if (missing(nodes)) {
rownames(groups) <- NULL
gmat <- GmatrixG(groups)
totalts <- sum(Mlevel(gmat))
if (ncol(fcasts) != totalts) {
stop("Argument fcasts requires all the forecasts.", call. = FALSE)
}
fcasts <- t(fcasts)
if (alg == "chol") {
smat <- Smatrix(gmat)
if (!is.null(w.1)) {
weights <- as.matrix.csr(w.1)
else if (is.null(nodes)) {
rownames(groups) <- NULL
gmat <- GmatrixG(groups)
totalts <- sum(Mlevel(gmat))
if (ncol(fcasts) != totalts) {
stop("Argument fcasts requires all the forecasts.", call. = FALSE)
}
allf <- CHOL(fcasts = fcasts, S = smat, weights = weights)
}
else {
smat <- SmatrixM(gmat)
if (!is.null(w.1)) {
weights <- methods::as(w.1, "sparseMatrix")
fcasts <- t(fcasts)
if (alg == "chol") {
smat <- Smatrix(gmat)
if (!is.null(w.1)) {
weights <- as.matrix.csr(w.1)
}
allf <- CHOL(fcasts = fcasts, S = smat, weights = weights, allow.changes = FALSE)
}
if (alg == "lu") {
allf <- LU(fcasts = fcasts, S = smat, weights = weights)
else {
smat <- SmatrixM(gmat)
if (!is.null(w.1)) {
weights <- methods::as(w.1, "sparseMatrix")
}
if (alg == "lu") {
allf <- LU(fcasts = fcasts, S = smat, weights = weights, allow.changes = FALSE)
}
else if (alg == "cg") {
allf <- CG(fcasts = fcasts, S = smat, weights = weights, allow.changes = FALSE)
}
}
else if (alg == "cg") {
allf <- CG(fcasts = fcasts, S = smat, weights = weights)
if (keep == "all") {
out <- t(allf)
}
else {
bottom <- totalts - (ncol(smat):1L) + 1L
bf <- t(allf[bottom, ])
if (keep == "gts") {
colnames(bf) <- cnames[bottom]
bf <- ts(bf, start = tspx[1L], frequency = tspx[3L])
out <- suppressMessages(gts(bf, groups = groups))
}
else {
out <- bf
}
}
}
if (keep == "all") {
out <- t(allf)
} else {
if (any(fcasts < 0)) {
fcasts[fcasts < 0] <- 0
warning("Negative base forecasts are truncated to zero.")
}
else {
bottom <- totalts - (ncol(smat):1L) + 1L
bf <- t(allf[bottom, ])
if (keep == "gts") {
colnames(bf) <- cnames[bottom]
bf <- ts(bf, start = tspx[1L], frequency = tspx[3L])
out <- suppressMessages(gts(bf, groups = groups))

lst.fc <- split(fcasts, row(fcasts))
if (parallel) {
if (is.null(num.cores)) {
num.cores <- detectCores()
}
else {
cl <- makeCluster(num.cores)
bf <- parSapplyLB(cl = cl, X = lst.fc, MinTbpv, nodes = nodes, groups = groups, res = res, covar = covar, alg = alg, control.nn = control.nn, simplify = TRUE)
stopCluster(cl = cl)
} else {
bf <- sapply(lst.fc, MinTbpv, nodes = nodes, groups = groups, res = res, covar = covar, alg = alg, control.nn = control.nn)
}
bf <- ts(t(bf), start = tspx[1L], frequency = tspx[3L])
if (is.null(groups)) {
if (keep == "bottom") {
out <- bf
} else {
out <- suppressMessages(hts(bf, nodes = nodes))
if (keep == "all") {
out <- aggts(out)
}
}
} else {
if (keep == "bottom") {
out <- bf
} else {
colnames(bf) <- tail(cnames, ncol(bf))
out <- suppressMessages(gts(bf, groups = groups))
if (keep == "all") {
out <- aggts(out)
}
}
}
}
Expand Down
Loading