Preprocessing

The tables below show the data and the summary of its variables with the missing values replaced by the mean of the variable they belong to.

Data Frame Summary

finalPro_data

Dimensions: 58 x 9
Duplicates: 0

No	Variable	Stats / Values	Freqs (% of Valid)	Graph	Valid	Missing
1	Country [factor]	1. Afghanistan 2. Albania 3. Algeria 4. American Samoa 5. Andorra 6. Angola 7. Antigua and Barbuda 8. Arab World 9. Argentina 10. Armenia [ 253 others ]	0 ( 0.0% ) 0 ( 0.0% ) 0 ( 0.0% ) 0 ( 0.0% ) 0 ( 0.0% ) 0 ( 0.0% ) 0 ( 0.0% ) 0 ( 0.0% ) 0 ( 0.0% ) 0 ( 0.0% ) 58 ( 100.0% )		58 (100.0%)	0 (0.0%)
2	Code [factor]	1. ABW 2. AFG 3. AGO 4. ALB 5. AND 6. ARB 7. ARE 8. ARG 9. ARM 10. ASM [ 253 others ]	0 ( 0.0% ) 0 ( 0.0% ) 0 ( 0.0% ) 0 ( 0.0% ) 0 ( 0.0% ) 0 ( 0.0% ) 0 ( 0.0% ) 0 ( 0.0% ) 0 ( 0.0% ) 0 ( 0.0% ) 58 ( 100.0% )		58 (100.0%)	0 (0.0%)
3	Year [numeric]	Mean (sd) : 1988.5 (16.9) min ≤ med ≤ max: 1960 ≤ 1988.5 ≤ 2017 IQR (CV) : 28.5 (0)	58 distinct values		58 (100.0%)	0 (0.0%)
4	GDP [numeric]	Mean (sd) : 939447947 (615865973) min ≤ med ≤ max: 112155599 ≤ 934787385 ≤ 2195599557 IQR (CV) : 1061657866 (0.7)	58 distinct values		58 (100.0%)	0 (0.0%)
5	Growth [numeric]	Mean (sd) : 1.2 (6.3) min ≤ med ≤ max: -36.7 ≤ 2 ≤ 9.5 IQR (CV) : 4.5 (5.2)	58 distinct values		58 (100.0%)	0 (0.0%)
6	CPI [numeric]	Mean (sd) : 73.6 (23.9) min ≤ med ≤ max: 36.8 ≤ 73.6 ≤ 186.9 IQR (CV) : 14.3 (0.3)	37 distinct values		58 (100.0%)	0 (0.0%)
7	Imports [numeric]	Mean (sd) : 30.8 (7.1) min ≤ med ≤ max: 18 ≤ 29.7 ≤ 46.9 IQR (CV) : 11.6 (0.2)	58 distinct values		58 (100.0%)	0 (0.0%)
8	Exports [numeric]	Mean (sd) : 20.7 (5.9) min ≤ med ≤ max: 10.7 ≤ 21.7 ≤ 34.3 IQR (CV) : 8.9 (0.3)	58 distinct values		58 (100.0%)	0 (0.0%)
9	Population [numeric]	Mean (sd) : 2974790 (1052736) min ≤ med ≤ max: 1503508 ≤ 2842456 ≤ 4659080 IQR (CV) : 1977168 (0.4)	58 distinct values		58 (100.0%)	0 (0.0%)

Generated by summarytools 1.0.1 (R version 4.2.0)
2022-06-24

Data Visualization

The time plot in Figure 1 shows some non-stationarity, with an overall decline. There appears to be an improvement in 1994 but it is followed by further economic decline after 1997. The plot does not show any strong evidence of changing variance, so doing a Log or Box-Cox transformation is not necessary. In addition, it appears that the time series is non-stationary due to trends and changing levels. However, this is further proven below.

ACF and PACF

The autocorrelation function (ACF), \(\hat \rho (h)\), which is the correlation between time series with a lagged version of itself, and partial autocorrelation function (PACF) or the conditional correlation of the time series with a lag of itself with the linear dependence of all the lags between them removed, both assume stationarity of the underlying time series. By using the ACF and PACF, an appropriate forecasting method or model can be selected.

Using the ACF shown in Figure 2 below, it is evident that the time series data Exports is non-stationary since the ACF is slowly decreasing and mostly remains above the significance range, which is shown as the blue dashed line in the plot. The first PACF of Exports \(\hat\rho_{11} \approx 1\), which is indicative of non-stationarity. Thus, it can be concluded that Exports is non-stationary using the ACF and PACF plots below and methods such as differencing and transformations can be implemented to transform it into a stationary form.

Figure 2: Time plot, ACF, and PACF plots for the Central African Republic Exports

Stationarity

In order to further determine if the time series Exports is stationary, the Augmented Dickey-Fuller (ADF) test can be performed. The following can be used to identify if the time series is indeed stationary.

• p-value \(> 0.05\): Fail to reject the null hypothesis (\(H_0\)), the data has a unit root and is non-stationary

• p-value \(\le 0.05\): Reject the null hypothesis (\(H_o\)) and accept the alternative hypothesis (\(H_a\)), the data does not have a unit root and is stationary

Since the p-value \(= 0.1006 \gt 0.05\) using ADF test, it can be concluded that Exports is non-stationary.

Another numerical method that can be used to check for the stationarity of Exports is the Kwiatkowski-Phillips-Schmidt-Shin (KPSS) Test for Stationarity.

• p-value \(< 0.05\): Fail to reject the null hypothesis (\(H_0\)), the data has a unit root and is non-stationary

• p-value \(\ge 0.05\): Reject the null hypothesis (\(H_o\)) and accept the alternative hypthesis (\(H_a\)), the data does not have a unit root and is stationary

Since the p-value \(= 0.01 \lt 0.05\) using KPSS test, it can be concluded that Exports is non-stationary.

## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_df[, "Exports"]
## Dickey-Fuller = -3.1744, Lag order = 3, p-value = 0.1006
## alternative hypothesis: stationary

## 
##  KPSS Test for Level Stationarity
## 
## data:  ts_df[, "Exports"]
## KPSS Level = 1.2824, Truncation lag parameter = 3, p-value = 0.01

Transformations

In the previous part, it is determined that transformations are not necessary for the time series Exports since the time plot in Figure 1 does not show any evidence of changing variance. In order to eliminate any non-constant variance or heteroskedasticity, performing transformation is necessary. Therefore, it must be further investigated as shown below.

Log and Box-cox Transformation

The time plot of transformed and original Exports is shown below in Figure 3. This plot shows that transformation does not have a significant effect on the time series Exports since the three plots appear similar. Therefore, it can be concluded that transformation is not necessary. Additionally, the ACF and PACF of Exports transformed using log and Box-cox are similar to the ACF and PACF of non-transformed Exports, which further support this conclusion.

## optimal lambda: -0.3650884

Figure 3: Transformed Exports of the Central African Republic

Log ACF and PACF

Figure 4: Time plot, ACF, and PACF plots for the transformed Central African Republic Exports

Box-cox ACF and PACF

Figure 5: Time plot, ACF, and PACF plots for the transformed Central African Republic Exports

Decompose

The time series Exports has only one seasonal cycle since it is a yearly data. Therefore, decomposing it is not necessary. However, if the seasonal cycle of Exports is 2 the trend, seasonality and error would appear as shown in Figure 6 below. For this project, the data is analyzed as non-seasonal.

Figure 6: Decomposition of Additive Exports

Differencing

Firs-order difference

In order to address the non-stationarity in Exports identified using the methods in the previous part, the first-order difference of the data can be used such that \(y'_t = y_t - y_{t-1}\), which is is the change between consecutive observations in the original series Exports. The differenced time series Exports is shown in Figure 7 below, which only includes \(T-1\) values since it is not possible to compute the difference \(y'_t\) for the first observation.

Figure 7 shows that Exports now appear to be stationary. The ACF below shows exponential decay, which is indicative of a stationary time series. The PACF shown in Figure 7 suggests an AR(2) model for the Exports since the PACF is tailing off at lag 2 where there are two spikes outside the threshold limit (blue dashed line). Hence, an ARIMA(2,1,0) is an initial candidate model. Meanwhile, the ACF suggests an MA(3) model since the ACF is cutting off sharply at lag 3 and negative at lag 1. Therefore, an alternative candidate is an ARIMA(0,1,3).

Figure 7: Time plot, ACF, and PACF plots for the differenced Central African Republic Exports

By conducting the KPSS test on the differenced Exports time series data, it can be further proven whether the differenced time series is indeed stationary. Based on the results shown below using the KPSS test, p-value = 0.1, which is greater than the significance level, 0.05, so the alternative hypothesis is accepted and it can be concluded that the differenced time series data is stationary.

## 
##  KPSS Test for Level Stationarity
## 
## data:  diff(ts_df[, "Exports"])
## KPSS Level = 0.092232, Truncation lag parameter = 3, p-value = 0.1

ARIMA(2,1,0) with drift

The results of fitting the proposed ARIMA(2,1,0) model with dift or mean is shown below. The ACF plot in Figure 8.1 of the residuals from the ARIMA(2,1,0) model shows that most autocorrelations are within the threshold limits, which indicates that the residuals are behaving like white noise. However, the ACF has a somewhat significant spike at lag 8. Therefore, this model may not be good when compared to other possible models. In addition, the Ljung-Box test, which is a test for a group of autocorrelations (portmanteau test), returns large values of \(Q^* = 10.7\) and p-value \(=0.1198 \gt 0.05\), also suggesting that the residuals are white noise since we fail to reject the null hypothesis of the test and conclude that the data values are independent. The ARIMA(2,1,0) with non-zero mean can be defined as follows.

\[y_t = c - 0.5230y_{t-1}-0.3065y_{t-2} + w_t\] where \(c=-0.2120[1-(- 0.5230)-(-0.3065)]= -0.387854\) and \(w_t\) is white noise with standard deviation \(\sqrt{6.675}\approx2.58\)

## Series: ts_df[, "Exports"] 
## ARIMA(2,1,0) with drift 
## 
## Coefficients:
##           ar1      ar2    drift
##       -0.5230  -0.3065  -0.2120
## s.e.   0.1262   0.1248   0.1841
## 
## sigma^2 = 6.675:  log likelihood = -133.63
## AIC=275.25   AICc=276.02   BIC=283.43
## 
## Training set error measures:
##                      ME     RMSE      MAE      MPE     MAPE      MASE
## Training set 0.02598873 2.492864 1.844491 -1.00703 8.930269 0.8422792
##                    ACF1
## Training set 0.05853289

Figure 8.1: Residual plots for the ARIMA(2,1,0) model with drift

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,0) with drift
## Q* = 11.459, df = 7, p-value = 0.1198
## 
## Model df: 3.   Total lags used: 10

Using the python package statsmodels to fit the ARIMA(2,1,0) with drift model, the following results are generated. The coefficient of the \(AR_1\) term is -0.5230 and the p-value in P>|z| column is highly significant since it is less than significance level 0.05. Similarly, the coefficient of the \(AR_2\) term, which is -0.3066, is also significant but less significant than \(AR_1\) term with p-value slightly less than \(0.05\). However, the drift or \(X1\) in the table below is insignificant with p-value = 0.309, which is much greater than 0.05.Therefore, the model can be fitted without the mean or drift.

##                                SARIMAX Results                                
## ==============================================================================
## Dep. Variable:                      y   No. Observations:                   58
## Model:                 ARIMA(2, 1, 0)   Log Likelihood                -133.627
## Date:                Fri, 24 Jun 2022   AIC                            275.254
## Time:                        14:00:40   BIC                            283.426
## Sample:                             0   HQIC                           278.430
##                                  - 58                                         
## Covariance Type:                  opg                                         
## ==============================================================================
##                  coef    std err          z      P>|z|      [0.025      0.975]
## ------------------------------------------------------------------------------
## x1            -0.2120      0.209     -1.016      0.309      -0.621       0.197
## ar.L1         -0.5230      0.155     -3.381      0.001      -0.826      -0.220
## ar.L2         -0.3066      0.151     -2.033      0.042      -0.602      -0.011
## sigma2         6.3232      1.086      5.821      0.000       4.194       8.452
## ===================================================================================
## Ljung-Box (L1) (Q):                   0.21   Jarque-Bera (JB):                 5.10
## Prob(Q):                              0.65   Prob(JB):                         0.08
## Heteroskedasticity (H):               0.38   Skew:                             0.57
## Prob(H) (two-sided):                  0.04   Kurtosis:                         3.92
## ===================================================================================
## 
## Warnings:
## [1] Covariance matrix calculated using the outer product of gradients (complex-step).

ARIMA(2,1,0)

Based on the results shown below, fitting the ARIMA(2,1,0) model without the constant lowers the AIC, AICc, and BIC. However, it increased the log-likelihood and standard deviation \(\sigma\). It also increased the errors: RMSE, MAE, MAPE, MASE. Despite this, it removed the spike at lag 8 in the ACF, which is now inside the threshold limit. The residuals are also still normally distributed. The p-values of the two \(AR\) terms are also both significant. Thus, the constant can be excluded from the model. The ARIMA(2,1,0) with zero constant can be defined as follows.

\[y_t = -0.5050y_{t-1}-0.2897y_{t-2} + w_t\]
where \(w_t\) is white noise with standard deviation \(\sqrt{6.706}\approx 2.59\).

## Series: ts_df[, "Exports"] 
## ARIMA(2,1,0) 
## 
## Coefficients:
##           ar1      ar2
##       -0.5050  -0.2897
## s.e.   0.1266   0.1254
## 
## sigma^2 = 6.706:  log likelihood = -134.27
## AIC=274.54   AICc=274.99   BIC=280.67
## 
## Training set error measures:
##                      ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -0.3425155 2.521754 1.867171 -2.968058 9.136747 0.8526358
##                    ACF1
## Training set 0.04169984

Figure 8.2: Residual plots for the ARIMA(2,1,0) model

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,0)
## Q* = 10.698, df = 8, p-value = 0.2194
## 
## Model df: 2.   Total lags used: 10

##                                SARIMAX Results                                
## ==============================================================================
## Dep. Variable:                      y   No. Observations:                   58
## Model:                 ARIMA(2, 1, 0)   Log Likelihood                -134.268
## Date:                Fri, 24 Jun 2022   AIC                            274.537
## Time:                        14:00:41   BIC                            280.666
## Sample:                             0   HQIC                           276.919
##                                  - 58                                         
## Covariance Type:                  opg                                         
## ==============================================================================
##                  coef    std err          z      P>|z|      [0.025      0.975]
## ------------------------------------------------------------------------------
## ar.L1         -0.5050      0.148     -3.413      0.001      -0.795      -0.215
## ar.L2         -0.2897      0.150     -1.937      0.053      -0.583       0.003
## sigma2         6.4708      1.118      5.787      0.000       4.279       8.662
## ===================================================================================
## Ljung-Box (L1) (Q):                   0.09   Jarque-Bera (JB):                 5.04
## Prob(Q):                              0.77   Prob(JB):                         0.08
## Heteroskedasticity (H):               0.44   Skew:                             0.56
## Prob(H) (two-sided):                  0.08   Kurtosis:                         3.92
## ===================================================================================
## 
## Warnings:
## [1] Covariance matrix calculated using the outer product of gradients (complex-step).

In Figure 8.3 below, The two red dots represent the roots of the polynomials \(\phi(B)\). The fitted ARIMA(2,1,0) model is stationary since all the red dots are inside the unit circle. Additonally, the roots are far away from the unit circle, which suggests that they are numerically stable and would be useful for forecasting.

Figure 8.3: Inverse characteristic rootsfor the ARIMA(2,1,0) model

ARIMA(0,1,3) with drift

The results of fitting the proposed ARIMA(0,1,3) model is also shown below. The ACF plot of the residuals in Figure 9.1 from the ARIMA(0,1,3) model shows that all autocorrelations are within the threshold limits, which indicates that the residuals are behaving like white noise. Moreover, the Ljung-Box test returns a large p-value= 0.4368, also suggesting that the residuals are white noise. The log-likelihood, AIC, and BIC are approximately equivalent to the fitted ARIMA(2,1,0) model above. The ARIMA(0,1,3) with non-zero mean can be defined as follows.

\[y_t = c -0.4537 w_{t-1} + 0.0922w_{t-2} + 0.2677w_{t-3} + w_t\]

where \(c=-0.1999\) and \(w_t\) is white noise with standard deviation \(\sqrt{6.611}\approx2.57\)

## Series: ts_df[, "Exports"] 
## ARIMA(0,1,3) with drift 
## 
## Coefficients:
##           ma1     ma2     ma3    drift
##       -0.4537  0.0922  0.2677  -0.1999
## s.e.   0.1319  0.1532  0.1354   0.2946
## 
## sigma^2 = 6.611:  log likelihood = -132.9
## AIC=275.8   AICc=276.97   BIC=286.01
## 
## Training set error measures:
##                         ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -0.0001457588 2.457915 1.819249 -0.838046 8.820526 0.8307526
##                     ACF1
## Training set -0.03570582

Figure 9.1: Residual plots for the ARIMA(0,1,3) model with drift

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,3) with drift
## Q* = 5.88, df = 6, p-value = 0.4368
## 
## Model df: 4.   Total lags used: 10

The following results are generated using statsmodels to fit the ARIMA(0,1,3) model with drift. The coefficient of the \(MA_1\) term is -0.4537 and the p-value in P>|z| column is highly significant since it is less than significance level 0.05. However, the coefficient of the \(MA_2\) term, which is 0.0922, is not significant with p-value = 0.514, which is greater than 0.05. The \(MA_3\) term is also significant with coefficient = 0.2677 and p-value=0.05. However, the drift or \(X1\) in the table below is insignificant with p-value = 0.523, which is much greater than 0.05.Therefore, the model can be fitted with two \(MA\) terms and without the mean or drift.

##                                SARIMAX Results                                
## ==============================================================================
## Dep. Variable:                      y   No. Observations:                   58
## Model:                 ARIMA(0, 1, 3)   Log Likelihood                -132.898
## Date:                Fri, 24 Jun 2022   AIC                            275.797
## Time:                        14:00:42   BIC                            286.012
## Sample:                             0   HQIC                           279.767
##                                  - 58                                         
## Covariance Type:                  opg                                         
## ==============================================================================
##                  coef    std err          z      P>|z|      [0.025      0.975]
## ------------------------------------------------------------------------------
## x1            -0.1999      0.313     -0.639      0.523      -0.813       0.413
## ma.L1         -0.4537      0.132     -3.439      0.001      -0.712      -0.195
## ma.L2          0.0922      0.141      0.653      0.514      -0.185       0.369
## ma.L3          0.2677      0.136      1.962      0.050       0.000       0.535
## sigma2         6.1473      1.134      5.421      0.000       3.925       8.370
## ===================================================================================
## Ljung-Box (L1) (Q):                   0.08   Jarque-Bera (JB):                 2.49
## Prob(Q):                              0.78   Prob(JB):                         0.29
## Heteroskedasticity (H):               0.42   Skew:                             0.41
## Prob(H) (two-sided):                  0.07   Kurtosis:                         3.62
## ===================================================================================
## 
## Warnings:
## [1] Covariance matrix calculated using the outer product of gradients (complex-step).

ARIMA(0,1,2)

Based on the results shown below, fitting the ARIMA(0,1,2) model without the constant lowers the AIC, AICc, and BIC. However, it increased the log-likelihood, standard deviation \(\sigma\), RMSE, MAE, MAPE, and MASE. It also resulted in a spike at lag 8 in the ACF, which is now very close to the threshold limit but still inside it. The residuals are also still normally distributed. The p-values of the two \(MA\) terms are also both significant. Therefore, the constant can be excluded from the model and the number of \(MA\) terms can be reduced to two. The ARIMA(0,1,2) with zero constant can be defined as follows.

\[y_t = -0.5847w_{t-1}+0.2873w_{t-2} + w_t\]

where \(w_t\) is white noise with standard deviation \(\sqrt{6.837}\approx 2.61\).

## Series: ts_df[, "Exports"] 
## ARIMA(0,1,2) 
## 
## Coefficients:
##           ma1     ma2
##       -0.5847  0.2873
## s.e.   0.1419  0.1626
## 
## sigma^2 = 6.837:  log likelihood = -134.85
## AIC=275.7   AICc=276.15   BIC=281.83
## 
## Training set error measures:
##                      ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -0.2728479 2.546162 1.901068 -2.470839 9.502337 0.8681151
##                    ACF1
## Training set 0.04374423

Figure 9.2: Residual plots for the ARIMA(0,1,2) model

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,2)
## Q* = 12.165, df = 8, p-value = 0.144
## 
## Model df: 2.   Total lags used: 10

##                                SARIMAX Results                                
## ==============================================================================
## Dep. Variable:                      y   No. Observations:                   58
## Model:                 ARIMA(0, 1, 2)   Log Likelihood                -134.848
## Date:                Fri, 24 Jun 2022   AIC                            275.697
## Time:                        14:00:43   BIC                            281.826
## Sample:                             0   HQIC                           278.079
##                                  - 58                                         
## Covariance Type:                  opg                                         
## ==============================================================================
##                  coef    std err          z      P>|z|      [0.025      0.975]
## ------------------------------------------------------------------------------
## ma.L1         -0.5847      0.149     -3.936      0.000      -0.876      -0.294
## ma.L2          0.2873      0.126      2.272      0.023       0.039       0.535
## sigma2         6.5966      1.097      6.015      0.000       4.447       8.746
## ===================================================================================
## Ljung-Box (L1) (Q):                   0.10   Jarque-Bera (JB):                 2.12
## Prob(Q):                              0.75   Prob(JB):                         0.35
## Heteroskedasticity (H):               0.51   Skew:                             0.28
## Prob(H) (two-sided):                  0.15   Kurtosis:                         3.76
## ===================================================================================
## 
## Warnings:
## [1] Covariance matrix calculated using the outer product of gradients (complex-step).

The three red dots in Figure 11.2 below correspond to the roots of the polynomials \(\theta(B)\). The fitted ARIMA(0,1,3) model is invertible since all the roots are inside the unit circle.

Figure 9.3: Inverse characteristic rootsfor the ARIMA(0,1,2) model

Auto selection

Using automated model selection, the best model identified is ARIMA(2,1,0) with drift, which is one of the candidate models identified by looking at the ACF and PACF of the differenced time series Exports in the previous section. The results for fitting the ARIMA(2,1,0) model is also shown below, which is similar to the plot shown in Figure 8.1. In figure 10, the ACF plot of the residuals from the ARIMA(2,1,0) model shows that all autocorrelations are within the threshold limits, which indicates that the residuals are behaving like white noise. In addition, the Ljung-Box test returns a large p-value \(= 0.1198 > 0.05\), also suggesting that the residuals are white noise. The log-likelihood, AIC, and BIC are similar to the fitted ARIMA(0,1,3) model above.

## Series: ts_df[, "Exports"] 
## ARIMA(2,1,0) with drift 
## 
## Coefficients:
##           ar1      ar2    drift
##       -0.5230  -0.3065  -0.2120
## s.e.   0.1262   0.1248   0.1841
## 
## sigma^2 = 6.675:  log likelihood = -133.63
## AIC=275.25   AICc=276.02   BIC=283.43
## 
## Training set error measures:
##                      ME     RMSE      MAE      MPE     MAPE      MASE
## Training set 0.02598873 2.492864 1.844491 -1.00703 8.930269 0.8422792
##                    ACF1
## Training set 0.05853289

Figure 10: Residual plots for the autoselected ARIMA(2,1,0) with drift model

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,0) with drift
## Q* = 11.459, df = 7, p-value = 0.1198
## 
## Model df: 3.   Total lags used: 10

Stepwise selection

ARIMA(2,1,2) with drift

Using stepwise selection, the best model identified is ARIMA(2,1,2) with drift. Similar to the ACF plots above, the ACF plot for ARIMA(2,1,2) shows that the residuals are white noise since all the autocorrelations are within the threshold limits (blue dashed line), which is further supported by the large p-value obtained using Ljung-Box test, which is greater than significance level = \(0.05\). The log-likelihood, AIC, and BIC are also similar to the fitted models above.The ARIMA(2,1,2) with non-zero mean can be defined as follows.

\[y_t = c -0.6722y_{t-1}-0.6989y_{t-2}+0.2273 w_{t-1} + 0.4558w_{t-2} + w_t\]

where \(c=-0.2099[1-(-0.6722)-(-0.6989)]=-0.49769389\) and \(w_t\) is white noise with standard deviation \(\sqrt{6.446}\approx2.54\).

## Series: ts_df[, "Exports"] 
## ARIMA(2,1,2) with drift 
## 
## Coefficients:
##           ar1      ar2     ma1     ma2    drift
##       -0.6722  -0.6989  0.2273  0.4558  -0.2099
## s.e.   0.1897   0.2010  0.2570  0.2486   0.2288
## 
## sigma^2 = 6.446:  log likelihood = -131.69
## AIC=275.38   AICc=277.06   BIC=287.64
## 
## Training set error measures:
##                      ME     RMSE      MAE        MPE    MAPE      MASE
## Training set 0.01876558 2.403994 1.796644 -0.8271656 8.75963 0.8204299
##                     ACF1
## Training set -0.02800852

Figure 11.1: Residual plots for the ARIMA(2,1,2) model

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,2) with drift
## Q* = 4.3964, df = 5, p-value = 0.4939
## 
## Model df: 5.   Total lags used: 10

The following results are generated using statsmodels to fit the ARIMA(2,1,2) model with drift. The drift or \(X1\) in the table below is insignificant with p-value=0.409, which is much greater than 0.05. Thus, it can be excluded from the model. The coefficient of the \(AR_1\) term is -0.6721 with p-value=0.002 so it is significant. The \(AR_2\) term is also significant. Both the \(MA\) terms appear to be insignificant with p-values greater than 0.05. Therefore, the model can be fitted with two \(MA\) terms and without the mean or drift. Based on these results, an ARIMA model without drift and without \(MA\) terms can be fitted. This model would be an ARIMA(2,1,0) without drift, which is fitted in the previous section.

##                                SARIMAX Results                                
## ==============================================================================
## Dep. Variable:                      y   No. Observations:                   58
## Model:                 ARIMA(2, 1, 2)   Log Likelihood                -131.689
## Date:                Fri, 24 Jun 2022   AIC                            275.378
## Time:                        14:00:46   BIC                            287.636
## Sample:                             0   HQIC                           280.142
##                                  - 58                                         
## Covariance Type:                  opg                                         
## ==============================================================================
##                  coef    std err          z      P>|z|      [0.025      0.975]
## ------------------------------------------------------------------------------
## x1            -0.2101      0.254     -0.826      0.409      -0.708       0.288
## ar.L1         -0.6721      0.215     -3.123      0.002      -1.094      -0.250
## ar.L2         -0.6989      0.204     -3.434      0.001      -1.098      -0.300
## ma.L1          0.2273      0.227      1.001      0.317      -0.218       0.672
## ma.L2          0.4558      0.324      1.406      0.160      -0.180       1.091
## sigma2         5.8805      0.986      5.961      0.000       3.947       7.814
## ===================================================================================
## Ljung-Box (L1) (Q):                   0.05   Jarque-Bera (JB):                 6.71
## Prob(Q):                              0.83   Prob(JB):                         0.03
## Heteroskedasticity (H):               0.32   Skew:                             0.64
## Prob(H) (two-sided):                  0.02   Kurtosis:                         4.10
## ===================================================================================
## 
## Warnings:
## [1] Covariance matrix calculated using the outer product of gradients (complex-step).

The two red dots in the left and right hand plot of Figure 11.2 below correspond to the roots of the polynomials \(\phi(B)\) and \(\theta(B)\), respectively. As expected, all the red dots are inside the unit circle. Hence, the fitted model is both stationary and invertible.

Figure 11.2: Inverse characteristic rootsfor the ARIMA(2,1,2) model

ARIMA(2,1,2)

Using stepwise selection without including the mean, the best model identified is ARIMA(2,1,2) without drift. Similar to the ACF plots above, the ACF plot for ARIMA(2,1,2) shows that the normally-distributed residuals are white noise since all the autocorrelations are within the threshold limits, which is further supported by the large p-value obtained using Ljung-Box test, which is greater than significance level = \(0.05\). The log-likelihood, standard deviation, AIC, AICc, and BIC are lower than the fitted ARIMA(2,1,2) model with drift above.The ARIMA(2,1,2) with zero mean can be defined as follows.

\[y_t = -0.6741y_{t-1} -0.7142y_{t-2}+0.2468w_{t-1} + 0.4831w_{t-2} + w_t\] where \(w_t\) is white noise with standard deviation \(\sqrt{6.416}\approx2.53\).

## Series: ts_df[, "Exports"] 
## ARIMA(2,1,2) 
## 
## Coefficients:
##           ar1      ar2     ma1     ma2
##       -0.6741  -0.7142  0.2468  0.4831
## s.e.   0.1821   0.2037  0.2531  0.2576
## 
## sigma^2 = 6.416:  log likelihood = -132.1
## AIC=274.2   AICc=275.37   BIC=284.41
## 
## Training set error measures:
##                      ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -0.2645835 2.421275 1.821659 -2.327502 8.929045 0.8318532
##                     ACF1
## Training set -0.04168949

Figure 11.1: Residual plots for the ARIMA(2,1,2) model

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,2)
## Q* = 4.122, df = 6, p-value = 0.6602
## 
## Model df: 4.   Total lags used: 10

Using a different package (fable) to choose the best model by stepwise selection, the same model is obtained as shown below.

## Series: value 
## Model: ARIMA(2,1,2) 
## 
## Coefficients:
##           ar1      ar2     ma1     ma2
##       -0.6741  -0.7142  0.2468  0.4831
## s.e.   0.1821   0.2037  0.2531  0.2576
## 
## sigma^2 estimated as 6.416:  log likelihood=-132.1
## AIC=274.2   AICc=275.37   BIC=284.41

Fitting the ARIMA(2,1,2) model without drift below, shows that both \(MA\) terms are still insignificant with p-values greater than 0.05, which further supports that the \(MA\) terms can be excluded and ARIMA(2,1,0) is an ideal model.

##                                SARIMAX Results                                
## ==============================================================================
## Dep. Variable:                      y   No. Observations:                   58
## Model:                 ARIMA(2, 1, 2)   Log Likelihood                -132.098
## Date:                Fri, 24 Jun 2022   AIC                            274.197
## Time:                        14:00:48   BIC                            284.412
## Sample:                             0   HQIC                           278.167
##                                  - 58                                         
## Covariance Type:                  opg                                         
## ==============================================================================
##                  coef    std err          z      P>|z|      [0.025      0.975]
## ------------------------------------------------------------------------------
## ar.L1         -0.6740      0.202     -3.329      0.001      -1.071      -0.277
## ar.L2         -0.7141      0.196     -3.637      0.000      -1.099      -0.329
## ma.L1          0.2468      0.216      1.141      0.254      -0.177       0.671
## ma.L2          0.4831      0.315      1.531      0.126      -0.135       1.101
## sigma2         5.9653      0.994      6.003      0.000       4.018       7.913
## ===================================================================================
## Ljung-Box (L1) (Q):                   0.12   Jarque-Bera (JB):                 6.72
## Prob(Q):                              0.73   Prob(JB):                         0.03
## Heteroskedasticity (H):               0.35   Skew:                             0.63
## Prob(H) (two-sided):                  0.03   Kurtosis:                         4.10
## ===================================================================================
## 
## Warnings:
## [1] Covariance matrix calculated using the outer product of gradients (complex-step).

Model comparisons

Using the errors in the table below, the best model can be identified based on accuracy. It appears that ARIMA(2,1,2) with drift is the best model with generally the largest log-likelihood and smallest errors (RMSE, MAE, MAPE,MASE). ARIMA(2,1,2) without drift has the lowest standard deviation and AIC. However, the constant and \(MA\) terms in these two models are not statistically significant, which reduces them to an ARIMA(2,1,0) model. Thus, ARIMA(2,1,0) appears to be the best model with lower AICc and BIC compared to the other models.

Table 1: Model Comparisons