Skip to content

Commit

Permalink
fix dealing with erronous seasonal = <numeric>; sync arima0() code …
Browse files Browse the repository at this point in the history
…& doc

git-svn-id: https://svn.r-project.org/R/trunk@87520 00db46b3-68df-0310-9c12-caf00c1e9a41
  • Loading branch information
maechler committed Jan 3, 2025
1 parent 48ab7b4 commit d12cd3e
Show file tree
Hide file tree
Showing 6 changed files with 74 additions and 54 deletions.
6 changes: 5 additions & 1 deletion doc/NEWS.Rd
Original file line number Diff line number Diff line change
Expand Up @@ -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 = <wrong-vector>)} correctly errors
now, ditto for \code{arima0()}, thanks to \I{Norbert Kuder}'s report
on the R-devel list.
}
}
}

Expand Down
23 changes: 12 additions & 11 deletions src/library/stats/R/arima.R
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
48 changes: 24 additions & 24 deletions src/library/stats/R/arma0.R
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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,
Expand All @@ -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
Expand All @@ -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")
Expand All @@ -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)
Expand Down Expand Up @@ -167,15 +169,15 @@ 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) {
ind <- arma[1L] + 1L:arma[2L]
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))
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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),
Expand Down Expand Up @@ -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")
Expand Down
11 changes: 6 additions & 5 deletions src/library/stats/man/arima.Rd
Original file line number Diff line number Diff line change
@@ -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}
Expand All @@ -24,21 +24,22 @@ 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.}

\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.}

Expand Down
29 changes: 16 additions & 13 deletions src/library/stats/man/arima0.Rd
Original file line number Diff line number Diff line change
@@ -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}
Expand All @@ -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),
Expand All @@ -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}.}
Expand All @@ -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"}.}

Expand All @@ -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.}

Expand Down Expand Up @@ -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}.

Expand All @@ -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.
Expand Down
11 changes: 11 additions & 0 deletions tests/reg-tests-1e.R
Original file line number Diff line number Diff line change
Expand Up @@ -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 = <numeric>)
(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())

0 comments on commit d12cd3e

Please sign in to comment.