[Rd] Fix for bug in arima function

Martin Maechler maechler at stat.math.ethz.ch
Thu May 28 10:17:23 CEST 2015


>>>>> Patrick Perry <pperry at stern.nyu.edu>
>>>>>     on Wed, 27 May 2015 23:19:09 -0400 writes:

 {@PP, you forgot this part:}
 >>>>> peter dalgaard <pdalgd at gmail.com>
 >>>>>     on Thu, 21 May 2015 14:36:03 +0200 writes:

    >> I suspect that what we really need is  
    >> 
    >> fitI <- lm(x ~ xreg - 1, na.action = na.omit)
    >> fit <- if(length(dx) > ncol(dxreg))
    >>     lm(dx ~ dxreg - 1, na.action = na.omit)
    >> else list(rank = 0L)
    >> if(fit$rank == 0L) {
    >>    ## Degenerate model. Proceed anyway so as not to break old code
    >>    fit <- fitI
    >> }
    >> n.used <- sum(!is.na(resid(fitI))) - length(Delta)
    >> init0 <- c(init0, coef(fit))
    >> 
    >> At least that would be the conservative change to get n.used indentical to what it was in 3.0.1

(Sorry for not taking up the thread ..)
That's definitely conservative and hence safest from that point
of view.  

On the other hand, to me, it did look a bit strange or "ugly" to always
perform to different lm()s and the coefficients of one and the
residuals of the other.

    > Along the same lines, here’s a solution that avoids the extra call to lm:

    > fit <- if(length(dx) > ncol(dxreg))
    >    lm(dx ~ dxreg - 1, na.action = na.omit)
    > else list(rank = 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)
    > }
    > isna <- apply(is.na(xreg), 1, any) | is.na(x)
    > n.used <- sum(!isna) - length(Delta)
    > init0 <- c(init0, coef(fit))

That is indeed nicer ... and with logic closer to the current
code. I have very slightly changed it, e.g., using  anyNA(.),
and tested it with some more examples... and this does look good
to me in the sense that it is "internally more consistent".

To get this going and "exposed to CRAN", I'm committing it
(with the regression tests, and other necessary "entries") to
R-devel only, but with the intent to port to R-patched in a
couple of days.

Martin



More information about the R-devel mailing list