Seasonal ARIMA – Constraints on the Coefficients and Possible Software Bug in ITSM

arimamodel selectionsoftwaretime series

I am attempting to fit a seasonal ARIMA models using ITSM software. The following is the model.

ARIMA$(1,1,0)\times(1,1,0)_{12}$: $\phi(B) \Phi(B^{12}) = (1-\phi B)(1-\Phi B^{12})=1-\Phi B^{12}-\phi_{1}B+\phi_{1} \Phi B^{13}$

Unfortunately, the output from ITSM has a negative coefficient for the 13th term, despite specifying in the software that the 13th term is multiplicative.

Method: Maximum Likelihood
ARMA Model: 
X(t) = - .3485 X(t-1) + .0000 X(t-2) + .0000 X(t-3) + .0000 X(t-4)
       + .0000 X(t-5) + .0000 X(t-6) + .0000 X(t-7) + .0000 X(t-8)
       + .0000 X(t-9) + .0000 X(t-10) + .0000 X(t-11) - .4565 X(t-12)
       - .1591 X(t-13)
     + Z(t)

WN Variance = .001421

AR Coefficients
      -.348543       .000000       .000000       .000000
       .000000       .000000       .000000       .000000
       .000000       .000000       .000000      -.456469
      -.159099

Standard Error of AR Coefficients
       .085875       .000000       .000000       .000000
       .000000       .000000       .000000       .000000
       .000000       .000000       .000000       .082635
       .000000

(Residual SS)/N = .00142069

AICC = -.433387E+03   
BIC  = -.434130E+03   
FPE  = .001469

Any chance that this output is valid in terms of the AICC and the BIC? When I ran an ARIMA$(0,1,1)\times(0,1,1)_{12}$, all of the coefficients where correct in terms of their signs, and the AICC and the BIC was not all that different from this, AICC = -.440119E+03,
BIC = -.444948E+03 .

Here is the data:

112 118 132 129 121 135 148 148 136 119 104 118 115 126 141 135 125 149 170 170 158 133 114 140 145 150 178 163 172 178 199 199 184 162 146 166 171 180 193 181 183 218 230 242 209 191 172 194 196 196 236 235 229 243 264 272 237 211 180 201 204 188 235 227 234 264 302 293 259 229 203 229 242 233 267 269 270 315 364 347 312 274 237 278 284 277 317 313 318 374 413 405 355 306 271 306 315 301 356 348 355 422 465 467 404 347 305 336 340 318 362 348 363 435 491 505 404 359 310 337 360 342 406 396 420 472 548 559 463 407 362 405

Note: I took the log of the data, and used Yule-Walker to get the initial parameter estimates.

Thanks!

Best Answer

I don't have the software ITSM but I could do some inspection upon the data and output that you posted. I fitted the model with different possible mistakes to see which one gives results closer to your output. I used a simplified version of the conditional sum of squares to fit this particular model.

My conclusion is that the code of ITSM that makes the computations is probably correct. The error arises most likely when reporting the results and printing the output. Below I give some details based on some examples that I run in the R software.

I fitted the ARIMA(1,1,0)(1,1,0) model written in expanded form $1 - \phi_{1}B - \Phi B^{12} + \phi_{1} \Phi B^{13}$ by minimizing the squared sum of residuals of the following equation with two parameters ($\phi_1$ and $\Phi$):

$$ y_t = \phi_1 y_{t-1} + \Phi y_{t-12} - \phi_1\Phi y_{t-13} + \hbox{error}_t \,. $$

y <- log(window(AirPassengers, end = c(1959, 12)))
# compute residual sum of squares for ARIMA(1,1,0)(1,1,0) with monthly data
# and minimize using BFGS
# this approach and code is intended just for testing purposes
ssq <- function(coefs, y)
{
  n <- length(y) - 13
  z <- diff(diff(y), 12)
  e <- rep(NA, n-13)
  for (i in seq.int(14, length(z)))
    e[i-13] <- z[i] - coefs[1]*z[i-1] - coefs[2]*z[i-12] + prod(coefs)*z[i-13]
  sum(e^2)
}
# minimize the sum of the squared residuals
res <- optim(ssq, par=c(0,0), method="BFGS", y=y)
# estimated AR coefficients (phi1 and Phi, respectively)
res$par
# [1] -0.3942546 -0.4516274
res$par[1] * res$par[2]
# [1] 0.1780562

The fitted model is therefore (in the expanded form and moving the terms to the right-hand-side):

$$ y_t = 0.3942546 y_{t-1} + 0.4516274 y_{t-12} - 0.1780562 y_{t-13} + \hbox{error}_t \,. $$

As the OP said, the sign of the AR coeffcient related to $B^{13}$ should be the opposite of those for the other coefficients and this is what we get here.

Now, I make intentionally a mistake in the sign of the $B^{13}$ coefficient by typing - prod(coefs)*z[i-13] instead of + prod(coefs)*z[i-13]. Then the following is obtained:

ssq.bug <- function(coefs, y)
{
  n <- length(y) - 13
  z <- diff(diff(y), 12)
  e <- rep(NA, n-13)
  for (i in seq.int(14, length(z)))
    e[i-13] <- z[i] - coefs[1]*z[i-1] - coefs[2]*z[i-12] - prod(coefs)*z[i-13]
  sum(e^2)
}
res.bug <- optim(ssq.bug, par=c(0,0), method="BFGS", y=y)
res.bug$par
# -0.2459668 -0.3511647
res.bug$par[1] * res.bug$par[2]
# [1] 0.08637484

The results are far from the values reported by the OP. The results obtained before are much closer to the output returned by ITSM (except for the sign related to $B^{13}$). Thus, I would say there is no bug in the code that makes the computations.

We could also suspect that the restriction was actually not enforced, but that cannot be the case because .1591 is exactly equal to the product .3485*.4565, so I don't think it is the case either.

It's interesting to see that if the restriction for $B^{13}$ is not enforced and the AR(13) is fitted (with zeros for the AR coefficients of order 2 to 11), then the estimated parameters are very close those obtained above for the restricted model (this may suggest that, given a model of order 13, the multiplicative seasonal ARIMA model is appropriate for the data).

The unrestricted model contains three coefficients, $\phi_1$, $\Phi$ and $\phi_{13}$. The optimization algorithm searches for values of theses coefficients that minimize the sum of squared residuals:

$$ y_t = \phi_1 y_{t-1} + \Phi y_{t-12} - \phi_{13} y_{t-13} + \hbox{error}_t \,. $$

ssq.ur <- function(coefs, y)
{
  n <- length(y) - 13
  z <- diff(diff(y), 12)
  e <- rep(NA, n-13)
  for (i in seq.int(14, length(z)))
    e[i-13] <- z[i] - coefs[1]*z[i-1] - coefs[2]*z[i-12] + coefs[3]*z[i-13]
  sum(e^2)
}
res.ur <- optim(ssq.ur, par=c(0,0,0), method="BFGS", y = y)
res.ur$par
# [1] -0.3951237 -0.4539019  0.1493490
res.ur$par[1] * res.ur$par[2]
# [1] 0.1793474

As we would expect, the parameters estimated in the unrestricted model do not satisfy the restriction exactly, but the last coefficient turns out to be close to the product of the two others. If we ran a test, the restriction would probably not be rejected.

Related Question