[Rd] Fix for bug in arima function

[Martin Maechler]

> I noticed that the 3.2.1 release cycle is about to start.  Is there any
> chance that this fix will make it into the next version of R?
> 
> This bug is fairly serious: getting the wrong variance estimate leads to
> the wrong log-likelihood and the wrong AIC, BIC etc, which can and does
> lead to suboptimal model selection.  If it's not fixed, this issue will
> affect every student taking our time series course in Fall 2015 (and
> probably lots of other students in other time series courses).  When I
> taught time series in Spring 2015, I had to teach students how to work
> around the bug, which wasted class time and shook student confidence in R.
> It'd be great if we didn’t have to deal with this issue next semester.
> 
> Again, the fix is trivial:
> 
> --- a/src/library/stats/R/arima.R
> +++ b/src/library/stats/R/arima.R
> @@ -211,8 +211,10 @@ arima <- function(x, order = c(0L, 0L, 0L),
>          if(fit$rank == 0L) {
>              ## Degenerate model. Proceed anyway so as not to break old code
>              fit <- lm(x ~ xreg - 1, na.action = na.omit)
> +            n.used <- sum(!is.na(resid(fit))) - length(Delta)
> +        } else {
> +            n.used <- sum(!is.na(resid(fit)))
>          }
> -        n.used <- sum(!is.na(resid(fit))) - length(Delta)
>          init0 <- c(init0, coef(fit))
>          ses <- summary(fit)$coefficients[, 2L]
>          parscale <- c(parscale, 10 * ses)
> 

Yes, such a change *is* small in the source code.
But we have to be sure about  its desirability.

In another post about this you mention "REML", and I think we
really are discussing if variance estimates should use a
denominator of  'n'  or 'n - p' in this case.


> The patch that introduced the bug (
> https://github.com/wch/r-source/commit/32f633885a903bc422537dc426644f743cc645e0
> ) was designed to change the initialization for the optimization routine.

> The proposed fix leaves the deliberate part of the patch unchanged (it
> preserves the value of "init0").

I can confirm this... a change introduced in R 3.0.2.

I'm about to commit changes ... after also adding a proper
regression test.

Martin Maechler, ETH Zurich