ARIMA Modelling of Simulated Time Series

Box–Jenkins identification, estimation and diagnostic checking

Introduction

This notebook analyses ten simulated time series, Y01–Y10, each of length n = 10{,}000, stored in data_arima.csv. Every series was generated from a non-seasonal ARIMA(p,d,q) model without drift. the analyses are done by using the Box–Jenkins methodology:

1. Identification: inspect the sample ACF and PACF, test for stationarity, difference the series if necessary, and read tentative orders (p,q) from the ACF and PACF of the stationary series.
2. Estimation: fit the tentative model by maximum likelihood and test the significance of every coefficient.
3. Diagnostic checking: verify that the residuals behave like Gaussian white noise.

For every series, the manually identified model is also compared with the model selected automatically by forecast::auto.arima(). All tests use a significance level of \alpha = 5\%.

Identification rules. With n = 10{,}000, the approximate 95% bands of the sample ACF and PACF are \pm 1.96/\sqrt{n} \approx \pm 0.0196 (dashed lines in the plots).

Process	ACF	PACF
AR(p)	tails off	cuts off after lag p
MA(q)	cuts off after lag q	tails off
ARMA(p,q)	tails off	tails off

Libraries

library(tidyverse)
library(forecast)
library(tseries)
library(patchwork)
library(lmtest)

Dataset

The data set contains a time index t and the ten series.

df <- read_csv("data_arima.csv")
head(df)

Each column Y01–Y10 is one series observed at t = 1, \dots, 10{,}000.

df %>%
  pivot_longer(cols = Y01:Y10, names_to = "series", values_to = "value") %>%
  ggplot(aes(x = t, y = value)) +
    geom_line(colour = "steelblue") +
    facet_wrap(~ series, ncol = 2, scales = "free_y") +
    labs(x = "t", y = NULL) +
    theme_bw()

Three series, Y01, Y08 and Y09, wander away from their starting level and never return to a fixed mean. This is the typical appearance of a stochastic trend, i.e. non-stationarity in the mean. The other seven series fluctuate around constant levels (about 20, 75, 30, 100, 10, 60 and 50 for Y02–Y07 and Y10) with roughly constant variance, which suggests stationarity.

ARIMA Modelling

Each subsection follows the same steps: ACF and PACF of the original series, stationarity tests, differencing (when needed), identification, estimation, the fitted model equation, comparison with auto.arima(), and diagnostic checking.

Model 01

ACF & PACF Series 01

We start with the sample ACF of the original series, shown up to lag 5,000 to display its complete decay, and the PACF up to lag 36.

p_acf  <- ggAcf(df$Y01, lag.max = 5000) + ggtitle("ACF Y01") 
p_pacf <- ggPacf(df$Y01, lag.max = 36) + ggtitle("PACF Y01") 

p_acf | p_pacf

The ACF starts near 1 and decays almost linearly: it stays positive and significant for more than 3,000 lags before turning negative. The PACF has a single dominant spike close to 1 at lag 1. Such an extremely slow decay of the ACF is the classic signature of a non-stationary, so the ACF and PACF of the levels cannot be used to identify p and q.

Stationarity Test

Three tests are applied to the original series. The Augmented Dickey–Fuller (ADF) and Phillips–Perron (PP) tests have H_0: the series has a unit root (non-stationary). The KPSS test has the reverse null, H_0: the series is level-stationary.

adf.test(df$Y01)

    Augmented Dickey-Fuller Test

data:  df$Y01
Dickey-Fuller = -2.5227, Lag order = 21, p-value = 0.3569
alternative hypothesis: stationary

ADF: Dickey–Fuller = -2.5227, p-value = 0.3569 > 0.05, so the unit-root hypothesis is not rejected.

pp.test(df$Y01)

    Phillips-Perron Unit Root Test

data:  df$Y01
Dickey-Fuller Z(alpha) = -8.2347, Truncation lag parameter = 12, p-value = 0.6507
alternative hypothesis: stationary

PP: p-value = 0.6507 > 0.05, again no evidence against a unit root.

kpss.test(df$Y01)

    KPSS Test for Level Stationarity

data:  df$Y01
KPSS Level = 67.718, Truncation lag parameter = 12, p-value = 0.01

KPSS: statistic = 67.718, p-value \le 0.01, so level stationarity is rejected. All three tests agree: Y01 is non-stationary in the mean and has to be differenced.

Differencing

We take the first difference W_t = (1-B)Y_t = Y_t - Y_{t-1} and compare the original and the differenced series.

ts_01 <- ts(data=df$Y01)
diff_ts_01 <- diff(ts_01, lag = 1)
autoplot(ts_01) | autoplot(diff_ts_01)

The original series (left) drifts between roughly -100 and 450 without a fixed mean, whereas the differenced series (right) fluctuates around zero with constant variance. The ADF test confirms this.

adf.test(diff_ts_01)

    Augmented Dickey-Fuller Test

data:  diff_ts_01
Dickey-Fuller = -18.95, Lag order = 21, p-value = 0.01
alternative hypothesis: stationary

Dickey–Fuller = -18.95, p-value \le 0.01: the unit root is rejected, so W_t is stationary and d = 1. The ACF and PACF of W_t are now used to identify p and q.

p01_acf_diff  <- ggAcf(diff_ts_01, lag.max = 36) + ggtitle("ACF diff(Y01)") 
p01_pacf_diff <- ggPacf(diff_ts_01, lag.max = 36) + ggtitle("PACF diff(Y01)") 

p01_acf_diff | p01_pacf_diff

The ACF of W_t has significant spikes at lags 1 to 8 and then cuts off abruptly: from lag 9 onward all autocorrelations lie inside the bands. The PACF, in contrast, dies down gradually in an irregular, oscillating way. A cut-off ACF combined with a tailing-off PACF points to a pure MA process, and the last significant ACF lag gives q = 8. The tentative model is therefore ARIMA(0,1,8), without a constant because the differenced series has mean zero (no drift).

Parameter Estimation

ARIMA(0,1,8) is estimated by maximum likelihood with Arima().

m01 <- Arima(ts_01, order = c(0,1,8))
summary(m01)

Series: ts_01 
ARIMA(0,1,8) 

Coefficients:
         ma1     ma2     ma3      ma4     ma5     ma6     ma7     ma8
      0.4000  0.2042  0.4099  -0.5064  0.3076  0.3090  0.2052  0.4695
s.e.  0.0088  0.0095  0.0093   0.0098  0.0096  0.0091  0.0094  0.0087

sigma^2 = 0.9913:  log likelihood = -14142.54
AIC=28303.08   AICc=28303.09   BIC=28367.97

Training set error measures:
                      ME     RMSE       MAE      MPE     MAPE      MASE         ACF1
Training set -0.01413526 0.995185 0.7960329 9.455486 15.21431 0.6921013 -0.003393213

All eight MA coefficients are estimated precisely (standard errors below 0.01). The residual variance is \hat\sigma_a^2 = 0.9913, with AIC = 28303.08, AICc = 28303.09 and BIC = 28367.97. Next, each coefficient is tested with H_0: \theta_j = 0.

coeftest(m01)

z test of coefficients:

      Estimate Std. Error z value  Pr(>|z|)    
ma1  0.4000239  0.0088346  45.279 < 2.2e-16 ***
ma2  0.2041584  0.0094758  21.545 < 2.2e-16 ***
ma3  0.4098869  0.0092919  44.112 < 2.2e-16 ***
ma4 -0.5064212  0.0097596 -51.890 < 2.2e-16 ***
ma5  0.3075549  0.0095837  32.092 < 2.2e-16 ***
ma6  0.3090082  0.0091218  33.876 < 2.2e-16 ***
ma7  0.2051990  0.0094095  21.808 < 2.2e-16 ***
ma8  0.4694732  0.0087348  53.748 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Every coefficient is highly significant (|z| \ge 21.5, p-value < 2.2\times10^{-16}), so no term can be dropped.

Fitted Model Equation

The general form of ARIMA(0,1,8) is

(1-B)\,Y_t = (1+\theta_1B+\theta_2B^2+\cdots+\theta_8B^8)\,a_t .

Substituting the estimates gives

\begin{aligned} (1-B)\,Y_t ={}& (1 + 0.4000B + 0.2042B^2 + 0.4099B^3 - 0.5064B^4 \\ &+ 0.3076B^5 + 0.3090B^6 + 0.2052B^7 + 0.4695B^8)\,a_t , \end{aligned}

and, written explicitly in terms of past observations and shocks,

\begin{aligned} Y_t ={}& Y_{t-1} + a_t + 0.4000\,a_{t-1} + 0.2042\,a_{t-2} + 0.4099\,a_{t-3} - 0.5064\,a_{t-4} \\ &+ 0.3076\,a_{t-5} + 0.3090\,a_{t-6} + 0.2052\,a_{t-7} + 0.4695\,a_{t-8}, \qquad a_t \sim \mathcal{N}(0,\ 0.9913). \end{aligned}

Each new value equals the previous value plus a shock that is a weighted sum of the current and the last eight innovations.

Comparison with auto.arima()

As a benchmark, auto.arima() selects a model automatically, allowing up to ten AR and ten MA terms.

aa01 <- auto.arima(ts_01, max.p = 10, max.q = 10)
aa01

Series: ts_01 
ARIMA(9,1,6) 

Coefficients:
         ar1     ar2      ar3      ar4     ar5     ar6      ar7      ar8     ar9     ma1      ma2     ma3      ma4     ma5     ma6
      0.2759  0.7792  -0.3439  -0.5522  0.3629  0.2653  -0.2345  -0.1193  0.1334  0.1203  -0.6988  0.3862  -0.0672  0.0550  0.4472
s.e.  0.0354  0.0287   0.0441   0.0260  0.0358  0.0252   0.0282   0.0208  0.0198  0.0347   0.0247  0.0437   0.0212  0.0201  0.0175

sigma^2 = 1.003:  log likelihood = -14195.62
AIC=28423.23   AICc=28423.29   BIC=28538.6

auto.arima() selects ARIMA(9,1,6), a model with 15 coefficients. It finds the correct order of differencing (d = 1), but its AICc (28423.29) is 120.2 units larger than the AICc of the manually identified ARIMA(0,1,8) (28303.09), and its BIC is 170.6 units larger. The automatic procedure therefore did not even find the model that is best by its own criterion: for a long series it searches stepwise using an approximate likelihood, and here it stopped in an over-parameterised local optimum.

coeftest(aa01)

z test of coefficients:

     Estimate Std. Error  z value  Pr(>|z|)    
ar1  0.275936   0.035350   7.8057 5.916e-15 ***
ar2  0.779181   0.028739  27.1124 < 2.2e-16 ***
ar3 -0.343920   0.044075  -7.8031 6.040e-15 ***
ar4 -0.552176   0.026017 -21.2236 < 2.2e-16 ***
ar5  0.362873   0.035785  10.1404 < 2.2e-16 ***
ar6  0.265342   0.025232  10.5162 < 2.2e-16 ***
ar7 -0.234477   0.028198  -8.3155 < 2.2e-16 ***
ar8 -0.119303   0.020786  -5.7397 9.483e-09 ***
ar9  0.133430   0.019826   6.7299 1.698e-11 ***
ma1  0.120336   0.034652   3.4727 0.0005152 ***
ma2 -0.698816   0.024715 -28.2750 < 2.2e-16 ***
ma3  0.386185   0.043749   8.8273 < 2.2e-16 ***
ma4 -0.067231   0.021216  -3.1688 0.0015306 ** 
ma5  0.054953   0.020058   2.7397 0.0061494 ** 
ma6  0.447165   0.017529  25.5104 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

All 15 coefficients are individually significant, but significance alone does not make this a better model: it is far less parsimonious and fits worse than ARIMA(0,1,8) according to both AICc and BIC. The manual model is preferred.

Diagnostic Checking

An adequate model leaves residuals that behave like white noise. checkresiduals() plots the residuals, their ACF and their histogram, and runs the Ljung–Box test with H_0: no autocorrelation up to the chosen lag (here 11 lags, of which 8 degrees of freedom are absorbed by the MA coefficients).

checkresiduals(m01)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,8)
Q* = 1.7283, df = 3, p-value = 0.6307

Model df: 8.   Total lags used: 11

Q^* = 1.7283 with 3 degrees of freedom and p-value = 0.6307 > 0.05: the residuals show no remaining autocorrelation. The residual ACF stays within the bands and the histogram is bell-shaped. Normality is checked formally with the Jarque–Bera test (H_0: the residuals are normally distributed).

res_01 <- m01$residuals
jarque.bera.test(res_01)

    Jarque Bera Test

data:  res_01
X-squared = 0.70053, df = 2, p-value = 0.7045

X^2 = 0.70053, p-value = 0.7045 > 0.05: normality is not rejected.

Conclusion. ARIMA(0,1,8) is adequate for Y01: all coefficients are significant and the residuals behave like Gaussian white noise.

Model 02

ACF & PACF Series 02

Y02 looks stationary in the overview plot, so its ACF and PACF are inspected directly (up to lag 36).

p_acf  <- ggAcf(df$Y02, lag.max = 36) + ggtitle("ACF Y02") 
p_pacf <- ggPacf(df$Y02, lag.max = 36) + ggtitle("PACF Y02") 

p_acf | p_pacf

The ACF decays exponentially from about 0.71 at lag 1 and becomes negligible after roughly lag 10, while the PACF has one large spike at lag 1 (about 0.71) and is insignificant afterwards. A tailing-off ACF together with a PACF that cuts off after lag 1 is the pattern of an AR(1) process, provided the series is stationary.

Stationarity Test

The same three tests are applied (ADF and PP: H_0 unit root; KPSS: H_0 level-stationary).

adf.test(df$Y02)

    Augmented Dickey-Fuller Test

data:  df$Y02
Dickey-Fuller = -20.106, Lag order = 21, p-value = 0.01
alternative hypothesis: stationary

ADF: Dickey–Fuller = -20.106, p-value \le 0.01, so the unit root is rejected.

pp.test(df$Y02)

    Phillips-Perron Unit Root Test

data:  df$Y02
Dickey-Fuller Z(alpha) = -2930.6, Truncation lag parameter = 12, p-value = 0.01
alternative hypothesis: stationary

PP: Z(\alpha) = -2930.6, p-value \le 0.01, so the unit root is rejected again.

kpss.test(df$Y02)

    KPSS Test for Level Stationarity

data:  df$Y02
KPSS Level = 0.19092, Truncation lag parameter = 12, p-value = 0.1

KPSS: statistic = 0.19092, p-value \ge 0.10, so level stationarity is not rejected. All tests agree that Y02 is stationary, so d = 0 and no differencing is needed. Because the series fluctuates around about 20, the tentative model is ARIMA(1,0,0) with a non-zero mean.

Parameter Estimation

ts_02 <- ts(data=df$Y02)
m02 <- Arima(ts_02, order = c(1,0,0))
summary(m02)

Series: ts_02 
ARIMA(1,0,0) with non-zero mean 

Coefficients:
         ar1    mean
      0.7076  19.957
s.e.  0.0071   0.034

sigma^2 = 0.9888:  log likelihood = -14132.25
AIC=28270.5   AICc=28270.5   BIC=28292.13

Training set error measures:
                       ME      RMSE       MAE        MPE     MAPE      MASE        ACF1
Training set 1.938135e-05 0.9942682 0.7979377 -0.2518663 4.028156 0.9253059 0.003414323

The estimates are \hat\phi_1 = 0.7076 and \hat\mu = 19.957. The value labelled mean (intercept in coeftest()) is the process mean, not the regression constant. The residual variance is \hat\sigma_a^2 = 0.9888, with AIC = AICc = 28270.5 and BIC = 28292.13.

coeftest(m02)

z test of coefficients:

            Estimate Std. Error z value  Pr(>|z|)    
ar1        0.7075777  0.0070623  100.19 < 2.2e-16 ***
intercept 19.9570229  0.0339929  587.09 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Both parameters are highly significant (z = 100.19 for \phi_1 and z = 587.09 for \mu).

Fitted Model Equation

The general form of ARIMA(1,0,0) with mean \mu is

(1-\phi_1B)\,(Y_t-\mu) = a_t .

With the estimates,

(1 - 0.7076B)\,(Y_t - 19.9570) = a_t ,

which expands to

Y_t = 5.8359 + 0.7076\,Y_{t-1} + a_t,\qquad a_t \sim \mathcal{N}(0,\ 0.9888),

where 5.8359 = 19.9570 \times (1 - 0.7076). Since |\hat\phi_1| < 1, the fitted model is stationary: the effect of a shock shrinks by a factor of 0.7076 every period.

Comparison with auto.arima()

aa02 <- auto.arima(ts_02, max.p = 10, max.q = 10)
aa02

Series: ts_02 
ARIMA(1,0,0) with non-zero mean 

Coefficients:
         ar1    mean
      0.7076  19.957
s.e.  0.0071   0.034

sigma^2 = 0.9888:  log likelihood = -14132.25
AIC=28270.5   AICc=28270.5   BIC=28292.13

auto.arima() selects exactly the same model, ARIMA(1,0,0) with a non-zero mean, with identical estimates and information criteria. Manual and automatic identification agree.

Diagnostic Checking

checkresiduals(m02)

    Ljung-Box test

data:  Residuals from ARIMA(1,0,0) with non-zero mean
Q* = 3.5031, df = 9, p-value = 0.941

Model df: 1.   Total lags used: 10

Ljung–Box: Q^* = 3.5031 with 9 degrees of freedom, p-value = 0.941: the residuals are white noise. The residual ACF lies within the bands and the histogram is symmetric and bell-shaped.

res_02 <- m02$residuals
jarque.bera.test(res_02)

    Jarque Bera Test

data:  res_02
X-squared = 1.9585, df = 2, p-value = 0.3756

Jarque–Bera: X^2 = 1.9585, p-value = 0.3756: normality is not rejected.

Conclusion. ARIMA(1,0,0) with a non-zero mean is adequate for Y02.

Model 03

ACF & PACF Series 03

p_acf  <- ggAcf(df$Y03, lag.max = 36) + ggtitle("ACF Y03") 
p_pacf <- ggPacf(df$Y03, lag.max = 36) + ggtitle("PACF Y03") 

p_acf | p_pacf

The ACF has significant spikes at lags 1 to 8 (positive at lags 1, 7 and 8, negative at lags 2 to 6) and cuts off sharply after lag 8. The PACF decays slowly in an irregular, oscillating way that is still visible around lag 20. This is the pattern of an MA(8) process, provided the series is stationary.

Stationarity Test

adf.test(df$Y03)

    Augmented Dickey-Fuller Test

data:  df$Y03
Dickey-Fuller = -26.085, Lag order = 21, p-value = 0.01
alternative hypothesis: stationary

ADF: Dickey–Fuller = -26.085, p-value \le 0.01, so the unit root is rejected.

pp.test(df$Y03)

    Phillips-Perron Unit Root Test

data:  df$Y03
Dickey-Fuller Z(alpha) = -3190.9, Truncation lag parameter = 12, p-value = 0.01
alternative hypothesis: stationary

PP: Z(\alpha) = -3190.9, p-value \le 0.01, so the unit root is rejected.

kpss.test(df$Y03)

    KPSS Test for Level Stationarity

data:  df$Y03
KPSS Level = 0.078839, Truncation lag parameter = 12, p-value = 0.1

KPSS: statistic = 0.078839, p-value \ge 0.10, so level stationarity is not rejected. Y03 is stationary (d = 0), and the tentative model is ARIMA(0,0,8) with a non-zero mean.

Parameter Estimation

ts_03 <- ts(data=df$Y03)
m03 <- Arima(ts_03, order = c(0,0,8))
summary(m03)

Series: ts_03 
ARIMA(0,0,8) with non-zero mean 

Coefficients:
         ma1      ma2      ma3      ma4      ma5      ma6     ma7     ma8     mean
      0.5026  -0.3882  -0.3044  -0.3055  -0.2031  -0.4057  0.2022  0.3961  74.9950
s.e.  0.0092   0.0103   0.0103   0.0105   0.0102   0.0105  0.0102  0.0092   0.0098

sigma^2 = 3.934:  log likelihood = -21034.82
AIC=42089.64   AICc=42089.66   BIC=42161.74

Training set error measures:
                        ME     RMSE      MAE         MPE     MAPE      MASE        ACF1
Training set -0.0001267169 1.982482 1.582264 -0.07077227 2.114717 0.6129917 0.002989087

The eight MA coefficients and the mean (\hat\mu = 74.9950) are estimated precisely. The residual variance is \hat\sigma_a^2 = 3.934, with AIC = 42089.64, AICc = 42089.66 and BIC = 42161.74.

coeftest(m03)

z test of coefficients:

            Estimate Std. Error  z value  Pr(>|z|)    
ma1        0.5026460  0.0092477   54.354 < 2.2e-16 ***
ma2       -0.3882372  0.0102587  -37.845 < 2.2e-16 ***
ma3       -0.3043659  0.0103447  -29.422 < 2.2e-16 ***
ma4       -0.3054866  0.0104822  -29.143 < 2.2e-16 ***
ma5       -0.2031159  0.0101965  -19.920 < 2.2e-16 ***
ma6       -0.4056971  0.0104551  -38.804 < 2.2e-16 ***
ma7        0.2021850  0.0101626   19.895 < 2.2e-16 ***
ma8        0.3961114  0.0091550   43.267 < 2.2e-16 ***
intercept 74.9949528  0.0097968 7655.048 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

All coefficients are highly significant (|z| \ge 19.9).

Fitted Model Equation

The general form of ARIMA(0,0,8) with mean \mu is

Y_t-\mu = (1+\theta_1B+\theta_2B^2+\cdots+\theta_8B^8)\,a_t .

With the estimates,

\begin{aligned} Y_t - 74.9950 ={}& (1 + 0.5026B - 0.3882B^2 - 0.3044B^3 - 0.3055B^4 \\ &- 0.2031B^5 - 0.4057B^6 + 0.2022B^7 + 0.3961B^8)\,a_t , \end{aligned}

that is, explicitly,

\begin{aligned} Y_t ={}& 74.9950 + a_t + 0.5026\,a_{t-1} - 0.3882\,a_{t-2} - 0.3044\,a_{t-3} - 0.3055\,a_{t-4} \\ &- 0.2031\,a_{t-5} - 0.4057\,a_{t-6} + 0.2022\,a_{t-7} + 0.3961\,a_{t-8}, \qquad a_t \sim \mathcal{N}(0,\ 3.934). \end{aligned}

For a pure MA model the constant equals the mean, and observations more than eight periods apart are uncorrelated.

Comparison with auto.arima()

aa03 <- auto.arima(ts_03, max.p = 10, max.q = 10)
aa03

Series: ts_03 
ARIMA(10,0,2) with non-zero mean 

Coefficients:
         ar1     ar2      ar3      ar4      ar5      ar6     ar7      ar8      ar9    ar10     ma1      ma2     mean
      0.2448  0.1229  -0.2527  -0.0840  -0.1607  -0.2612  0.2524  -0.0301  -0.1835  0.0105  0.2477  -0.6464  74.9947
s.e.  0.0356  0.0193   0.0126   0.0179   0.0147   0.0159  0.0195   0.0115   0.0116  0.0146  0.0341   0.0328   0.0091

sigma^2 = 4.097:  log likelihood = -21235.51
AIC=42499.02   AICc=42499.06   BIC=42599.96

auto.arima() selects ARIMA(10,0,2) with a non-zero mean. This model fits clearly worse: \hat\sigma_a^2 = 4.097 instead of 3.934, AICc = 42499.06 (409.4 units larger than for MA(8)) and BIC = 42599.96 (438.2 units larger). An MA(8) process can be written exactly only as an AR model of infinite order; with at most ten AR terms plus two MA terms, the automatic model is just a rough approximation of it, and the stepwise search never reached the pure MA(8) structure. The manual MA(8) is far better.

Diagnostic Checking

checkresiduals(m03)

    Ljung-Box test

data:  Residuals from ARIMA(0,0,8) with non-zero mean
Q* = 3.1399, df = 3, p-value = 0.3706

Model df: 8.   Total lags used: 11

Ljung–Box: Q^* = 3.1399 with 3 degrees of freedom, p-value = 0.3706: no remaining autocorrelation.

res_03 <- m03$residuals
jarque.bera.test(res_03)

    Jarque Bera Test

data:  res_03
X-squared = 0.22751, df = 2, p-value = 0.8925

Jarque–Bera: X^2 = 0.22751, p-value = 0.8925: normality is not rejected.

Conclusion. ARIMA(0,0,8) with a non-zero mean is adequate for Y03.

Model 04

ACF & PACF Series 04

p_acf  <- ggAcf(df$Y04, lag.max = 36) + ggtitle("ACF Y04") 
p_pacf <- ggPacf(df$Y04, lag.max = 36) + ggtitle("PACF Y04") 

p_acf | p_pacf

Only the first two autocorrelations are significant (about 0.66 at lag 1 and 0.18 at lag 2), so the ACF cuts off after lag 2. The PACF alternates in sign (positive, negative, positive, …) and dies down within about eight lags. This is the pattern of an MA(2) process.

Stationarity Test

adf.test(df$Y04)

    Augmented Dickey-Fuller Test

data:  df$Y04
Dickey-Fuller = -20.116, Lag order = 21, p-value = 0.01
alternative hypothesis: stationary

ADF: Dickey–Fuller = -20.116, p-value \le 0.01, so the unit root is rejected.

pp.test(df$Y04)

    Phillips-Perron Unit Root Test

data:  df$Y04
Dickey-Fuller Z(alpha) = -2502.1, Truncation lag parameter = 12, p-value = 0.01
alternative hypothesis: stationary

PP: Z(\alpha) = -2502.1, p-value \le 0.01, so the unit root is rejected.

kpss.test(df$Y04)

    KPSS Test for Level Stationarity

data:  df$Y04
KPSS Level = 0.2533, Truncation lag parameter = 12, p-value = 0.1

KPSS: statistic = 0.2533, p-value \ge 0.10, so level stationarity is not rejected. Y04 is stationary (d = 0), and the tentative model is ARIMA(0,0,2) with a non-zero mean.

Parameter Estimation

ts_04 <- ts(data=df$Y04)
m04 <- Arima(ts_04, order = c(0,0,2))
summary(m04)

Series: ts_04 
ARIMA(0,0,2) with non-zero mean 

Coefficients:
         ma1     ma2     mean
      1.2064  0.5018  30.0235
s.e.  0.0086  0.0087   0.0407

sigma^2 = 2.259:  log likelihood = -18263.33
AIC=36534.67   AICc=36534.67   BIC=36563.51

Training set error measures:
                       ME     RMSE      MAE        MPE     MAPE      MASE         ACF1
Training set 6.009103e-05 1.502776 1.202647 -0.2550741 4.050067 0.7452819 -0.003071841

The estimates are \hat\theta_1 = 1.2064, \hat\theta_2 = 0.5018 and \hat\mu = 30.0235, with \hat\sigma_a^2 = 2.259, AIC = AICc = 36534.67 and BIC = 36563.51.

coeftest(m04)

z test of coefficients:

            Estimate Std. Error z value  Pr(>|z|)    
ma1        1.2063581  0.0086365 139.682 < 2.2e-16 ***
ma2        0.5017600  0.0087254  57.505 < 2.2e-16 ***
intercept 30.0235158  0.0406936 737.794 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

All parameters are highly significant (z = 139.68, 57.51 and 737.79).

Fitted Model Equation

The general form of ARIMA(0,0,2) with mean \mu is

Y_t-\mu = (1+\theta_1B+\theta_2B^2)\,a_t .

With the estimates,

Y_t - 30.0235 = (1 + 1.2064B + 0.5018B^2)\,a_t ,

that is, explicitly,

Y_t = 30.0235 + a_t + 1.2064\,a_{t-1} + 0.5018\,a_{t-2},\qquad a_t \sim \mathcal{N}(0,\ 2.259).

Comparison with auto.arima()

aa04 <- auto.arima(ts_04, max.p = 10, max.q = 10)
aa04

Series: ts_04 
ARIMA(0,0,2) with non-zero mean 

Coefficients:
         ma1     ma2     mean
      1.2064  0.5018  30.0235
s.e.  0.0086  0.0087   0.0407

sigma^2 = 2.259:  log likelihood = -18263.33
AIC=36534.67   AICc=36534.67   BIC=36563.51

auto.arima() selects the same model, ARIMA(0,0,2) with a non-zero mean, with identical estimates. Both approaches agree.

Diagnostic Checking

checkresiduals(m04)

    Ljung-Box test

data:  Residuals from ARIMA(0,0,2) with non-zero mean
Q* = 2.9641, df = 8, p-value = 0.9366

Model df: 2.   Total lags used: 10

Ljung–Box: Q^* = 2.9641 with 8 degrees of freedom, p-value = 0.9366: the residuals are white noise.

res_04 <- m04$residuals
jarque.bera.test(res_04)

    Jarque Bera Test

data:  res_04
X-squared = 5.9529, df = 2, p-value = 0.05097

Jarque–Bera: X^2 = 5.9529, p-value = 0.05097, so normality is (narrowly) not rejected at the 5% level. With n = 10{,}000 the Jarque–Bera test reacts to very small departures; the histogram is bell-shaped, so the residuals can be treated as normal.

Conclusion. ARIMA(0,0,2) with a non-zero mean is adequate for Y04.

Model 05

ACF & PACF Series 05

p_acf  <- ggAcf(df$Y05, lag.max = 36) + ggtitle("ACF Y05") 
p_pacf <- ggPacf(df$Y05, lag.max = 36) + ggtitle("PACF Y05") 

p_acf | p_pacf

The ACF decays slowly in a wave-like pattern: large positive values up to about lag 8 are followed by small negative values between lags 10 and 20. The PACF has six significant spikes (lags 1 to 6, the last two negative) and cuts off after lag 6. This is the pattern of an AR(6) process.

Stationarity Test

adf.test(df$Y05)

    Augmented Dickey-Fuller Test

data:  df$Y05
Dickey-Fuller = -20.009, Lag order = 21, p-value = 0.01
alternative hypothesis: stationary

ADF: Dickey–Fuller = -20.009, p-value \le 0.01, so the unit root is rejected.

pp.test(df$Y05)

    Phillips-Perron Unit Root Test

data:  df$Y05
Dickey-Fuller Z(alpha) = -9051, Truncation lag parameter = 12, p-value = 0.01
alternative hypothesis: stationary

PP: Z(\alpha) = -9051, p-value \le 0.01, so the unit root is rejected.

kpss.test(df$Y05)

    KPSS Test for Level Stationarity

data:  df$Y05
KPSS Level = 0.18261, Truncation lag parameter = 12, p-value = 0.1

KPSS: statistic = 0.18261, p-value \ge 0.10, so level stationarity is not rejected. Y05 is stationary (d = 0), and the tentative model is ARIMA(6,0,0) with a non-zero mean.

Parameter Estimation

ts_05 <- ts(data=df$Y05)
m05 <- Arima(ts_05, order = c(6,0,0))
summary(m05)

Series: ts_05 
ARIMA(6,0,0) with non-zero mean 

Coefficients:
         ar1     ar2     ar3     ar4      ar5      ar6     mean
      0.0966  0.3931  0.3987  0.3931  -0.2997  -0.2931  99.9288
s.e.  0.0096  0.0091  0.0091  0.0091   0.0091   0.0096   0.0967

sigma^2 = 9.068:  log likelihood = -25210.7
AIC=50437.4   AICc=50437.41   BIC=50495.08

Training set error measures:
                       ME     RMSE      MAE        MPE     MAPE      MASE        ACF1
Training set 8.236513e-06 3.010228 2.408564 -0.0912889 2.417534 0.6537646 0.005604284

The six AR coefficients and the mean (\hat\mu = 99.9288) are estimated with standard errors below 0.01 (0.0967 for the mean). The residual variance is \hat\sigma_a^2 = 9.068, with AIC = 50437.40, AICc = 50437.41 and BIC = 50495.08.

coeftest(m05)

z test of coefficients:

            Estimate Std. Error  z value  Pr(>|z|)    
ar1        0.0966306  0.0095593   10.109 < 2.2e-16 ***
ar2        0.3931050  0.0091289   43.062 < 2.2e-16 ***
ar3        0.3986794  0.0091321   43.657 < 2.2e-16 ***
ar4        0.3931380  0.0091266   43.076 < 2.2e-16 ***
ar5       -0.2996752  0.0091363  -32.800 < 2.2e-16 ***
ar6       -0.2930599  0.0095665  -30.634 < 2.2e-16 ***
intercept 99.9287715  0.0967230 1033.144 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

All coefficients are significant; the smallest test statistic is z = 10.11 for \phi_1.

Fitted Model Equation

The general form of ARIMA(6,0,0) with mean \mu is

(1-\phi_1B-\phi_2B^2-\cdots-\phi_6B^6)\,(Y_t-\mu) = a_t .

With the estimates,

(1 - 0.0966B - 0.3931B^2 - 0.3987B^3 - 0.3931B^4 + 0.2997B^5 + 0.2931B^6)\,(Y_t - 99.9288) = a_t ,

which expands to

\begin{aligned} Y_t ={}& 31.0960 + 0.0966\,Y_{t-1} + 0.3931\,Y_{t-2} + 0.3987\,Y_{t-3} + 0.3931\,Y_{t-4} \\ &- 0.2997\,Y_{t-5} - 0.2931\,Y_{t-6} + a_t, \qquad a_t \sim \mathcal{N}(0,\ 9.068), \end{aligned}

where 31.0960 = 99.9288 \times (1 - 0.6888) and 0.6888 is the sum of the six AR coefficients.

Comparison with auto.arima()

aa05 <- auto.arima(ts_05, max.p = 10, max.q = 10)
aa05

Series: ts_05 
ARIMA(7,0,0) with non-zero mean 

Coefficients:
         ar1     ar2     ar3     ar4      ar5      ar6     ar7     mean
      0.1022  0.3986  0.3911  0.3859  -0.3070  -0.2947  0.0185  99.9288
s.e.  0.0100  0.0096  0.0099  0.0100   0.0099   0.0096  0.0100   0.0985

sigma^2 = 9.066:  log likelihood = -25208.99
AIC=50435.98   AICc=50436   BIC=50500.88

auto.arima() selects ARIMA(7,0,0), i.e. one more AR term. The extra coefficient is tiny (\hat\phi_7 = 0.0185 with standard error 0.0100), and the gain is negligible: AICc is only 1.41 units lower (50436.00 vs 50437.41), while BIC is 5.80 units higher (50500.88 vs 50495.08).

coeftest(aa05)

z test of coefficients:

            Estimate Std. Error   z value Pr(>|z|)    
ar1        0.1022347  0.0099994   10.2241  < 2e-16 ***
ar2        0.3985785  0.0096067   41.4897  < 2e-16 ***
ar3        0.3911033  0.0099422   39.3375  < 2e-16 ***
ar4        0.3858599  0.0099606   38.7386  < 2e-16 ***
ar5       -0.3069778  0.0099450  -30.8677  < 2e-16 ***
ar6       -0.2946904  0.0096129  -30.6556  < 2e-16 ***
ar7        0.0184744  0.0100009    1.8473  0.06471 .  
intercept 99.9287867  0.0985301 1014.1960  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The z-test confirms that ar7 is not significant at the 5% level (z \approx 1.85). An AICc difference below 2 means the two models are practically equivalent, and BIC, which penalises each parameter by \ln n \approx 9.21 instead of 2, prefers AR(6). By the principle of parsimony, the manual AR(6) is preferred.

Diagnostic Checking

checkresiduals(m05)

    Ljung-Box test

data:  Residuals from ARIMA(6,0,0) with non-zero mean
Q* = 4.0466, df = 4, p-value = 0.3997

Model df: 6.   Total lags used: 10

Ljung–Box: Q^* = 4.0466 with 4 degrees of freedom, p-value = 0.3997: no remaining autocorrelation.

res_05 <- m05$residuals
jarque.bera.test(res_05)

    Jarque Bera Test

data:  res_05
X-squared = 0.62972, df = 2, p-value = 0.7299

Jarque–Bera: X^2 = 0.62972, p-value = 0.7299: normality is not rejected.

Conclusion. ARIMA(6,0,0) with a non-zero mean is adequate for Y05.

Model 06

ACF & PACF Series 06

p_acf  <- ggAcf(df$Y06, lag.max = 36) + ggtitle("ACF Y06") 
p_pacf <- ggPacf(df$Y06, lag.max = 36) + ggtitle("PACF Y06") 

p_acf | p_pacf

The ACF has a single significant, negative spike at lag 1 (about -0.43) and is insignificant afterwards. The PACF is negative at the first lags and decays exponentially towards zero. This is the pattern of an MA(1) process with a negative coefficient.

Stationarity Test

adf.test(df$Y06)

    Augmented Dickey-Fuller Test

data:  df$Y06
Dickey-Fuller = -22.722, Lag order = 21, p-value = 0.01
alternative hypothesis: stationary

ADF: Dickey–Fuller = -22.722, p-value \le 0.01, so the unit root is rejected.

pp.test(df$Y06)

    Phillips-Perron Unit Root Test

data:  df$Y06
Dickey-Fuller Z(alpha) = -11880, Truncation lag parameter = 12, p-value = 0.01
alternative hypothesis: stationary

PP: Z(\alpha) = -11880, p-value \le 0.01, so the unit root is rejected.

kpss.test(df$Y06)

    KPSS Test for Level Stationarity

data:  df$Y06
KPSS Level = 0.046055, Truncation lag parameter = 12, p-value = 0.1

KPSS: statistic = 0.046055, p-value \ge 0.10, so level stationarity is not rejected. Y06 is stationary (d = 0), and the tentative model is ARIMA(0,0,1) with a non-zero mean.

Parameter Estimation

ts_06 <- ts(data=df$Y06)
m06 <- Arima(ts_06, order = c(0,0,1))
summary(m06)

Series: ts_06 
ARIMA(0,0,1) with non-zero mean 

Coefficients:
          ma1    mean
      -0.5936  9.9999
s.e.   0.0081  0.0041

sigma^2 = 1.023:  log likelihood = -14303.21
AIC=28612.43   AICc=28612.43   BIC=28634.06

Training set error measures:
                       ME     RMSE       MAE       MPE     MAPE      MASE        ACF1
Training set 0.0001789569 1.011426 0.8086233 -1.056863 8.270721 0.5086325 0.004490465

The estimates are \hat\theta_1 = -0.5936 and \hat\mu = 9.9999, with \hat\sigma_a^2 = 1.023, AIC = AICc = 28612.43 and BIC = 28634.06. The implied lag-1 autocorrelation, \hat\theta_1/(1+\hat\theta_1^2) = -0.439, matches the spike seen in the ACF.

coeftest(m06)

z test of coefficients:

            Estimate Std. Error  z value  Pr(>|z|)    
ma1       -0.5936216  0.0081425  -72.904 < 2.2e-16 ***
intercept  9.9999364  0.0041108 2432.584 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Both parameters are highly significant (z = -72.90 and z = 2432.58).

Fitted Model Equation

The general form of ARIMA(0,0,1) with mean \mu is

Y_t-\mu = (1+\theta_1B)\,a_t .

With the estimates,

Y_t - 9.9999 = (1 - 0.5936B)\,a_t ,

that is, explicitly,

Y_t = 9.9999 + a_t - 0.5936\,a_{t-1},\qquad a_t \sim \mathcal{N}(0,\ 1.023).

Since |\hat\theta_1| < 1, the model is invertible.

Comparison with auto.arima()

aa06 <- auto.arima(ts_06, max.p = 10, max.q = 10)
aa06

Series: ts_06 
ARIMA(1,0,1) with non-zero mean 

Coefficients:
         ar1      ma1    mean
      0.0125  -0.6018  9.9999
s.e.  0.0170   0.0136  0.0041

sigma^2 = 1.023:  log likelihood = -14302.94
AIC=28613.88   AICc=28613.88   BIC=28642.72

auto.arima() selects ARIMA(1,0,1) with a non-zero mean, adding an AR term with \hat\phi_1 = 0.0125 (standard error 0.0170).

coeftest(aa06)

z test of coefficients:

           Estimate Std. Error   z value Pr(>|z|)    
ar1        0.012518   0.016968    0.7378   0.4606    
ma1       -0.601809   0.013609  -44.2227   <2e-16 ***
intercept  9.999930   0.004079 2451.5779   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The extra AR coefficient is clearly not significant (z \approx 0.74). Moreover, the exact AICc of ARIMA(1,0,1) (28613.88) is 1.45 units higher than that of MA(1), and its BIC is 8.66 units higher, so auto.arima() returned a model that is worse by its own criterion. The reason is that, for long series, auto.arima() ranks the candidate models with an approximate likelihood during the search and re-estimates only the chosen model by exact maximum likelihood; under exact likelihood, MA(1) is better. The manual MA(1) is preferred.

Diagnostic Checking

checkresiduals(m06)

    Ljung-Box test

data:  Residuals from ARIMA(0,0,1) with non-zero mean
Q* = 6.9066, df = 9, p-value = 0.6468

Model df: 1.   Total lags used: 10

Ljung–Box: Q^* = 6.9066 with 9 degrees of freedom, p-value = 0.6468: the residuals are white noise.

res_06 <- m06$residuals
jarque.bera.test(res_06)

    Jarque Bera Test

data:  res_06
X-squared = 6.1596, df = 2, p-value = 0.04597

Jarque–Bera: X^2 = 6.1596, p-value = 0.04597, so normality is rejected at the 5% level, but only marginally. The series was simulated with Gaussian innovations, so this is a false rejection (a type I error): with n = 10{,}000 the test is very powerful, and at \alpha = 5\% about one in twenty correctly specified models is flagged by chance. The histogram is bell-shaped, and normality is not required for consistent estimation of ARIMA parameters (it matters mainly for exact prediction intervals).

Conclusion. ARIMA(0,0,1) with a non-zero mean is adequate for Y06.

Model 07

ACF & PACF Series 07

p_acf  <- ggAcf(df$Y07, lag.max = 36) + ggtitle("ACF Y07") 
p_pacf <- ggPacf(df$Y07, lag.max = 36) + ggtitle("PACF Y07") 

p_acf | p_pacf

Both functions tail off. The ACF oscillates with alternating signs and decays over about 20 lags, and the PACF shows large spikes at the first lags that die down irregularly over about 15 lags. Neither function cuts off, so Y07 follows a mixed ARMA(p,q) process whose orders cannot be read directly from the plots.

Stationarity Test

adf.test(df$Y07)

    Augmented Dickey-Fuller Test

data:  df$Y07
Dickey-Fuller = -19.546, Lag order = 21, p-value = 0.01
alternative hypothesis: stationary

ADF: Dickey–Fuller = -19.546, p-value \le 0.01, so the unit root is rejected.

pp.test(df$Y07)

    Phillips-Perron Unit Root Test

data:  df$Y07
Dickey-Fuller Z(alpha) = -18395, Truncation lag parameter = 12, p-value = 0.01
alternative hypothesis: stationary

PP: Z(\alpha) = -18395, p-value \le 0.01, so the unit root is rejected.

kpss.test(df$Y07)

    KPSS Test for Level Stationarity

data:  df$Y07
KPSS Level = 0.084136, Truncation lag parameter = 12, p-value = 0.1

KPSS: statistic = 0.084136, p-value \ge 0.10, so level stationarity is not rejected. Y07 is stationary (d = 0).

Order Selection for the Mixed Model

Because the ACF and the PACF both tail off, the orders are chosen by fitting every ARMA(p,q) model with p \le 4 and q \le 6 by exact maximum likelihood and ranking them by BIC. This exhaustive search over a small grid is the manual counterpart of auto.arima().

ic_grid <- function(y, d, p_max = 4, q_max = 6) {
  grid <- expand.grid(p = 0:p_max, q = 0:q_max)
  grid$AICc <- NA_real_
  grid$BIC  <- NA_real_
  for (i in seq_len(nrow(grid))) {
    fit <- tryCatch(Arima(y, order = c(grid$p[i], d, grid$q[i])),
                    error = function(e) NULL)
    if (!is.null(fit)) {
      grid$AICc[i] <- fit$aicc
      grid$BIC[i]  <- fit$bic
    }
  }
  grid %>%
    mutate(model = sprintf("ARIMA(%d,%d,%d)", p, d, q)) %>%
    arrange(BIC) %>%
    select(model, AICc, BIC)
}

ts_07   <- ts(data=df$Y07)
grid_07 <- ic_grid(ts_07, d = 0)
head(grid_07, 5)

The smallest BIC is attained by ARIMA(3,0,5), and the smallest AICc by ARIMA(3,0,5). The tentative model is therefore ARIMA(3,0,5) with a non-zero mean.

Parameter Estimation

m07 <- Arima(ts_07, order = c(3,0,5))
summary(m07)

Series: ts_07 
ARIMA(3,0,5) with non-zero mean 

Coefficients:
          ar1     ar2     ar3     ma1     ma2     ma3      ma4     ma5     mean
      -0.4989  0.0896  0.4993  -0.399  0.6059  0.6050  -0.3026  0.4107  59.9868
s.e.   0.0126  0.0155  0.0122   0.013  0.0120  0.0115   0.0123  0.0115   0.0419

sigma^2 = 3.944:  log likelihood = -21049.72
AIC=42119.45   AICc=42119.47   BIC=42191.55

Training set error measures:
                        ME     RMSE      MAE        MPE     MAPE      MASE         ACF1
Training set -0.0001913709 1.985106 1.588324 -0.1130719 2.675219 0.1980213 0.0005221438

The three AR coefficients, five MA coefficients and the mean (\hat\mu = 59.9868) are estimated precisely. The residual variance is \hat\sigma_a^2 = 3.944, with AIC = 42119.45, AICc = 42119.47 and BIC = 42191.55.

coeftest(m07)

z test of coefficients:

           Estimate Std. Error   z value  Pr(>|z|)    
ar1       -0.498897   0.012605  -39.5805 < 2.2e-16 ***
ar2        0.089606   0.015477    5.7896 7.055e-09 ***
ar3        0.499349   0.012229   40.8321 < 2.2e-16 ***
ma1       -0.399026   0.012969  -30.7669 < 2.2e-16 ***
ma2        0.605871   0.012008   50.4536 < 2.2e-16 ***
ma3        0.604994   0.011511   52.5577 < 2.2e-16 ***
ma4       -0.302583   0.012259  -24.6827 < 2.2e-16 ***
ma5        0.410708   0.011482   35.7685 < 2.2e-16 ***
intercept 59.986830   0.041872 1432.6073 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

All coefficients are significant; the smallest test statistic is z = 5.79 for \phi_2 (p-value \approx 7\times10^{-9}).

Fitted Model Equation

The general form of ARIMA(3,0,5) with mean \mu is

(1-\phi_1B-\phi_2B^2-\phi_3B^3)\,(Y_t-\mu) = (1+\theta_1B+\theta_2B^2+\cdots+\theta_5B^5)\,a_t .

With the estimates,

\begin{aligned} (1 + 0.4989B - 0.0896B^2 - 0.4993B^3)\,(Y_t - 59.9868) ={}& (1 - 0.3990B + 0.6059B^2 + 0.6050B^3 \\ &- 0.3026B^4 + 0.4107B^5)\,a_t , \end{aligned}

which expands to

\begin{aligned} Y_t ={}& 54.5845 - 0.4989\,Y_{t-1} + 0.0896\,Y_{t-2} + 0.4993\,Y_{t-3} \\ &+ a_t - 0.3990\,a_{t-1} + 0.6059\,a_{t-2} + 0.6050\,a_{t-3} - 0.3026\,a_{t-4} + 0.4107\,a_{t-5}, \qquad a_t \sim \mathcal{N}(0,\ 3.944), \end{aligned}

where 54.5845 = 59.9868 \times (1 - 0.0901) and 0.0901 is the sum of the three AR coefficients.

Comparison with auto.arima()

aa07 <- auto.arima(ts_07, max.p = 10, max.q = 10)
aa07

Series: ts_07 
ARIMA(10,0,3) with non-zero mean 

Coefficients:
          ar1      ar2     ar3     ar4      ar5      ar6      ar7     ar8     ar9    ar10     ma1     ma2     ma3     mean
      -1.2787  -0.1822  1.2038  0.9888  -0.0655  -0.6036  -0.3578  0.0641  0.1723  0.0848  0.3828  0.1757  0.5486  59.9855
s.e.   0.0379   0.0440  0.0313  0.0641   0.0513   0.0421   0.0490  0.0267  0.0250  0.0174  0.0367  0.0214  0.0199   0.0430

sigma^2 = 3.95:  log likelihood = -21055.18
AIC=42140.35   AICc=42140.4   BIC=42248.51

auto.arima() selects ARIMA(10,0,3) with a non-zero mean, with 13 ARMA coefficients instead of 8. It is worse by both criteria: AICc = 42140.40 (20.9 units larger than for ARIMA(3,0,5)) and BIC = 42248.51 (57.0 units larger). The stepwise search, run on an approximate likelihood, drifted towards a long AR polynomial that only approximates the true structure, whereas the exhaustive grid search with exact likelihood found the better and more parsimonious ARIMA(3,0,5).

Diagnostic Checking

checkresiduals(m07)

    Ljung-Box test

data:  Residuals from ARIMA(3,0,5) with non-zero mean
Q* = 2.8719, df = 3, p-value = 0.4118

Model df: 8.   Total lags used: 11

Ljung–Box: Q^* = 2.8719 with 3 degrees of freedom, p-value = 0.4118: no remaining autocorrelation.

res_07 <- m07$residuals
jarque.bera.test(res_07)

    Jarque Bera Test

data:  res_07
X-squared = 0.62123, df = 2, p-value = 0.733

Jarque–Bera: X^2 = 0.62123, p-value = 0.733: normality is not rejected.

Conclusion. ARIMA(3,0,5) with a non-zero mean is adequate for Y07.

Model 08

ACF & PACF Series 08

As for Y01, the ACF of the levels is shown up to lag 5,000.

p_acf  <- ggAcf(df$Y08, lag.max = 5000) + ggtitle("ACF Y08") 
p_pacf <- ggPacf(df$Y08, lag.max = 36) + ggtitle("PACF Y08") 

p_acf | p_pacf

The ACF decays very slowly, with long waves that persist up to lag 5,000, and the PACF has a spike close to 1 at lag 1. Both are signs of a non-stationary series.

Stationarity Test

adf.test(df$Y08)

    Augmented Dickey-Fuller Test

data:  df$Y08
Dickey-Fuller = -2.2797, Lag order = 21, p-value = 0.4599
alternative hypothesis: stationary

ADF: Dickey–Fuller = -2.2797, p-value = 0.4599 > 0.05, so the unit root is not rejected.

pp.test(df$Y08)

    Phillips-Perron Unit Root Test

data:  df$Y08
Dickey-Fuller Z(alpha) = -23.055, Truncation lag parameter = 12, p-value = 0.03844
alternative hypothesis: stationary

PP: Z(\alpha) = -23.055, p-value = 0.03844 < 0.05, so this test rejects the unit root at the 5% level, in contrast with the ADF test.

kpss.test(df$Y08)

    KPSS Test for Level Stationarity

data:  df$Y08
KPSS Level = 10.338, Truncation lag parameter = 12, p-value = 0.01

KPSS: statistic = 10.338, p-value \le 0.01, so level stationarity is rejected. Two of the three tests, together with the plots, indicate non-stationarity. The PP test corrects for serial correlation non-parametrically (here with only 12 lags), and it is known to over-reject a true unit root when the differenced series has strong short-run dynamics, especially moving-average components; the ADF test handles such dynamics through its 21 augmentation lags. We therefore conclude that Y08 is non-stationary and difference it.

Differencing

ts_08 <- ts(data=df$Y08)
diff_ts_08 <- diff(ts_08, lag = 1)
autoplot(ts_08) | autoplot(diff_ts_08)

The original series wanders between about 80 and 250 without a fixed mean, while the differenced series W_t = Y_t - Y_{t-1} fluctuates around zero with constant variance.

adf.test(diff_ts_08)

    Augmented Dickey-Fuller Test

data:  diff_ts_08
Dickey-Fuller = -23.467, Lag order = 21, p-value = 0.01
alternative hypothesis: stationary

Dickey–Fuller = -23.467, p-value \le 0.01: W_t is stationary and d = 1.

p08_acf_diff  <- ggAcf(diff_ts_08, lag.max = 36) + ggtitle("ACF diff(Y08)") 
p08_pacf_diff <- ggPacf(diff_ts_08, lag.max = 36) + ggtitle("PACF diff(Y08)") 

p08_acf_diff | p08_pacf_diff

The ACF of W_t is a damped sine wave (alternating groups of positive and negative autocorrelations that fade out within about 20 lags), and the PACF also dies down: large spikes at lags 1 to 4 are followed by smaller significant spikes up to about lag 10. Both functions tail off, so W_t follows a mixed ARMA(p,q) process.

Order Selection for the Mixed Model

As for Y07, all models with p \le 4 and q \le 6 are compared, now with d = 1 (the function ic_grid() was defined in Model 07).

grid_08 <- ic_grid(ts_08, d = 1)
head(grid_08, 5)

The smallest BIC is attained by ARIMA(2,1,3), and the smallest AICc by ARIMA(4,1,5). The tentative model is ARIMA(2,1,3) without a constant.

Parameter Estimation

m08 <- Arima(ts_08, order = c(2,1,3))
summary(m08)

Series: ts_08 
ARIMA(2,1,3) 

Coefficients:
         ar1      ar2     ma1      ma2      ma3
      0.6846  -0.6997  0.6033  -0.1857  -0.4097
s.e.  0.0086   0.0083  0.0108   0.0137   0.0105

sigma^2 = 2.251:  log likelihood = -18245.18
AIC=36502.36   AICc=36502.37   BIC=36545.62

Training set error measures:
                       ME     RMSE      MAE         MPE      MAPE      MASE        ACF1
Training set -0.005375228 1.500045 1.201171 -0.00893472 0.8403576 0.3886642 0.003783089

The estimates are \hat\phi_1 = 0.6846, \hat\phi_2 = -0.6997, \hat\theta_1 = 0.6033, \hat\theta_2 = -0.1857 and \hat\theta_3 = -0.4097, with \hat\sigma_a^2 = 2.251, AIC = 36502.36, AICc = 36502.37 and BIC = 36545.62.

coeftest(m08)

z test of coefficients:

      Estimate Std. Error z value  Pr(>|z|)    
ar1  0.6846163  0.0085716  79.870 < 2.2e-16 ***
ar2 -0.6997123  0.0082530 -84.783 < 2.2e-16 ***
ma1  0.6033251  0.0107631  56.055 < 2.2e-16 ***
ma2 -0.1857441  0.0136597 -13.598 < 2.2e-16 ***
ma3 -0.4097409  0.0104906 -39.058 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

All coefficients are significant; the smallest test statistic is |z| = 13.60 for \theta_2.

Fitted Model Equation

The general form of ARIMA(2,1,3) is

(1-\phi_1B-\phi_2B^2)\,(1-B)\,Y_t = (1+\theta_1B+\theta_2B^2+\theta_3B^3)\,a_t .

With the estimates,

(1 - 0.6846B + 0.6997B^2)\,(1-B)\,Y_t = (1 + 0.6033B - 0.1857B^2 - 0.4097B^3)\,a_t .

Multiplying out the AR side, (1-\phi_1B-\phi_2B^2)(1-B) = 1-(1+\phi_1)B+(\phi_1-\phi_2)B^2+\phi_2B^3 = 1 - 1.6846B + 1.3843B^2 - 0.6997B^3, so the explicit form is

\begin{aligned} Y_t ={}& 1.6846\,Y_{t-1} - 1.3843\,Y_{t-2} + 0.6997\,Y_{t-3} \\ &+ a_t + 0.6033\,a_{t-1} - 0.1857\,a_{t-2} - 0.4097\,a_{t-3}, \qquad a_t \sim \mathcal{N}(0,\ 2.251). \end{aligned}

The coefficients of Y_{t-1}, Y_{t-2} and Y_{t-3} sum to one (1.6846 - 1.3843 + 0.6997 = 1), which reflects the unit root.

Comparison with auto.arima()

aa08 <- auto.arima(ts_08, max.p = 10, max.q = 10)
aa08

Series: ts_08 
ARIMA(7,1,3) 

Coefficients:
         ar1      ar2     ar3      ar4      ar5     ar6      ar7     ma1      ma2      ma3
      0.7348  -0.8512  0.1217  -0.0828  -0.0460  0.0474  -0.0403  0.5577  -0.1055  -0.3369
s.e.  0.0937   0.0398  0.0793   0.0588   0.0612  0.0289   0.0244  0.0936   0.1229   0.0745

sigma^2 = 2.249:  log likelihood = -18236.57
AIC=36495.13   AICc=36495.16   BIC=36574.44

auto.arima() selects ARIMA(7,1,3), with 10 coefficients instead of 5. Its AICc (36495.16) is 7.2 units lower than that of ARIMA(2,1,3), but its BIC (36574.44) is 28.8 units higher.

coeftest(aa08)

z test of coefficients:

     Estimate Std. Error  z value  Pr(>|z|)    
ar1  0.734817   0.093741   7.8388 4.549e-15 ***
ar2 -0.851172   0.039796 -21.3886 < 2.2e-16 ***
ar3  0.121681   0.079311   1.5342   0.12498    
ar4 -0.082842   0.058824  -1.4083   0.15904    
ar5 -0.046002   0.061197  -0.7517   0.45223    
ar6  0.047368   0.028917   1.6381   0.10141    
ar7 -0.040307   0.024363  -1.6544   0.09804 .  
ma1  0.557676   0.093550   5.9612 2.503e-09 ***
ma2 -0.105457   0.122889  -0.8581   0.39081    
ma3 -0.336931   0.074454  -4.5254 6.029e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Six of the ten coefficients (ar3 to ar7 and ma2) are not significant at the 5% level. AICc penalises each extra parameter by about 2, whereas BIC penalises it by \ln n \approx 9.21. The additional AR terms absorb the tiny residual autocorrelation discussed below and improve the likelihood slightly, but at the cost of six insignificant coefficients. Parsimony and BIC favour ARIMA(2,1,3), which is also the model the series was generated from; the automatic model over-fits.

Diagnostic Checking

checkresiduals(m08)

    Ljung-Box test

data:  Residuals from ARIMA(2,1,3)
Q* = 15.186, df = 5, p-value = 0.009597

Model df: 5.   Total lags used: 10

Ljung–Box: Q^* = 15.186 with 5 degrees of freedom, p-value = 0.0096, so white noise is rejected at the 5% level. The residual ACF shows why: only a handful of very small autocorrelations, the largest at lag 7 (about 0.024), lie just outside the \pm 0.0196 bands. With n = 10{,}000 the Ljung–Box test detects even negligible dependence (an autocorrelation of 0.024 explains less than 0.1% of the residual variance). Because the series was simulated from ARIMA(2,1,3), this rejection is a false alarm (a type I error). In a real application one would try adding terms, but the over-fitted auto.arima() model above shows that such terms are insignificant and worsen BIC.

res_08 <- m08$residuals
jarque.bera.test(res_08)

    Jarque Bera Test

data:  res_08
X-squared = 0.34746, df = 2, p-value = 0.8405

Jarque–Bera: X^2 = 0.34746, p-value = 0.8405: normality is not rejected.

Conclusion. ARIMA(2,1,3) is retained for Y08. All coefficients are significant and the residuals are normal; the Ljung–Box rejection reflects autocorrelations too small to matter in practice.

Model 09

ACF & PACF Series 09

p_acf  <- ggAcf(df$Y09, lag.max = 5000) + ggtitle("ACF Y09") 
p_pacf <- ggPacf(df$Y09, lag.max = 36) + ggtitle("PACF Y09") 

p_acf | p_pacf

The ACF of the levels decays slowly with long waves, and the PACF has a dominant spike close to 1 at lag 1 followed by smaller spikes at lags 2 to 5. The slow ACF decay again indicates non-stationarity.

Stationarity Test

adf.test(df$Y09)

    Augmented Dickey-Fuller Test

data:  df$Y09
Dickey-Fuller = -2.4453, Lag order = 21, p-value = 0.3897
alternative hypothesis: stationary

ADF: Dickey–Fuller = -2.4453, p-value = 0.3897 > 0.05, so the unit root is not rejected.

pp.test(df$Y09)

    Phillips-Perron Unit Root Test

data:  df$Y09
Dickey-Fuller Z(alpha) = -23.056, Truncation lag parameter = 12, p-value = 0.03843
alternative hypothesis: stationary

PP: Z(\alpha) = -23.056, p-value = 0.03843 < 0.05, so the PP test again rejects the unit root, unlike the ADF test.

kpss.test(df$Y09)

    KPSS Test for Level Stationarity

data:  df$Y09
KPSS Level = 35.883, Truncation lag parameter = 12, p-value = 0.01

KPSS: statistic = 35.883, p-value \le 0.01, so level stationarity is rejected. As for Y08, the ADF and KPSS tests and the plots agree, while the PP test is affected by the strong short-run dynamics of the differenced series. Y09 is non-stationary and is differenced.

Differencing

ts_09 <- ts(data=df$Y09)
diff_ts_09 <- diff(ts_09, lag = 1)
autoplot(ts_09) | autoplot(diff_ts_09)

The original series wanders without a fixed mean, while the differenced series fluctuates around zero with constant variance.

adf.test(diff_ts_09)

    Augmented Dickey-Fuller Test

data:  diff_ts_09
Dickey-Fuller = -23.161, Lag order = 21, p-value = 0.01
alternative hypothesis: stationary

Dickey–Fuller = -23.161, p-value \le 0.01: W_t is stationary and d = 1.

p09_acf_diff  <- ggAcf(diff_ts_09, lag.max = 36) + ggtitle("ACF diff(Y09)") 
p09_pacf_diff <- ggPacf(diff_ts_09, lag.max = 36) + ggtitle("PACF diff(Y09)") 

p09_acf_diff | p09_pacf_diff

The ACF of W_t is a slowly damped sine wave (with a period of about 8 lags) that is still visible beyond lag 30, whereas the PACF has four significant spikes (lags 1 to 4) and cuts off afterwards. A tailing-off ACF with a PACF that cuts off after lag 4 identifies an AR(4) process for W_t, so the tentative model is ARIMA(4,1,0) without a constant.

Parameter Estimation

m09 <- Arima(ts_09, order = c(4,1,0))
summary(m09)

Series: ts_09 
ARIMA(4,1,0) 

Coefficients:
         ar1      ar2      ar3      ar4
      0.3991  -0.2036  -0.2096  -0.3379
s.e.  0.0094   0.0100   0.0100   0.0094

sigma^2 = 1.035:  log likelihood = -14358.39
AIC=28726.79   AICc=28726.79   BIC=28762.84

Training set error measures:
                      ME     RMSE       MAE          MPE      MAPE      MASE         ACF1
Training set 0.004983071 1.017068 0.8110864 0.0007137312 0.1600883 0.6847773 0.0006931769

The estimates are \hat\phi_1 = 0.3991, \hat\phi_2 = -0.2036, \hat\phi_3 = -0.2096 and \hat\phi_4 = -0.3379, with \hat\sigma_a^2 = 1.035, AIC = AICc = 28726.79 and BIC = 28762.84.

coeftest(m09)

z test of coefficients:

      Estimate Std. Error z value  Pr(>|z|)    
ar1  0.3990913  0.0094106  42.409 < 2.2e-16 ***
ar2 -0.2036409  0.0100131 -20.337 < 2.2e-16 ***
ar3 -0.2096380  0.0100126 -20.937 < 2.2e-16 ***
ar4 -0.3378805  0.0094210 -35.864 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

All four AR coefficients are highly significant (|z| \ge 20.3).

Fitted Model Equation

The general form of ARIMA(4,1,0) is

(1-\phi_1B-\phi_2B^2-\phi_3B^3-\phi_4B^4)\,(1-B)\,Y_t = a_t .

With the estimates,

(1 - 0.3991B + 0.2036B^2 + 0.2096B^3 + 0.3379B^4)\,(1-B)\,Y_t = a_t .

Multiplying out the AR side gives

\begin{aligned} &1-(1+\phi_1)B+(\phi_1-\phi_2)B^2+(\phi_2-\phi_3)B^3+(\phi_3-\phi_4)B^4+\phi_4B^5 \\ &= 1 - 1.3991B + 0.6027B^2 + 0.0060B^3 + 0.1282B^4 - 0.3379B^5 , \end{aligned}

so the explicit form is

\begin{aligned} Y_t ={}& 1.3991\,Y_{t-1} - 0.6027\,Y_{t-2} - 0.0060\,Y_{t-3} - 0.1282\,Y_{t-4} + 0.3379\,Y_{t-5} + a_t, \qquad a_t \sim \mathcal{N}(0,\ 1.035). \end{aligned}

As for Y08, the five autoregressive coefficients sum to one (up to rounding), which reflects the unit root.

Comparison with auto.arima()

aa09 <- auto.arima(ts_09, max.p = 10, max.q = 10)
aa09

Series: ts_09 
ARIMA(4,1,0) 

Coefficients:
         ar1      ar2      ar3      ar4
      0.3991  -0.2036  -0.2096  -0.3379
s.e.  0.0094   0.0100   0.0100   0.0094

sigma^2 = 1.035:  log likelihood = -14358.39
AIC=28726.79   AICc=28726.79   BIC=28762.84

coeftest(aa09)

z test of coefficients:

      Estimate Std. Error z value  Pr(>|z|)    
ar1  0.3990913  0.0094106  42.409 < 2.2e-16 ***
ar2 -0.2036409  0.0100131 -20.337 < 2.2e-16 ***
ar3 -0.2096380  0.0100126 -20.937 < 2.2e-16 ***
ar4 -0.3378805  0.0094210 -35.864 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

auto.arima() selects exactly the same model, ARIMA(4,1,0), with identical estimates and significance tests. Both approaches agree.

Diagnostic Checking

checkresiduals(m09)

    Ljung-Box test

data:  Residuals from ARIMA(4,1,0)
Q* = 3.7051, df = 6, p-value = 0.7165

Model df: 4.   Total lags used: 10

Ljung–Box: Q^* = 3.7051 with 6 degrees of freedom, p-value = 0.7165: the residuals are white noise.

res_09 <- m09$residuals
jarque.bera.test(res_09)

    Jarque Bera Test

data:  res_09
X-squared = 0.11057, df = 2, p-value = 0.9462

Jarque–Bera: X^2 = 0.11057, p-value = 0.9462: normality is not rejected.

Conclusion. ARIMA(4,1,0) is adequate for Y09.

Model 10

ACF & PACF Series 10

p_acf  <- ggAcf(df$Y10, lag.max = 36) + ggtitle("ACF Y10") 
p_pacf <- ggPacf(df$Y10, lag.max = 36) + ggtitle("PACF Y10") 

p_acf | p_pacf

The ACF is a damped sine wave: positive at lags 1 and 2, negative at lags 3 to 6, positive again around lags 8 to 10, and negligible after about lag 12. The PACF has two significant spikes (about 0.76 at lag 1 and -0.60 at lag 2) and cuts off afterwards. This is the pattern of an AR(2) process with complex roots.

Stationarity Test

adf.test(df$Y10)

    Augmented Dickey-Fuller Test

data:  df$Y10
Dickey-Fuller = -21.278, Lag order = 21, p-value = 0.01
alternative hypothesis: stationary

ADF: Dickey–Fuller = -21.278, p-value \le 0.01, so the unit root is rejected.

pp.test(df$Y10)

    Phillips-Perron Unit Root Test

data:  df$Y10
Dickey-Fuller Z(alpha) = -1452.6, Truncation lag parameter = 12, p-value = 0.01
alternative hypothesis: stationary

PP: Z(\alpha) = -1452.6, p-value \le 0.01, so the unit root is rejected.

kpss.test(df$Y10)

    KPSS Test for Level Stationarity

data:  df$Y10
KPSS Level = 0.09477, Truncation lag parameter = 12, p-value = 0.1

KPSS: statistic = 0.09477, p-value \ge 0.10, so level stationarity is not rejected. Y10 is stationary (d = 0), and the tentative model is ARIMA(2,0,0) with a non-zero mean.

Parameter Estimation

ts_10 <- ts(data=df$Y10)
m10 <- Arima(ts_10, order = c(2,0,0))
summary(m10)

Series: ts_10 
ARIMA(2,0,0) with non-zero mean 

Coefficients:
         ar1      ar2     mean
      1.2131  -0.6018  49.9794
s.e.  0.0080   0.0080   0.0520

sigma^2 = 4.094:  log likelihood = -21235.86
AIC=42479.73   AICc=42479.73   BIC=42508.57

Training set error measures:
                       ME     RMSE      MAE        MPE     MAPE      MASE        ACF1
Training set 9.326962e-05 2.022957 1.621585 -0.1660213 3.268997 0.7485863 0.003651908

The estimates are \hat\phi_1 = 1.2131, \hat\phi_2 = -0.6018 and \hat\mu = 49.9794, with \hat\sigma_a^2 = 4.094, AIC = AICc = 42479.73 and BIC = 42508.57.

coeftest(m10)

z test of coefficients:

            Estimate Std. Error z value  Pr(>|z|)    
ar1        1.2130912  0.0079855 151.911 < 2.2e-16 ***
ar2       -0.6018124  0.0079897 -75.323 < 2.2e-16 ***
intercept 49.9794164  0.0520412 960.381 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

All parameters are highly significant (z = 151.91, -75.32 and 960.38).

Fitted Model Equation

The general form of ARIMA(2,0,0) with mean \mu is

(1-\phi_1B-\phi_2B^2)\,(Y_t-\mu) = a_t .

With the estimates,

(1 - 1.2131B + 0.6018B^2)\,(Y_t - 49.9794) = a_t ,

which expands to

Y_t = 19.4281 + 1.2131\,Y_{t-1} - 0.6018\,Y_{t-2} + a_t,\qquad a_t \sim \mathcal{N}(0,\ 4.094),

where 19.4281 = 49.9794 \times (1 - 1.2131 + 0.6018).

Comparison with auto.arima()

aa10 <- auto.arima(ts_10, max.p = 10, max.q = 10)
aa10

Series: ts_10 
ARIMA(2,0,0) with non-zero mean 

Coefficients:
         ar1      ar2     mean
      1.2131  -0.6018  49.9794
s.e.  0.0080   0.0080   0.0520

sigma^2 = 4.094:  log likelihood = -21235.86
AIC=42479.73   AICc=42479.73   BIC=42508.57

auto.arima() selects the same model, ARIMA(2,0,0) with a non-zero mean, with identical estimates. Both approaches agree.

Diagnostic Checking

checkresiduals(m10)

    Ljung-Box test

data:  Residuals from ARIMA(2,0,0) with non-zero mean
Q* = 7.8159, df = 8, p-value = 0.4517

Model df: 2.   Total lags used: 10

Ljung–Box: Q^* = 7.8159 with 8 degrees of freedom, p-value = 0.4517: the residuals are white noise.

res_10 <- m10$residuals
jarque.bera.test(res_10)

    Jarque Bera Test

data:  res_10
X-squared = 1.3388, df = 2, p-value = 0.512

Jarque–Bera: X^2 = 1.3388, p-value = 0.512: normality is not rejected.

Conclusion. ARIMA(2,0,0) with a non-zero mean is adequate for Y10.

Manual vs auto.arima()

The table below collects, for every series, the manually identified model and the model chosen by auto.arima(), together with their AICc and BIC, the number of insignificant coefficients (p-value > 0.05) in the automatic model, and the Ljung–Box p-values computed with the same lag choice as checkresiduals().

manual <- list(Y01 = m01, Y02 = m02, Y03 = m03, Y04 = m04, Y05 = m05,
               Y06 = m06, Y07 = m07, Y08 = m08, Y09 = m09, Y10 = m10)
auto   <- list(Y01 = aa01, Y02 = aa02, Y03 = aa03, Y04 = aa04, Y05 = aa05,
               Y06 = aa06, Y07 = aa07, Y08 = aa08, Y09 = aa09, Y10 = aa10)

model_label <- function(fit) {
  o <- arimaorder(fit)
  sprintf("ARIMA(%d,%d,%d)", o[1], o[2], o[3])
}
n_insignificant <- function(fit) {
  ct <- coeftest(fit)
  sum(ct[!rownames(ct) %in% c("intercept", "drift"), 4] > 0.05)
}
ljung_box_p <- function(fit) {
  k <- sum(fit$arma[1:2])                       # number of AR + MA coefficients
  Box.test(residuals(fit), lag = max(k + 3, 10), type = "Ljung-Box", fitdf = k)$p.value
}

comparison <- tibble(
  Series                       = names(manual),
  Manual                       = map_chr(manual, model_label),
  `auto.arima`                 = map_chr(auto, model_label),
  Agree                        = Manual == `auto.arima`,
  `AICc manual`                = map_dbl(manual, "aicc"),
  `AICc auto`                  = map_dbl(auto, "aicc"),
  `BIC manual`                 = map_dbl(manual, "bic"),
  `BIC auto`                   = map_dbl(auto, "bic"),
  `Insignificant coef. (auto)` = map_int(auto, n_insignificant),
  `Ljung-Box p (manual)`       = map_dbl(manual, ljung_box_p),
  `Ljung-Box p (auto)`         = map_dbl(auto, ljung_box_p)
) %>%
  mutate(across(where(is.double), ~ round(.x, 3)))

comparison

Agreement. auto.arima() reproduces the manual model exactly for four series: Y02 (AR(1)), Y04 (MA(2)), Y09 (ARIMA(4,1,0)) and Y10 (AR(2)). These are low-order models with simple ACF/PACF patterns, and both approaches give identical estimates. The order of differencing chosen by auto.arima() is correct for all ten series, so its KPSS-based choice of d agrees with the unit-root tests used here.

One extra, insignificant term. For Y05 and Y06, auto.arima() adds one parameter that is not significant (ar7 in AR(7), z \approx 1.85; ar1 in ARMA(1,1), z \approx 0.74). For Y05, the AICc of AR(7) is only 1.41 units lower, a difference below 2 that does not indicate a real improvement, and BIC prefers AR(6). For Y06, the selected ARMA(1,1) even has a higher exact AICc (by 1.45) and BIC (by 8.66) than MA(1).

Several extra terms. For Y08, ARIMA(7,1,3) has a lower AICc than ARIMA(2,1,3) (by 7.21) but a higher BIC (by 28.82) and six insignificant coefficients. The extra AR terms model a tiny amount of residual autocorrelation that is statistically detectable only because n is so large.

Over-parameterised local optima. For Y01, Y03 and Y07, auto.arima() returns large models, ARIMA(9,1,6), ARIMA(10,0,2) and ARIMA(10,0,3), whose AICc is worse than that of the manual model by 120.2, 409.4 and 20.9 units, respectively (BIC is worse by 170.6, 438.2 and 57.0 units). Here the automatic search failed by its own standard, because it did not find the model with the smallest AICc.

Conclusion. The manual Box–Jenkins approach recovered the data-generating model for all ten series and produced parsimonious models with significant coefficients and well-behaved residuals. It combined cut-off patterns in the ACF and PACF, significance tests, residual diagnostics, and an exhaustive information-criterion search reserved for the genuinely mixed models. auto.arima() is a valuable and fast benchmark, especially for low-order models, but its output should be checked rather than accepted blindly. A more thorough automatic search is available through auto.arima(y, stepwise = FALSE, approximation = FALSE, max.order = 15, ic = "bic"), at a much higher computational cost for series of this length.

Summary