There are a number of reasons:
Restricting the search space limits the computation time required. This is an important factor if you are modeling and forecasting many time series.
More complex ARIMA models are hard to interpret. Intepretability may not be high on our list of desirable qualities if all we want to do is forecast, but often we need explain our models to a non-statistician user. I'd rather not try to explain even second differences, let alone third or higher ones.
Forecasters' experience has been that most time series are quite adequately described by ARIMA models of low orders. Let's model the 819 nonseasonal M3 series using auto.arima(,max.p=10,max.q=10,max.d=3):
library(Mcomp)
M3.nonseasonal <- M3[sapply(M3,"[[","period")%in%c("YEARLY","OTHER")]
models <- matrix(NA,nrow=length(M3.nonseasonal),ncol=3,
dimnames=list(names(M3.nonseasonal),c("p","d","q")))
pb <- winProgressBar(max=length(M3.nonseasonal))
for ( ii in seq_along(M3.nonseasonal) ) {
setWinProgressBar(pb,ii,paste(ii,"of",length(M3.nonseasonal)))
fit <- auto.arima(M3.nonseasonal[[ii]]$x,max.p=10,max.q=10,max.d=3)
models[ii,] <- fit$arma[c(1,6,2)]
}
close(pb)
sort(table(
paste0( (",apply(models,1,paste0,collapse=","),")")),decreasing=TRUE)
summary(models)
The output:
> sort(table(paste0("(",apply(models,1,paste0,collapse=","),")")),decreasing=TRUE)
(0,1,0) (0,1,1) (0,2,0) (0,2,1) (1,0,0) (1,1,0) (0,0,0) (0,0,1) (2,0,0) (1,2,0)
413 80 65 60 31 29 23 19 15 13
(1,1,1) (1,2,1) (2,1,0) (2,1,2) (0,3,0) (1,1,2) (2,2,1) (0,1,2) (0,2,2) (1,0,1)
6 6 6 6 5 5 5 4 3 3
(1,3,0) (1,0,2) (1,2,2) (2,0,1) (2,1,1) (2,2,2) (3,0,0) (0,3,1) (1,3,1) (2,2,0)
3 2 2 2 2 2 2 1 1 1
(2,3,1) (3,0,1) (3,1,0) (3,2,3)
1 1 1 1
> summary(models)
p d q
Min. :0.0000 Min. :0.0 Min. :0.0000
1st Qu.:0.0000 1st Qu.:1.0 1st Qu.:0.0000
Median :0.0000 Median :1.0 Median :0.0000
Mean :0.2393 Mean :1.1 Mean :0.2906
3rd Qu.:0.0000 3rd Qu.:1.0 3rd Qu.:1.0000
Max. :3.0000 Max. :3.0 Max. :3.0000
In some few cases, auto.arima() opts for $d=3$, but I'm skeptical about integration of order 3, both because it is very hard to interpret and because it would lead to cubic trends, which rarely make sense. $p$ and $q$ never exceed 3.
Finally, more complex ARIMA models rarely are more accurate. As a matter of fact, the simplest possible ARIMA(0,0,0) model - i.e., white noise, with the optimal forecast being simply the historical mean - often outperforms more complex ARIMA models.
More recently, Petropoulos et al. (2022) found that constraining the "total order" of ARIMA models (i.e., $p+q+P+Q$) to some small number on the order of 4 can lead to more accurate forecasts than letting this go even up to the orders that auto.arima() supports. And at the same time, of course, runtimes decreased drastically compared to the larger possible model space.
As far as I am aware, the specific defaults chosen don't stem from any specific research, and any such research would of course be of questionable generality. If you are really interested in how the package authors came up with these numbers, you could ask them (perhaps report here what you heard back?), but I strongly suspect that the answer will be "experience".
