Bug in arima (sigma2 calculation)

arima-bug

# There is currently a bug in the arima function. The innovation variance
# estimate uses the wrong denominator whenever xreg is non-null.  This bug was
# introduced in the patch displayed here:
#
# https://github.com/wch/r-source/commit/32f633885a903bc422537dc426644f743cc645e0
#
# The fix is very simple:
#
# https://bugs.r-project.org/bugzilla3/show_bug.cgi?id=16278#c4  (comment 4).
#
# The following script demonstrates the bug.
#

testcase <- function(d, D, period=2, n=10)
{
    # generate data
    x <- x0 <- rnorm(n)
    xreg <- xreg0 <- rep(1, n)
    if (d > 0) {
        x <- diffinv(x, 1, d)
        xreg <- diffinv(xreg, 1, d)
    }
    if (D > 0) {
        x <- diffinv(x, period, D)
        xreg <- diffinv(xreg, period, D)
    }

    # compute expected ss, sigma2
    ss <- sum((x0 - mean(x0))^2)
    expected.sigma2 <- ss / n

    # fit via arima0
    fit0 <- arima0(x, c(0,d,0), list(order=c(0,D,0), period=period),
                   xreg=xreg, include.mean=FALSE)
    arima0.sigma2 <- fit0$sigma2
    arima0.n.used <- round(ss / fit0$sigma2)

    # fit via arima
    fit1 <- arima(x, c(0,d,0), list(order=c(0,D,0), period=period), xreg=xreg,
                  include.mean=FALSE, method="ML")
    arima.sigma2 <- fit1$sigma2
    arima.n.used <- round(ss / fit1$sigma2)

    data.frame(d, D, period, n, expected.sigma2, arima0.n.used,  arima0.sigma2,
               arima.n.used, arima.sigma2)
}


set.seed(0)
res <- do.call("rbind",
               apply(data.frame(d = rep(0:3, 4), D = rep(0:3, each=4)), 1,
                     function(x) testcase(x[["d"]], x[["D"]])))
print(res)


# Without proposed fix:

##    d D period  n expected.sigma2 arima0.n.used arima0.sigma2 arima.n.used arima.sigma2
## 1  0 0      2 10       1.3075522            10     1.3075522           10    1.3075522
## 2  1 0      2 10       0.4147800            10     0.4147800            9    0.4608667
## 3  2 0      2 10       0.4374619            10     0.4374619            8    0.5468274
## 4  3 0      2 10       0.5358400            10     0.5358400            7    0.7654858
## 5  0 1      2 10       1.1225139            10     1.1225139            8    1.4031423
## 6  1 1      2 10       1.3653702            10     1.3653702            7    1.9505289
## 7  2 1      2 10       0.3569846            10     0.3569846            6    0.5949744
## 8  3 1      2 10       0.4977978            10     0.4977978            5    0.9955956
## 9  0 2      2 10       0.5496022            10     0.5496022            6    0.9160037
## 10 1 2      2 10       0.6129431            10     0.6129431            5    1.2258861
## 11 2 2      2 10       1.0540106            10     1.0540106            4    2.6350265
## 12 3 2      2 10       0.6235336            10     0.6235336            3    2.0784453
## 13 0 3      2 10       0.9345655            10     0.9345655            4    2.3364138
## 14 1 3      2 10       0.8364767            10     0.8364767            3    2.7882558
## 15 2 3      2 10       0.4359491            10     0.4359491            2    2.1797453
## 16 3 3      2 10       1.1016849            10     1.1016849            1   11.0168459



# With proposed fix:

##    d D period  n expected.sigma2 arima0.n.used arima0.sigma2 arima.n.used arima.sigma2
## 1  0 0      2 10       1.3075522            10     1.3075522           10    1.3075522
## 2  1 0      2 10       0.4147800            10     0.4147800           10    0.4147800
## 3  2 0      2 10       0.4374619            10     0.4374619           10    0.4374619
## 4  3 0      2 10       0.5358400            10     0.5358400           10    0.5358400
## 5  0 1      2 10       1.1225139            10     1.1225139           10    1.1225139
## 6  1 1      2 10       1.3653702            10     1.3653702           10    1.3653702
## 7  2 1      2 10       0.3569846            10     0.3569846           10    0.3569846
## 8  3 1      2 10       0.4977978            10     0.4977978           10    0.4977978
## 9  0 2      2 10       0.5496022            10     0.5496022           10    0.5496022
## 10 1 2      2 10       0.6129431            10     0.6129431           10    0.6129431
## 11 2 2      2 10       1.0540106            10     1.0540106           10    1.0540106
## 12 3 2      2 10       0.6235336            10     0.6235336           10    0.6235336
## 13 0 3      2 10       0.9345655            10     0.9345655           10    0.9345655
## 14 1 3      2 10       0.8364767            10     0.8364767           10    0.8364767
## 15 2 3      2 10       0.4359491            10     0.4359491           10    0.4359491
## 16 3 3      2 10       1.1016849            10     1.1016849           10    1.1016846