From d12cd3e2885de8ce3947b214381c762fff635d81 Mon Sep 17 00:00:00 2001 From: maechler Date: Fri, 3 Jan 2025 21:33:24 +0000 Subject: [PATCH] fix dealing with erronous `seasonal = `; sync arima0() code & doc git-svn-id: https://svn.r-project.org/R/trunk@87520 00db46b3-68df-0310-9c12-caf00c1e9a41 --- doc/NEWS.Rd | 6 ++++- src/library/stats/R/arima.R | 23 ++++++++-------- src/library/stats/R/arma0.R | 48 ++++++++++++++++----------------- src/library/stats/man/arima.Rd | 11 ++++---- src/library/stats/man/arima0.Rd | 29 +++++++++++--------- tests/reg-tests-1e.R | 11 ++++++++ 6 files changed, 74 insertions(+), 54 deletions(-) diff --git a/doc/NEWS.Rd b/doc/NEWS.Rd index e030cdb7682..2fbe0260262 100644 --- a/doc/NEWS.Rd +++ b/doc/NEWS.Rd @@ -515,7 +515,11 @@ \item \code{dev.capabilities() $ events} now reports \code{"Idle"} if the device provides it, fixing \PR{18836}, thanks to \I{Trevor Davis}. - } + + \item \code{arima(.., seasonal = )} correctly errors + now, ditto for \code{arima0()}, thanks to \I{Norbert Kuder}'s report + on the R-devel list. + } } } diff --git a/src/library/stats/R/arima.R b/src/library/stats/R/arima.R index 6d9a8f3d86d..55ce8a7d8b6 100644 --- a/src/library/stats/R/arima.R +++ b/src/library/stats/R/arima.R @@ -1,7 +1,7 @@ # File src/library/stats/R/arima.R # Part of the R package, https://www.R-project.org # -# Copyright (C) 2002-2015 The R Core Team +# Copyright (C) 2002-2025 The R Core Team # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -52,7 +52,7 @@ arima <- function(x, order = c(0L, 0L, 0L), .Call(C_ARIMA_Like, y, mod, 0L, TRUE) } - ## the objective function called by optim() + ## the objective function called by optim(); using {coef, mask, arma, mod, ....} armafn <- function(p, trans) { par <- coef @@ -124,10 +124,11 @@ arima <- function(x, order = c(0L, 0L, 0L), if(!is.numeric(seasonal$order) || length(seasonal$order) != 3L || any(seasonal$order < 0L)) stop("'seasonal$order' must be a non-negative numeric vector of length 3") - } else if(is.numeric(order)) { - if(length(order) == 3L) seasonal <- list(order=seasonal) - else ("'seasonal' is of the wrong length") - } else stop("'seasonal' must be a list with component 'order'") + } else if(is.numeric(seasonal)) { # meant to be seasonal$order + if(length(seasonal) != 3L || any(seasonal < 0)) + stop("if not a list, 'seasonal' must be a non-negative numeric vector of length 3") + seasonal <- list(order=seasonal) + } else stop("'seasonal' is neither a list with component 'order' nor a numeric vector of length 3") if (is.null(seasonal$period) || is.na(seasonal$period) || seasonal$period == 0) seasonal$period <- frequency(x) @@ -290,14 +291,15 @@ arima <- function(x, order = c(0L, 0L, 0L), } } trarma <- .Call(C_ARIMA_transPars, init, arma, transform.pars) - mod <- makeARIMA(trarma[[1L]], trarma[[2L]], Delta, kappa, SSinit) res <- if(no.optim) list(convergence = 0, par = numeric(), value = armafn(numeric(), as.logical(transform.pars))) - else + else { + mod <- makeARIMA(trarma[[1L]], trarma[[2L]], Delta, kappa, SSinit) optim(init[mask], armafn, method = optim.method, hessian = TRUE, control = optim.control, trans = as.logical(transform.pars)) + } if(res$convergence > 0) warning(gettextf("possible convergence problem: optim gave code = %d", res$convergence), domain = NA) @@ -333,9 +335,8 @@ arima <- function(x, order = c(0L, 0L, 0L), } else var <- if(no.optim) numeric() else solve(res$hessian * n.used) trarma <- .Call(C_ARIMA_transPars, coef, arma, FALSE) mod <- makeARIMA(trarma[[1L]], trarma[[2L]], Delta, kappa, SSinit) - val <- if(ncxreg > 0L) - arimaSS(x - xreg %*% coef[narma + (1L:ncxreg)], mod) - else arimaSS(x, mod) + val <- arimaSS(if(ncxreg > 0L) x - xreg %*% coef[narma + (1L:ncxreg)] + else x, mod) sigma2 <- val[[1L]][1L]/n.used } value <- 2 * n.used * res$value + n.used + n.used * log(2 * pi) diff --git a/src/library/stats/R/arma0.R b/src/library/stats/R/arma0.R index 773303a75ca..df02a62fe8c 100644 --- a/src/library/stats/R/arma0.R +++ b/src/library/stats/R/arma0.R @@ -1,7 +1,7 @@ # File src/library/stats/R/arma0.R # Part of the R package, https://www.R-project.org # -# Copyright (C) 1999-2019 The R Core Team +# Copyright (C) 1999-2025 The R Core Team # # This program is free software; you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -16,8 +16,8 @@ # A copy of the GNU General Public License is available at # https://www.R-project.org/Licenses/ -arima0 <- function(x, order = c(0, 0, 0), - seasonal = list(order = c(0, 0, 0), period = NA), +arima0 <- function(x, order = c(0L, 0L, 0L), + seasonal = list(order = c(0L, 0L, 0L), period = NA), xreg = NULL, include.mean = TRUE, delta = 0.01, transform.pars = TRUE, fixed = NULL, init = NULL, method = c("ML", "CSS"), n.cond, @@ -41,7 +41,7 @@ arima0 <- function(x, order = c(0, 0, 0), { ## polyroot can't cope with leading zero. q <- length(ma) - q0 <- max(which(c(1,ma) != 0)) - 1 + q0 <- max(which(c(1,ma) != 0)) - 1L if(!q0) return(ma) roots <- polyroot(c(1, ma[1L:q0])) ind <- Mod(roots) < 1 @@ -55,9 +55,10 @@ arima0 <- function(x, order = c(0, 0, 0), } series <- deparse1(substitute(x)) - if(NCOL(x) > 1) + if(NCOL(x) > 1L) stop("only implemented for univariate time series") method <- match.arg(method) + x <- as.ts(x) if(!is.numeric(x)) stop("'x' must be numeric") @@ -72,12 +73,13 @@ arima0 <- function(x, order = c(0, 0, 0), if(is.null(seasonal$order)) stop("'seasonal' must be a list with component 'order'") if(!is.numeric(seasonal$order) || length(seasonal$order) != 3L - || any(seasonal$order < 0)) + || any(seasonal$order < 0L)) stop("'seasonal$order' must be a non-negative numeric vector of length 3") - } else if(is.numeric(order)) { - if(length(order) == 3) seasonal <- list(order=seasonal) - else ("'seasonal' is of the wrong length") - } else stop("'seasonal' must be a list with component 'order'") + } else if(is.numeric(seasonal)) { # meant to be seasonal$order + if(length(seasonal) != 3L || any(seasonal < 0)) + stop("if not a list, 'seasonal' must be a non-negative numeric vector of length 3") + seasonal <- list(order=seasonal) + } else stop("'seasonal' is neither a list with component 'order' nor a numeric vector of length 3") if(is.null(seasonal$period) || is.na(seasonal$period) || seasonal$period == 0) seasonal$period <- frequency(x) @@ -167,7 +169,7 @@ arima0 <- function(x, order = c(0, 0, 0), if(!arCheck(init[1L:arma[1L]])) stop("non-stationary AR part") if(arma[3L] > 0) - if(!arCheck(init[sum(arma[1L:2]) + 1L:arma[3L]])) + if(!arCheck(init[sum(arma[1L:2L]) + 1L:arma[3L]])) stop("non-stationary seasonal AR part") ## enforce invertibility if(arma[2L] > 0) { @@ -175,7 +177,7 @@ arima0 <- function(x, order = c(0, 0, 0), init[ind] <- maInvert(init[ind]) } if(arma[4L] > 0) { - ind <- sum(arma[1L:3]) + 1L:arma[4L] + ind <- sum(arma[1L:3L]) + 1L:arma[4L] init[ind] <- maInvert(init[ind]) } init <- .Call(C_Invtrans, G, as.double(init)) @@ -210,12 +212,12 @@ arima0 <- function(x, order = c(0, 0, 0), class(resid) <- "ts" n.used <- sum(!is.na(resid)) nm <- NULL - if(arma[1L] > 0) nm <- c(nm, paste0("ar", 1L:arma[1L])) - if(arma[2L] > 0) nm <- c(nm, paste0("ma", 1L:arma[2L])) - if(arma[3L] > 0) nm <- c(nm, paste0("sar", 1L:arma[3L])) - if(arma[4L] > 0) nm <- c(nm, paste0("sma", 1L:arma[4L])) + if(arma[1L] > 0L) nm <- c(nm, paste0("ar", 1L:arma[1L])) + if(arma[2L] > 0L) nm <- c(nm, paste0("ma", 1L:arma[2L])) + if(arma[3L] > 0L) nm <- c(nm, paste0("sar", 1L:arma[3L])) + if(arma[4L] > 0L) nm <- c(nm, paste0("sma", 1L:arma[4L])) fixed[mask] <- coef - if(ncxreg > 0) { + if(ncxreg > 0L) { nm <- c(nm, cn) if(!orig.xreg) { ind <- narma + 1L:ncxreg @@ -255,11 +257,9 @@ print.arima0 <- function(x, digits = max(3L, getOption("digits") - 3L), print.default(coef, print.gap = 2) cm <- x$call$method if(is.null(cm) || cm != "CSS") - cat("\nsigma^2 estimated as ", - format(x$sigma2, digits = digits), - ": log likelihood = ", format(round(x$loglik,2)), - ", aic = ", format(round(x$aic,2)), - "\n", sep = "") + cat("\nsigma^2 estimated as ", format(x$sigma2, digits = digits), + ": log likelihood = ", format(round(x$loglik, 2L)), + ", aic = ", format(round(x$aic, 2L)), "\n", sep = "") else cat("\nsigma^2 estimated as ", format(x$sigma2, digits = digits), @@ -294,12 +294,12 @@ predict.arima0 <- xm <- drop(as.matrix(newxreg) %*% coefs[-(1L:narma)]) } else xm <- 0 ## check invertibility of MA part(s) - if(arma[2L] > 0) { + if(arma[2L] > 0L) { ma <- coefs[arma[1L] + 1L:arma[2L]] if(any(Mod(polyroot(c(1, ma))) < 1)) warning("MA part of model is not invertible") } - if(arma[4L] > 0) { + if(arma[4L] > 0L) { ma <- coefs[sum(arma[1L:3L]) + 1L:arma[4L]] if(any(Mod(polyroot(c(1, ma))) < 1)) warning("seasonal MA part of model is not invertible") diff --git a/src/library/stats/man/arima.Rd b/src/library/stats/man/arima.Rd index c996ba36157..653bee2ee08 100644 --- a/src/library/stats/man/arima.Rd +++ b/src/library/stats/man/arima.Rd @@ -1,6 +1,6 @@ % File src/library/stats/man/arima.Rd % Part of the R package, https://www.R-project.org -% Copyright 1995-2024 R Core Team +% Copyright 1995-2025 R Core Team % Distributed under GPL 2 or later \name{arima} @@ -24,13 +24,13 @@ arima(x, order = c(0L, 0L, 0L), \arguments{ \item{x}{a univariate time series} - \item{order}{A specification of the non-seasonal part of the ARIMA + \item{order}{a specification of the non-seasonal part of the ARIMA model: the three integer components \eqn{(p, d, q)} are the AR order, the degree of differencing, and the MA order.} - \item{seasonal}{A specification of the seasonal part of the ARIMA + \item{seasonal}{a specification of the seasonal part of the ARIMA model, plus the period (which defaults to \code{frequency(x)}). - This may be a list with components \code{order} and + This may be a \code{\link{list}} with components \code{order} and \code{period}, or just a numeric vector of length 3 which specifies the seasonal \code{order}. In the latter case the default period is used.} @@ -38,7 +38,8 @@ arima(x, order = c(0L, 0L, 0L), \item{xreg}{Optionally, a vector or matrix of external regressors, which must have the same number of rows as \code{x}.} - \item{include.mean}{Should the ARMA model include a mean/intercept term? The + \item{include.mean}{logical indicating if the ARMA model should include a + mean/intercept term. The default is \code{TRUE} for undifferenced series, and it is ignored for ARIMA models with differencing.} diff --git a/src/library/stats/man/arima0.Rd b/src/library/stats/man/arima0.Rd index 5289b9354e5..9d79258843a 100644 --- a/src/library/stats/man/arima0.Rd +++ b/src/library/stats/man/arima0.Rd @@ -1,6 +1,6 @@ % File src/library/stats/man/arima0.Rd % Part of the R package, https://www.R-project.org -% Copyright 1995-2024 R Core Team +% Copyright 1995-2025 R Core Team % Distributed under GPL 2 or later \name{arima0} @@ -12,6 +12,7 @@ \description{ Fit an ARIMA model to a univariate time series, and forecast from the fitted model. + For new projects, consider using \code{\link{arima}()} instead. } \usage{ arima0(x, order = c(0, 0, 0), @@ -25,16 +26,16 @@ arima0(x, order = c(0, 0, 0), \arguments{ \item{x}{a univariate time series} - \item{order}{A specification of the non-seasonal part of the ARIMA - model: the three components \eqn{(p, d, q)} are the AR order, the + \item{order}{a specification of the non-seasonal part of the ARIMA + model: the three integer components \eqn{(p, d, q)} are the AR order, the degree of differencing, and the MA order.} - \item{seasonal}{A specification of the seasonal part of the ARIMA + \item{seasonal}{a specification of the seasonal part of the ARIMA model, plus the period (which defaults to \code{frequency(x)}). - This should be a list with components \code{order} and - \code{period}, but a specification of just a numeric vector of - length 3 will be turned into a suitable list with the specification - as the \code{order}.} + This may be a \code{\link{list}} with components \code{order} and + \code{period}, or just a numeric vector of length 3 which + specifies the seasonal \code{order}. In the latter case the + default period is used.} \item{xreg}{Optionally, a vector or matrix of external regressors, which must have the same number of rows as \code{x}.} @@ -47,7 +48,7 @@ arima0(x, order = c(0, 0, 0), \item{delta}{A value to indicate at which point \sQuote{fast recursions} should be used. See the \sQuote{Details} section.} - \item{transform.pars}{Logical. If true, the AR parameters are + \item{transform.pars}{logical; if true, the AR parameters are transformed to ensure that they remain in the region of stationarity. Not used for \code{method = "CSS"}.} @@ -62,10 +63,10 @@ arima0(x, order = c(0, 0, 0), regression coefficients. Values already specified in \code{fixed} will be ignored.} - \item{method}{Fitting method: maximum likelihood or minimize + \item{method}{fitting method: maximum likelihood or minimize conditional sum-of-squares. Can be abbreviated.} - \item{n.cond}{Only used if fitting by conditional-sum-of-squares: the + \item{n.cond}{only used if fitting by conditional-sum-of-squares: the number of initial observations to ignore. It will be ignored if less than the maximum lag of an AR term.} @@ -204,7 +205,7 @@ arima0(x, order = c(0, 0, 0), Harvey, A. C. and McKenzie, C. R. (1982). Algorithm AS 182: An algorithm for finite sample prediction from ARIMA - processes. + processes. \emph{Applied Statistics}, \bold{31}, 180--187. \doi{10.2307/2347987}. @@ -216,7 +217,9 @@ arima0(x, order = c(0, 0, 0), } \note{ - This is a preliminary version, and will be replaced by \code{\link{arima}}. + This has been a preliminary version, and is mostly replaced by + \code{\link{arima}}, notably in the presence of missing values. + \code{arima0()} remains mostly for reproducibility reasons. The standard errors of prediction exclude the uncertainty in the estimation of the ARMA model and the regression coefficients. diff --git a/tests/reg-tests-1e.R b/tests/reg-tests-1e.R index 132a72f25a8..0837f712226 100644 --- a/tests/reg-tests-1e.R +++ b/tests/reg-tests-1e.R @@ -1712,6 +1712,17 @@ if(length(iLA) && nzchar(La_version())) { cat("sessionInfo - La_* checking: ") ## the "LAPACK: .." was entirely empty when si$LAPACK was "" +## arima(*, seasonal = ) +(m <- tryCmsg( arima(presidents, order=c(2,0,1), seasonal=c(1, 0)) )) +stopifnot(exprs = { + grepl("'seasonal'", m, fixed=TRUE) + !englishMsgs || + grepl("must be a non-negative numeric vector", m, fixed=TRUE) +}) +## gave solve.default() error (as wrong model failed fitting) + + + ## keep at end rbind(last = proc.time() - .pt, total = proc.time())