*! version 1.0.6 07oct2022

Title

arimaauto -- Finds the best [TS] arima model with the
help of a Stata-adjusted Hyndman-Khandakar (2008) algorithm through stepwise
traversing of the model space or a bulk estimation

Requirements

hegy: Hylleberg et al. (1990) seasonal unit-root tests
with OLS and GLS detrending (Rodrigues and Taylor 2007) and lag-length selection
based on MAIC, AIC, BIC, and the sequential method by Tomás del Barrio Castro,
Andrii Bodnar and Andreu Sansó, University of the Balearic Islands Palma, Spain.

kpss: Kwiatkowski et al. (1992) tests for stationarity
of a time series by Christopher F. Baum, Boston College, U.S.

Syntax

arimaauto
[varlist]
[if] [in]
[weight]
[, options]

Best model selection
ic(string) model selection criterion for [TS] arima:
"llf" (maximization of the LLF),
"aic" (minimization of the AIC) or "sic" (minimization
of the SIC) (by default: "aic")

stationary limit the model space to stationary models only,
ARIMA(p,0,q)(P,0,Q) which will disable the automatic
determination of #d with the help of [TS] dfgls
and kpss

noseasonal limit the model space to non-seasonal models only,
ARIMA(p,d,q) which will disable the automatic determination
of #D with the help of hegy

nostepwise switch to bulk estimation from the Hyndman-Khandakar
algorithm

Tests-related
level(#) confidence level, constrained by the critical values of
hegy: 90, 95 or 99
(by default: c(level))

mode(string) lag selection criterion for unit root tests:
"maic", "bic" or "seq" (by default: "maic")

hegy(...) options directly passed to hegy
except gls (enabled by default, to disable specify nogls),
mode(string) and maxlag(#) which will be ignored

dfgls(...) options directly passed
to [TS] dfgls except maxlag(#) which will
be ignored

kpss(...) options directly passed to kpss
except maxlag(#) which will be ignored

ARIMA-related
noconstant suppress constant term if both
arima(#p,#d,#q) and sarima(#P,#D,#Q,#s) are specified

arima(#p,#d,#q) specify ARIMA(p,d,q) model for dependent
variable which will disable the automatic determination of #d
with the help of [TS] dfgls and kpss

sarima(#P,#D,#Q,#s) specify period-#s multiplicative
seasonal ARIMA term which will disable the automatic determination
of #D with the help of hegy

arima(...) options directly passed to
[TS] arima except ar(numlist), ma(numlist),
mar(numlist, #s) and
mma(numlist, #s)} which will be
ignored

Limits and maximum values
max(numlist) limit of p and q
in ARIMA(p,d,q) (by default: 5 5)

mmax(numlist) limit of P and Q
in ARIMA(p,d,q)(P,D,Q) (by default: 2 2)

invroot(real) limit of inverse characteristic roots
in estat aroots (by default: 1/1.001)

maxlag(#) maximum lag order to be used in unit root tests
(by default: . aka unlimited)

maxmodels(#) maximum number of estimated models
(by default: . aka unlimited)

iterate(#) maximum number of iterations(influences convergence)
(by default: 100)

Reporting
trace(#) print hegy, [TS] dfgls and
kpss output if 1 and [TS] arima output
if 2 (by default: 0)

The user must tsset the data before using arimaauto;
see [TS] tsset.

depvar and indepvars may contain time-series operators;
see tsvarlist.

by, collect, fp, rolling, statsby, and xi
are allowed; see prefix.

iweights are allowed; see weights.

See [TS] arima postestimation for features
available after estimation.

Description

arimaauto is de facto an "augmented" Mata-written sister program
to Christopher F. Baum's ARMA-limited arimasel with mutually consistent
output, allowing for ARIMA(p,d,q) and multiplicative seasonal
ARIMA(p,d,q)(P,D,Q) models, selecting the best model based on
the LLF, AIC or SIC, and returning its estimates at the same time. However,
unlike arimasel, the selection is by default performed with the help of
the Hyndman-Khandakar algorithm, first implemented in the
auto.arima
function (part of the "forecast" package) in the R language.

arimaauto can therefore be considered an estimation command, built
on top of [TS] arima, returning the standard arimasel
output in r(...), performed tests in r(tests), and estimated models
in r(models). The user can also access r(ictests), r(icarima),
and hidden r(limits).

Stata-adjusted Hyndman-Khandakar algorithm:

The model selection algorithm described in Hyndman and Khandakar (2008) is
based on a combination of a modified Canova-Hansen seasonal unit root test
(with an empirical formula for calculation of its critical values) and of the
KPSS unit root test, aimed at avoiding (alleged) overdifferencing caused
by tests which assume unit root in their null hypothesis such as hegy
and [TS] dfuller. Since the Canova-Hansen test was unavailable
in Stata 17 and its implementation would have been a feat of its own,
the algorithm was "inverted" to work with more powerful GLS-based
hegy and [TS] dfgls unit root tests with a correction
by the KPSS unit root test to prevent the mentioned overdifferencing aka large
#d in ARIMA(p,d,q) and ARIMA(p,d,q)(P,D,Q) models. The
user can disable GLS in hegy and pass additional options to all the
three tests; see options (syntax).

Unit root tests:

hegy performs the Hylleberg, Engle, Granger, and Yoo (HEGY) (1990) test
for both monthly and quarterly data, allowing for generalized least-squares
(GLS) and ordinary least-squares (OLS) and detrending (controlled by the
nogls option). It automatically determines the lag length (order of
augmentation) based on the modified AIC (MAIC) (del Barrio Castro, Osborn, and
Taylor, 2016), Bayesian information criterion (BIC), and sequential t test
for the last augmentation lag (Hall, 1994; Ng and Perron, 2001). User may
control the deterministic terms with help of the det() option; see
hegy for syntax and examples. The critical values are constructed
from the empirically estimated response surface (del Barrio Castro,
Bodnar, and Sansó, 2015).

NB arimaauto only uses the hegy's t test for Nyquist
frequencies (t_S/2).

[TS] dfgls performs a modified Dickey-Fuller t test for a unit
root in which the series has been transformed by a generalized least-squares
regression. By default, a trend is included; see [TS] dfgls
for syntax and examples.

kpss performs the Kwiatkowski, Phillips, Schmidt, Shin (KPSS, 1992)
test for stationarity of a time series. This test differs from
[TS] dfgls, by having a null hypothesis of stationarity. The
test may be conducted under the null of either trend stationarity
(the default) or level stationarity; see kpss for syntax and
examples.

PS The maximum lag length in the unit root tests is calculated with the
help of the Schwert (1989) formula(s) (the default) or arimaauto's
maxlag(#) option.

Options

For hegy(...) see a list of
hegy options

For dfgls(...) see a list of
dfgls options

For kpss(...) see a list of
kpss options

For arima(...) see a list of
arima options

Methods and formulas

Bulk estimation:

The bulk estimation (activated with the nostepwise option in command
(see options (syntax)) is based on a large
model space generated from combinations of vectors of p, q,
P and Q with values lying in the range <0, limit>. For
example, the default non-seasonal "bulk" model space includes
36 models and the seasonal one 324 models already. Therefore, caution
and the use of maxmodels(#) is advisable.

NB Some models may take a long time to converge or the optimizer may
even become stuck on flat regions with repeated "(backed up)" messages
(if trace(2) was specified). The user is advised to press the
break key in such cases.

PS To match the arimasel command, the user should not forget
to increase the inverse characteristic root limit to 1 with the help of
invroot(1).

Stepwise traversing:

Both the Stata-adjusted and the original Hyndman-Khandakar algorithm consist of
two steps, the second of which is iterated.

Step 1: Four initial models are considered as the model space unless
arima(#p,#d,#q) and/or sarima(#P,#D,#Q,#s) are specified:

• ARIMA(2,d,2) if #s = 0 and
ARIMA(2,d,2)(1,D,1) if #s ≥ 4

• ARIMA(0,d,0) if #s = 0 and
ARIMA(0,d,0)(0,D,0) if #s ≥ 4

• ARIMA(1,d,0) if #s = 0 and
ARIMA(1,d,0)(1,D,0) if #s ≥ 4

• ARIMA(0,d,1) if #s = 0 and
ARIMA(0,d,1)(0,D,1) if #s ≥ 4

Otherwise, the algorithm starts with eventual combinations of the "specified"
and default terms or with a single model. If d + D ≤ 1, the model(s)
is(are) fitted with a constant or else the constant is omitted.

Step 2: Out of the model space, the model with the biggest LLF,
smallest AIC or smallest SIC (based on what is set in
options (syntax)) is selected and is called
the "current" model, of which thirteen variations are considered:

• where one of p, q, P and Q
varies by ±1 from the "current" model;

• where p and q both vary by ±1
from the "current" model;

• where P and Q both vary by ±1
from the "current" model;

• where the constant is excluded/included
if present/absent in the "current" model.

This step is iterated until no better "current" model can be found.

Default limits:

The default limits of the Hyndman-Khandakar algorithm are p ≤ 5,
q ≤ 5, P ≤ 2, Q ≤ 2, every characteristic root ≥ 1.001 (in
absolute value), and an error-free fit of the model, which can be changed with
the help of arimaauto's options (syntax).

Remarks

arimaauto is an interface for the Mata class ARIMAAuto (stored in
the larimaauto.mlib file) which can be a) extended via
inheritance and b) used as a base for third
programs. ARIMAAuto stores several types of variables in an associative
array (of which it is itself a child) declared with:

. mata: AA = ARIMAAuto()

Contents of variables are stored and read with:

. mata: AA.put("key", value)
. mata: AA.get("key")

where "key" is the name of the variable.

The variables are varlist - [varlist] (string scalar), ifin -
concatenated [if] [in] (string scalar), iw - [weight]
(string scalar), o_* - options (string scalar), f_* - flags
(real scalar), level - level (real scalar), mode and ic -
information criteria (string scalar), and finally L, T, and
MS - real matrices of limits, tests, and models.

NB The variables should be declared in the above given order.

L is a (1 x 8) matrix/vector of which r(limits)
is a transposed copy, i.e. r(limits) = AA.get("L")'

T is a (. x 6) matrix of which r(tests) is
a copy, i.e. r(tests) = AA.get("T")

MS is a (. x 12) matrix of which r(models) is
a reduced copy, i.e. r(tests) = AA.get("MS")[,(1,3,4,6,8,10,11,12)]. The
left out columns are #d, #D, #s, and constant (in the given
order).

NB The example of their use can be found in the arimaauto.ado file.

PS ARIMAAuto includes a virtual function get_cv_seas() returning
a real matrix of critical values (test x 0.01,0.05,0.10) for the
hegy test (second row is retrieved), which can be replaced in an
eventual child class when a different test is required (the second row rule must
be observed, for example, coding AA.put("T", J(1,6,.)) in the beginning)
together with the other variables). ARIMAAuto considers unit root to be
present if statistic > critical value, hence signs may need to be adjusted
for unit root tests other than hegy, [TS] dfgls, and
kpss.

PSS All estimations in ARIMAAuto are performed under version 15.

Examples

quarterly data:
. sysuse gnp96.dta, clear
. arimaauto gnp96
. arimaauto gnp96, nostep maxm(15)

monthly data:
. webuse air2.dta, clear
. arimaauto air, sarima(0,1,0,12)

Author

Ilya Bolotov

Prague University of Economics and Business

Prague, Czech Republic

ilya.bolotov@vse.cz

Thanks for citing this software and my works on the topic:

Bolotov, I. (2022). ARIMAAUTO: Stata module to find the best ARIMA model
with the help of a Stata-adjusted Hyndman-Khandakar (2008) algorithm.
Available from https://ideas.repec.org/c/boc/bocode/s459043.html.

References

del Barrio Castro, T., A. Bodnar, and A. Sansó. 2015. Numerical
distribution functions for seasonal unit root tests with OLS and GLS
detrending. Working Paper 73, DEA.
http://dea.uib.es/digitalAssets/353/353054_w73.pdf.

del Barrio Castro, T., D. R. Osborn, and A. M. R. Taylor. 2016. The
performance of lag selection and detrending methods for HEGY seasonal unit
root tests. Econometric Reviews 35: 122-168.

Hall, A. 1994. Testing for a unit root in time series with pretest data-based
model selection. Journal of Business and Economic Statistics 12: 461-470.

Hylleberg, S., R. F. Engle, C. W. J. Granger, and B. S. Yoo. 1990. Seasonal
integration and cointegration. Journal of Econometrics 44: 215-238.

Hyndman, R. J. and Y. Khandakar. 2008. Automatic time series forecasting: the
forecast package for R. Journal of statistical software, 27(1), 1-22.

Kwiatkowski, D., P. C. Phillips, P. Schmidt, and Y. Shin (1992). Testing the
null hypothesis of stationarity against the alternative of a unit root: How
sure are we that economic time series have a unit root?.
Journal of econometrics, 54(1-3), 159-178.

Ng, S., and P. Perron. 2001. Lag length selection and the construction of unit
root tests with good size and power. Econometrica 69: 1519-1554.

Rodrigues, P. M. M., and A. M. R. Taylor. 2007. Efficient tests of the
seasonal unit root hypothesis. Journal of Econometrics 141: 548-573.

Schwert, G. W. 1989. Tests for unit roots: A Monte Carlo investigation.
Journal of Business and Economic Statistics 2: 147-159.