Use of Teleconnections to Predict Western Australian Seasonal Rainfall Using ARIMAX Model

: Increased demand for engineering propositions to forecast rainfall events in an area or region has resulted in developing di ﬀ erent rainfall prediction models. Interestingly, rainfall is a very complicated natural system that requires consideration of various attributes. However, regardless of the predictability performance, easy to use models have always been welcomed over the complex and ambiguous alternatives. This study presents the development of Auto–Regressive Integrated Moving Average models with exogenous input (ARIMAX) to forecast autumn rainfall in the South West Division (SWD) of Western Australia (WA). Climate drivers such as Indian Ocean Dipole (IOD) and El Nino Southern Oscillation (ENSO) were used as predictors. Eight rainfall stations with 100 years of continuous data from two coastal regions (south coast and north coast) were selected. In the south coast region, Albany (0,1,1) with exogenous input DMI Oct –Nino3 Nov , and Northampton (0,1,1) with exogenous input DMI Jan –Nino3 Nov were able to forecast autumn rainfall 4 months and 2 months in advance, respectively. Statistical performance of the ARIMAX model was compared with the multiple linear regression (MLR) model, where for calibration and validation periods, the ARIMAX model showed signiﬁcantly higher correlations (0.60 and 0.80, respectively), compared to the MLR model (0.44 and 0.49, respectively). It was evident that the ARIMAX model can predict rainfall up to 4 months in advance, while the MLR has shown strict limitation of prediction up to 1 month in advance. For WA, the developed ARIMAX model can help to overcome the di ﬃ culty in seasonal rainfall prediction as well as its application can make an invaluable contribution to stakeholders’ economic preparedness plans.


Introduction
Long-term rainfall forecasting is one of the most challenging and demanding tasks. However, this is an important issue, which is directly related to the economy of a country. A reliable and accurate forecast can be beneficial for different authorities dealing with disaster management, water management, agricultural production, and flood management to set strategies, decision making, and taking precautionary measures before any natural calamities such as floods, drought, and bushfire could occur. Generally, two methods (statistical and dynamic modeling) have been used to forecast long-term rainfall [1]. A statistical method is less complex and requires less development time, but an uninterrupted and reliable data source is mandatory. On the other hand, the dynamic method is quite complex, requires more development time and expense [2,3]. Therefore, the application of statistical methods to forecast rainfall has become very popular among the researchers due to its simplicity, cost-effectiveness, and easy to implement characteristics. However, the use of the physically-based empirical model for seasonal rainfall prediction has also become popular since it overcomes the difficulties associated with conventional dynamic and statistical models [4,5]. To conduct this study, monthly rainfall data and climate index data were collected for eight different stations located in two different regions. The monthly rainfall data was retrieved from the Australian Bureau of Meteorology website (www.bom.gov.au/climate/data/) and climate index data was collected from the climate explorer website (http://climexp.knmi.nl/). All the rainfall stations were in the south-west division of WA and for each station, 100-years continuous rainfall data was collected. The location, description, and seasonal rainfall data for Summer (Dec-Jan-Feb), Autumn (Mar-Apr-May), Winter (Jun-Jul-Aug), and Spring (Sep-Oct-Nov) for these stations is presented in Table 1. Past studies suggest that climate indices namely, Dipole Mode Index (DMI), Nino3, Nino3.4, Nino4, Southern Oscillation Index (SOI), Southwest Australian Circulation (SWAC), Southern Annular Mode (SAM), and ENSO Modoki Index (EMI) cause influence on Western Australian rainfall To conduct this study, monthly rainfall data and climate index data were collected for eight different stations located in two different regions. The monthly rainfall data was retrieved from the Australian Bureau of Meteorology website (www.bom.gov.au/climate/data/) and climate index data was collected from the climate explorer website (http://climexp.knmi.nl/). All the rainfall stations were in the south-west division of WA and for each station, 100-years continuous rainfall data was collected. The location, description, and seasonal rainfall data for Summer (Dec-Jan-Feb), Autumn (Mar-Apr-May), Winter (Jun-Jul-Aug), and Spring (Sep-Oct-Nov) for these stations is presented in Table 1.  [26,33,36,59]. For this study, Nino3, Nino3.4, Nino4, SOI, EMI, and DMI were selected as potential climate indices due to the availability of 100 years of continuous data. However, SWAC and SAM were excluded from the analyses due to the non-availability of long-term continuous data.
All the climate indices have a complex but interesting formation mechanism. ENSO is measured by two types of indicators; SOI, which is the measure of Sea Level Pressure (SLP) difference between Darwin and Tahiti; and the rest are Nino3.4, Nino3, and Nino4. Nino3.4 is the average Sea Surface Temperature (SST) anomalies in the western pacific bounded by 5 • N to 5 • S, from 170 • W to 120 • W; Nino3 is the average SST anomalies in the eastern pacific bounded by 5 • S-5 • N, from 90 • W-150 • W; and Nino4 is the average SST anomalies in central pacific bounded by 5 • N to 5 • S, from 150 • W to 160 • E [26]. The EMI is the difference between El-Nino and traditional events, which have maximum warming in central Pacific and eastern Pacific, respectively. It is expressed by the following equation: [60]. DMI is the indicator of IOD, which is defined as the average SST anomalies between tropical western Indian Ocean bounded by 10 • N to 10 • S, from 50 • E to 70 • E and the tropical south-eastern Indian Ocean bounded by 110 • S to the equator from 90 • E to 110 • E [40].

Methodology
In this study, Time Series Analysis was used to develop forecasting models. Among various time series methods, ARIMA method has become very popular. For time series analysis, two types of ARIMA methods are mostly used: univariate ARIMA and multivariate ARIMA or transfer function model. In an ARIMA model, the inclusion of exogenous variable (predictors) or explanatory variables is termed as ARIMAX model [15]. In the ARIMAX model, ARIMA orders for dependent variables while transfer function orders for independent variables or predictors. Using ARIMAX model in particular, Kamruzzaman, Beecham [12] have illustrated a detailed way of including predictors, the possible number of predictors, and how the inclusion of the past value of climate indices can improve the forecasting ability. In this study, the ARIMAX model was used for the prediction of Western Australian autumn rainfall, where large-scale climate indices were used as independent variables. 100 years of data were split into two periods: the calibration period  and the validation period . Collected monthly rainfall data was then converted into seasonal autumn rainfall (Australian standard autumn months: March-April-May) data and climate indices data was collected and accumulated for the preceding months to perform analysis with lagged climate indices. Bivariate correlation was performed to select significant independent variables which were later included as transfer functions to develop the ARIMAX model. Results of the ARIMAX model were evaluated using standard statistical test parameters, profoundly used in past studies [46,47,[61][62][63]. These include Pearson correlation coefficient (r), root mean square error (RMSE), mean absolute percentage error (MAPE), mean absolute error (MAE), normalized Bayesian information criterion (BIC), and refined Willmot index of agreement (d r ) values. Additionally, ARIMAX produced results were compared with the results of multiple linear regression (MLR) analysis which were performed earlier in different investigations [25,61]. All these analyses were performed in the IBM SPSS Statistics 26 software package.

Multivariate Auto-Regressive Integrated Moving Average with Exogenous Input (ARIMAX)
ARIMA model is composed of 'AR', 'I', and 'MA' where 'AR' stands for Auto-Regressive, 'I' is for Integrated, a time series which needs to be differenced to make a non-stationary series to stationary, 'MA' stands for Moving Average. The general expression of the ARIMA model is ARIMA (p, d, q) *(P, D, Q). This contains two parts: one is non-seasonal and the other one is seasonal, where, (p, d, q) is the non-seasonal part, and (P, D, Q) is the seasonal part. Here, 'p' denotes non-seasonal auto-regressive order, 'd' denotes non-seasonal differencing, 'q' denotes non-seasonal moving average, whereas, 'P' denotes seasonal auto-regressive order, 'D' denotes seasonal differencing, and 'Q' denotes seasonal moving average. In this study, the time-series show no seasonality, therefore, only the non-seasonal part was used. In the ARIMAX model, two types of input orders are required: one is ARIMA order (dependent variable; autumn rainfall) and the other is transfer function order (predictors or independent variable; climate indices) [64]. A description of these input orders are presented as follows: ARIMA Order (p, d, q): • Autoregressive (p): it is defined as the number of autoregressive orders in the model. Autoregressive orders specify which previous values from the series are used to predict the present values. • Difference (d): it is defined as the order of differencing in the model. Differencing is necessary to make a non-stationary series into stationary series. As the ARIMA model is stationary, it is necessary to remove the non-stationary effects before estimating models. A first-order differencing is often proved as enough for linear trends, while a second-order differencing is required for quadratic trend.

•
Moving Average (q): it is defined as the number of moving average orders in the model. Moving average orders specify how deviation from the series mean of previous values (past errors) is used to predict current values. Therefore, it is the number of lagged forecast errors in the predicted equation.
Therefore, for example, a model described as ARIMA (0,1,1) refers to containing zero auto-regressive (p) and one moving average (q) parameters with the first order of differencing (d) of the data.
Transfer Function Orders: Transfer function order specifies which past values of independent variables are used to predict future values of dependent series.

•
Numerator: it is the number of orders of the numerator in the transfer function. It states which past values of independent series are used to predict the current value of the dependent series. Such as a numerator value of 1 states that the one-time period past value and the current value of an independent series are used to predict the current values of the dependent series.

•
Denominator: it is the number of orders of the denominator in the transfer function. It is specified to predict the current value of the dependent series and how deviations from the series mean for previous values of independent variables are used. For instance, a denominator number of 1 indicates that the deviation from the mean of the independent series one period in the past is considered to predict the current values of the dependent series. • Difference: it is the order of differencing. It is applied to make a non-stationary series into stationary before estimating the model.
An elaborative formation of the ARIMAX model equation has been discussed and presented in Box, Jenkins [65], and Hamilton [66]. The ARIMAX (p, d, q) model equation for time series Y t and exogeneous data X t is presented below in Equation (1): where, ϕ 1 , . . . , ϕ p and θ 1 , . . . , θ q are the parameters; ε t , ε t−1 are white noise error and β 1 , . . . , β m are the parameters of independent variables input X t and t is the time.

1.
Identification: in this stage, first, the raw data series is plotted to identify whether the data is stationary or not. If the raw data series is found to be non-stationary, differencing is required. After the first order differencing, correlograms of the autocorrelation function (ACF) and partial autocorrelation function (PACF) is investigated. From these plots, the order of AR and MA gets identified.

2.
Parameter Estimation and Selection: the number of AR depends on the lag of PACF cuts and the number of MA depends on the lag of the ACF plot. However, decision making on the order of AR and MA by looking at the cuts/spikes is not straightforward. Most of the time it required experimentation with several alternative orders of different models to choose the appropriate order. The following guidelines are usually followed during the selection of the AR and MA order: • If the ACF plot shows exponential decay and PACF spikes at lag-1, no correlation for other lags, in that case, one autoregressive parameter (p = 1) can be selected.

•
If the ACF plot shows a sine-wave shape pattern or a set of exponential decay and PACF spikes at lag-1 and lag-2, no correlation for other lags, in that case, two autoregressive parameters (p = 2) can be selected.

•
If the PACF plot shows exponential decay and ACF spikes at lag-1, no correlation for other lags, in that case, one moving average parameter (q = 1) can be selected.

•
If the PACF plot shows a sine-wave shape pattern or a set of exponential decay and ACF spikes at lag-1 and lag-2, no correlation for other lags, in that case, two moving average parameters (q = 2) can be selected.

•
One auto-regressive and one moving average parameter can be selected if both shows exponential decay starting at lag-1.

•
Sometimes, using both AR and MA orders in a model can cancel each other's impact. Therefore, it is often wise to use mixed AR and MA models with less number of orders.

3.
Diagnostics Check: the diagnostic check is required to verify the adequacy of the developed model. The residual of the developed model should be white noise (no autocorrelation). To check whether the residual is white noise or not, at first, an inspection of the residual ACF and PACF plot is required. If 95% of the spikes stay between the black lines, it indicates that the autocorrelation is white noise. If two or more spikes or more than 5% of spikes are located outside of the boundary line, then the series is not white noise. Another way of checking the model accuracy is to perform the Ljung-Box test. Such a test is conducted to verify the null hypothesis of being white noise of residual if the p-value is greater than 0.05 [68]. A p-value greater than 0.05 implies that lag autocorrelation among the residuals is zero and the developed model is adequate to fit the data set.

Statistical Parameters for ARIMAX Model
Among all the statistical parameters that were used to evaluate the ARIMAX model's performance, RMSE calculates prediction errors, the measure of how much a dependent series varies from its model-predicted level. MAE is the average of the absolute errors/residuals between observed and predicted value, while MAPE is the measure of prediction accuracy of a developed model. It is also known as the mean absolute percentage deviation. For both RMSE and MAE, a value of 0 indicates a perfect predictability performance. Thus, the lower the value of RMSE, MAE, and MAPE the better is the model's performance [69,70]. The equation for RMSE, MAE, and MAPE is presented in Equations (2)-(4): where O i is the observed value, P i is the predicted value and n is the number of observations. Normalized Bayesian information criterion (BIC) is the measure of the overall fit of the developed model. It measures the model's complexity based on mean square error, the number of parameters, and series length. A lower value indicates a low complexity of the model and reverses for the higher value. Thus, a model with a lower BIC is always preferred.
However, all these tests have some shortcomings, which can be overruled by the introduction of a refined index of agreement (d r ) developed by Willmott [71]. 'd r ' is the reformulation of Willmott's index of agreement (d). It specifies the sum of the magnitudes of the differences between the predicted and observed deviations from the observed mean relative to the sum of the magnitudes of the perfect model (P i = O i , for all i) and observed deviations from the observed mean. The refined index of agreement (d r ) can be calculated using the following Equation (5): where P i is the predicted value of the ith observation, O i is the observed value of the ith observation, O is the observed mean value, and n is the sample size. c value equals 2 as suggested in the equation.

Results
Eight rainfall stations from two coastal regions, namely, South Coast and North Coast in WA were selected to conduct this study. Lagged monthly climate indices values (June n−1 to Feb n ) were used to investigate the correlation with autumn rainfall (March-April-May). Here, 'n' is the predicted year, and 'n -1' is the previous year. At first, Bivariate correlation analysis was conducted to determine the significant correlation between rainfall and climate indices. Table 2 presents the bivariate correlation results for the selected rainfall stations.
Climate indices, which showed a significant correlation (r) were used as a predictor in the ARIMAX technique. Past studies have shown that a combination of two or more climate indices has improved forecasting ability [26,72]. Therefore, from earlier correlation analysis, a combination of influential climate indices was used as exogenous variables (Predictors) in ARIMAX modeling.
Once significant climate indices were selected, stationarity and seasonality of the data were assessed. In the identification stage, both rainfall and climate index data were found as non-stationary and rainfall pattern was found as non-seasonal. To make the data stationary, a first-order differencing was performed on both rainfall and climate indices data. Differencing was kept limited until first-order only as it deemed sufficient for any linear trend [64]. For the station Albany, Figure 2 depicts the non-stationary and stationary state of rainfall data before and after first order differencing, while Figures 3 and 4 illustrate the same for associated climate indices. A similar strategy has been followed for other stations analyzed in this study.
* Correlation is significant at the 0.05 level (two-tailed), ** Correlation is significant at the 0.01 level (two-tailed).
Hydrology 2020, 7, x FOR PEER REVIEW 9 of 25 * Correlation is significant at the 0.05 level (two-tailed), ** Correlation is significant at the 0.01 level (two-tailed).
(a) (b)    The next step involved parameter estimation, where AR and MA order was selected from the PACF and ACF plots of selected rainfall stations. PACF plot was used to select AR order, while the ACF plot was used to select MA order. Figure 5 presents ACF and PACF plots with respective lag numbers for the station Albany. From Figure 5a, it can be observed that the ACF plot contains only one spike located outside of the confidence limit, while the rest are close to zero. On the other hand, in Figure 5b, in the PACF plot, exponential decay is evident with four spikes being present outside the confidence limit. Considering such facts, ARIMAX (0,1,1) was selected for Albany. For the remaining rainfall stations, a similar strategy was followed to select the appropriate AR and MA orders. For most of the rainfall stations, ARIMAX (0,1,1) model was found as most suitable except two stations (Mount Barker and Busselton Shire). ARIMA orders and transfer function orders for all these stations are presented in Tables 3 and 4. The insertion of ARIMAX criteria for model and transfer functions is presented in Appendix A, Figure A1.  The next step involved parameter estimation, where AR and MA order was selected from the PACF and ACF plots of selected rainfall stations. PACF plot was used to select AR order, while the ACF plot was used to select MA order. Figure 5 presents ACF and PACF plots with respective lag numbers for the station Albany. From Figure 5a, it can be observed that the ACF plot contains only one spike located outside of the confidence limit, while the rest are close to zero. On the other hand, in Figure 5b, in the PACF plot, exponential decay is evident with four spikes being present outside the confidence limit. Considering such facts, ARIMAX (0,1,1) was selected for Albany. For the remaining rainfall stations, a similar strategy was followed to select the appropriate AR and MA orders. For most of the rainfall stations, ARIMAX (0,1,1) model was found as most suitable except two stations (Mount Barker and Busselton Shire). ARIMA orders and transfer function orders for all these stations are presented in Tables 3 and 4. The insertion of ARIMAX criteria for model and transfer functions is presented in Appendix A, Figure A1. The next step involved parameter estimation, where AR and MA order was selected from the PACF and ACF plots of selected rainfall stations. PACF plot was used to select AR order, while the ACF plot was used to select MA order. Figure 5 presents ACF and PACF plots with respective lag numbers for the station Albany.  The next step involved parameter estimation, where AR and MA order was selected from the PACF and ACF plots of selected rainfall stations. PACF plot was used to select AR order, while the ACF plot was used to select MA order. Figure 5 presents ACF and PACF plots with respective lag numbers for the station Albany. From Figure 5a, it can be observed that the ACF plot contains only one spike located outside of the confidence limit, while the rest are close to zero. On the other hand, in Figure 5b, in the PACF plot, exponential decay is evident with four spikes being present outside the confidence limit. Considering such facts, ARIMAX (0,1,1) was selected for Albany. For the remaining rainfall stations, a similar strategy was followed to select the appropriate AR and MA orders. For most of the rainfall stations, ARIMAX (0,1,1) model was found as most suitable except two stations (Mount Barker and Busselton Shire). ARIMA orders and transfer function orders for all these stations are presented in Tables 3 and 4. The insertion of ARIMAX criteria for model and transfer functions is presented in Appendix A, Figure A1. From Figure 5a, it can be observed that the ACF plot contains only one spike located outside of the confidence limit, while the rest are close to zero. On the other hand, in Figure 5b, in the PACF plot, exponential decay is evident with four spikes being present outside the confidence limit. Considering such facts, ARIMAX (0,1,1) was selected for Albany. For the remaining rainfall stations, a similar strategy was followed to select the appropriate AR and MA orders. For most of the rainfall stations, ARIMAX (0,1,1) model was found as most suitable except two stations (Mount Barker and Busselton Shire). ARIMA orders and transfer function orders for all these stations are presented in Tables 3 and 4. The insertion of ARIMAX criteria for model and transfer functions is presented in Appendix A, Figure A1. Several ARIMAX models were developed using different combinations such as DMI-Nino3, DMI-Nino4, DMI-Nino3.4, DMI-SOI, and DMI-EMI for both the south coast and north coast regions. Table 5 presents the total outcome of a combination set of different models.  Table 5 illustrates that the DMI-Nino3 model shows the highest significant correlations (r) for all rainfall stations in both regions. For Albany, the DMI-SOI model showed the highest correlation (r), however, the DMI-Nino3 model showed a consistently better correlation for all other stations.
Thus, the DMI-Nino3 model was selected as the best model. Table 6 presents the statistical performance of the developed DMI-Nino3 models for selected rainfall stations. Once the ARIMAX model got developed, a diagnostic check was performed to check the adequacy of the model. A Ljung-Box test was conducted to check the residuals, whether these are white noise or not. For all these developed models, the associated p-values are found to be greater than 0.05, which holds the null hypothesis of being white noise [68]. Another alternative way of checking the residual autocorrelation is to draw residual ACF and PACF plots as presented in Figure 6. Figure 6 demonstrates that all the spikes are located within the black boundary lines, depicting that no autocorrelation among the residuals exist. Similar results were found for other rainfall stations (refer to Appendix A, Figure A2).
(r), however, the DMI-Nino3 model showed a consistently better correlation for all other stations. Thus, the DMI-Nino3 model was selected as the best model. Table 6 presents the statistical performance of the developed DMI-Nino3 models for selected rainfall stations. Once the ARIMAX model got developed, a diagnostic check was performed to check the adequacy of the model. A Ljung-Box test was conducted to check the residuals, whether these are white noise or not. For all these developed models, the associated p-values are found to be greater than 0.05, which holds the null hypothesis of being white noise [68]. Another alternative way of checking the residual autocorrelation is to draw residual ACF and PACF plots as presented in Figure  6. Figure 6 demonstrates that all the spikes are located within the black boundary lines, depicting that no autocorrelation among the residuals exist. Similar results were found for other rainfall stations (refer to Appendix A, Figure A2). A validation test was performed with the same model input sets using ARIMAX analysis. In the validation period, the developed model showed an increased correlation (r) compared to the calibration period for all stations except Mount Barker, Northampton, and Nabawa. Refined Willmott A validation test was performed with the same model input sets using ARIMAX analysis. In the validation period, the developed model showed an increased correlation (r) compared to the calibration period for all stations except Mount Barker, Northampton, and Nabawa. Refined Willmott index of agreement (d r ) was also calculated for both calibration and validation periods. Based on the correlation and refined index of agreement statistics, a detailed comparison between the developed ARIMAX model and previously developed MLR models were made and presented in Table 7. From Table 7, it is observed that the ARIMAX model showed significantly higher correlations (0.56-0.82), compared to the MLR model (0.34-0.44) in the calibration period. For Albany and Northampton, the prediction performance of the ARIMAX and MLR model is presented in Figure 7. Similar plots for other stations are also presented in Appendix B, Figure A3. index of agreement ( ) was also calculated for both calibration and validation periods. Based on the correlation and refined index of agreement statistics, a detailed comparison between the developed ARIMAX model and previously developed MLR models were made and presented in Table 7. From Table 7, it is observed that the ARIMAX model showed significantly higher correlations (0.56-0.82), compared to the MLR model (0.34-0.44) in the calibration period. For Albany and Northampton, the prediction performance of the ARIMAX and MLR model is presented in Figure 7. Similar plots for other stations are also presented in Appendix B, Figure A3.  To evaluate the ARIMAX model's performance to capture extreme cases, peak and trough analysis was performed. The correlation coefficients (r) between observed and predicted peak and trough values for selected rainfall stations are presented in Table 8. The DMI-Nino3 model has shown its capability of capturing peaks with a correlation of 0.52-0.90, and 0.77-0.89 on the south coast and north coast, respectively. In contrast, for the trough, the correlation ranged from 0.43 to 0.68, and 0.51 to 0.66, respectively. Peak and trough plots for the station Albany and Northampton are presented in Figure 8, while similar plots for other stations are presented in Appendix B, Figure A4.  To evaluate the ARIMAX model's performance to capture extreme cases, peak and trough analysis was performed. The correlation coefficients (r) between observed and predicted peak and trough values for selected rainfall stations are presented in Table 8. The DMI-Nino3 model has shown its capability of capturing peaks with a correlation of 0.52-0.90, and 0.77-0.89 on the south coast and north coast, respectively. In contrast, for the trough, the correlation ranged from 0.43 to 0.68, and 0.51 to 0.66, respectively. Peak and trough plots for the station Albany and Northampton are presented in Figure 8, while similar plots for other stations are presented in Appendix B, Figure A4.
Finally, a comparative evaluation of previously developed models for WA was performed. A summary of the comparison is presented in Table 9.  Finally, a comparative evaluation of previously developed models for WA was performed. A summary of the comparison is presented in Table 9. The comparison made above (refer to Table 9) suggests that high Pearson correlation (r) values were obtained in the ARIMAX model, referring to its superiority over multiple linear regression (MLR) models. ARIMAX model was also found as more favorable due to its 4 months lagged prediction capability.

Discussion
This paper presents the inclusion of exogenous variables in the ARIMA model (termed as ARIMAX) and showed good prediction performance for WA rainfall variability. The inclusion of an exogenous variable is only possible if these predictors show a significant correlation with the dependent variable [64]. Correlation analyses were conducted between single climate index and autumn rainfall to identify statistically significant climate index for selected rainfall stations in WA. Except for EMI, all the selected climate indices showed low to medium Pearson correlation (r) with autumn rainfall in 5 months lagged period (October to February). However, EMI showed a significant correlation only for north coast rainfall stations that aligns with the findings of Taschetto and England [59]. SST based climate indices (Nino3.4, Nino3, and Nino4) showed a better correlation (r) compared to SLP based index (SOI). This has also been consistent with the findings of Drosdowsky and Chambers [29], where they showed SOI has less predictability skill for Western Australian autumn rainfall compared to SST indices with 1-3 months lag period.
Several ARIMAX models were developed using a combination of different significant climate indices. Model sets were developed with lagged DMI-Nino3, lagged DMI-Nino4, lagged DMI-Nino3.4, lagged DMI-SOI, and lagged DMI-EMI for both the south coast and north coast regions. All the model sets and their Pearson correlation (r) are presented in Table 5, from where the DMI-Nino3 model was selected as the best model for both the regions.  The comparison made above (refer to Table 9) suggests that high Pearson correlation (r) values were obtained in the ARIMAX model, referring to its superiority over multiple linear regression (MLR) models. ARIMAX model was also found as more favorable due to its 4 months lagged prediction capability.

Discussion
This paper presents the inclusion of exogenous variables in the ARIMA model (termed as ARIMAX) and showed good prediction performance for WA rainfall variability. The inclusion of an exogenous variable is only possible if these predictors show a significant correlation with the dependent variable [64]. Correlation analyses were conducted between single climate index and autumn rainfall to identify statistically significant climate index for selected rainfall stations in WA. Except for EMI, all the selected climate indices showed low to medium Pearson correlation (r) with autumn rainfall in 5 months lagged period (October to February). However, EMI showed a significant correlation only for north coast rainfall stations that aligns with the findings of Taschetto and England [59]. SST based climate indices (Nino3.4, Nino3, and Nino4) showed a better correlation (r) compared to SLP based index (SOI). This has also been consistent with the findings of Drosdowsky and Chambers [29], where they showed SOI has less predictability skill for Western Australian autumn rainfall compared to SST indices with 1-3 months lag period.
Several ARIMAX models were developed using a combination of different significant climate indices. Model sets were developed with lagged DMI-Nino3, lagged DMI-Nino4, lagged DMI-Nino3.4, lagged DMI-SOI, and lagged DMI-EMI for both the south coast and north coast regions. All the model sets and their Pearson correlation (r) are presented in Table 5, from where the DMI-Nino3 model was selected as the best model for both the regions.
The statistical performance parameters for the selected DMI-Nino3 model are presented in Table 6, wherein the calibration period, the Pearson correlation (r) for the rainfall magnitudes ranged from a minimum of 0.58 to a maximum of 0.67 in the south coast region. These developed models can predict seasonal autumn rainfall for 4 months in advance for the region. Similarly, for the north coast region, Pearson correlation (r) ranged from 0.56 to 0.82 in the calibration period. The station Mingenew and Northampton predicted seasonal autumn rainfall for 4 months and 2 months in advance. The remaining two stations, namely Nabawa and Ogilvie, showed the prediction for 1 month in advance only. Such capability of the model to predict up to 6 months in advance has also been justified in several past studies [3,24,46,47,61,73]. All these studies considered lagged climate indices to forecast Australian seasonal rainfall as the current study did. Schepen, Wang [3] conducted a study on the evidence of using lagged climate indices for forecasting Australian seasonal rainfall using one climate index. It has shown positive evidence to use climate indices to predict rainfall in several months advance (1-3 months or one season). A study conducted by Abbot and Marohasy [47] used artificial neural network (ANN) with 1-3 months lagged period input parameters. The finding of their study depicted that the monthly forecast in 3 months in advance is as skillful as 1 month in advance. Similar attempts made by Mekanik, Imteaz [46], Ghamariadyan, Imteaz [73], Rasel, Imteaz [35] showed that seasonal rainfall forecast can be possible from 3 to 6 months in advance for different regions in Australia.
The developed ARIMAX models for both south coast and north coast regions have also shown low values of RMSE, MAPE, MAE, and normalized BIC values (refer to Table 6). Low values of all these parameters indicate a good prediction performance of the developed models. Once the models were developed for the calibration period, validation tests were performed with the same model input sets using ARIMAX analysis (refer to Table 7). In the validation period, the developed model showed an increased correlation (r) compared to the calibration period for all stations except Mount Barker, Northampton, and Nabawa. For the south coast region, Albany showed the highest correlation (r) value of 0.60 in calibration and 0.80 in the validation period with a lag of 4 months. In contrast, for the north coast region, Northampton showed the highest Pearson correlation (r) value of 0.82 in calibration and 0.70 in the validation period with a lag of 2 months. Statistically, a Pearson correlation (r) value greater than 0.5 indicates a large effect [74]. A similar observation was made for Refined Willmott index of agreement (d r ), where the 'd r ' value was over 0.60 for all the stations, except for Mingenew. In both calibration and validation periods, the highest value of 'd r ' was reported for Albany on the south coast and Northampton on the north coast, respectively. Since all the stations depicted positive higher values of 'd r ', this indicates a good fit of the model [71].
To understand the effectiveness of the ARIMAX model, its statistical performance parameters were compared with previously developed MLR models for the same regions. From such comparison, a significant rise in Pearson correlation (r) values was observed in ARIMAX models for both calibration and validation periods. An increase in 'd r ' values was also obtained for the ARIMAX models. Moreover, MLR models showed relatively poor prediction performance with only 1 month lag, whereas the ARIMAX models showed better prediction with 4 months lag. In terms of reliability, ARIMAX outperformed its MLR counterparts with better efficiency and accuracy. Furthermore, ARIMAX models successfully captured some of the extreme events, and followed the same rainfall trend as observed rainfall, whereas, the MLR models failed to do so (refer to Figure 7). The finding of this study reinforced ARIMAX models' superiority over MLR models. Many past studies also emphasized ARIMAX models' superiority over other modeling approaches [54][55][56][57][58].

Conclusions
This study investigated the influence of climate drivers on Western Australian rainfall variability and developed a forecast model to predict seasonal autumn rainfall using the ARIMAX technique. In this attempt, climate drivers were used as transfer functions, while all other previous attempts considered conventional time series models only. From a statistical perspective, it was evident that the predictability performance of the ARIMAX model is much higher than the MLR models. The ARIMAX models are capable of predicting rainfall 4 months in advance while the conventional MLR model can predict the rainfall only 1 month in advance. The developed ARIMAX models have shown a strong correlation (r) as well as minimum errors, smaller BIC values, and higher refined index of agreement values (d r ) in both calibration and validation periods, depicting the model's non-erroneous prediction capability. It was also observed that all these ARIMAX models were successful in predicting some of the extreme rainfall events and droughts, while their ability to predict seasonal autumn rainfall in advance of 4 months has strengthened their acceptability. From the stakeholder's perspective, such flexibility offered in the developed model has greater importance, as a timely prediction can help in strategic decision making and reducing associated risks and damage potentials. Overall, the DMI Oct -Nino3 Nov model for the south coast and DMI Jan -Nino3 Nov model for the north coast showed exceptional performance with good prediction accuracy and can be recommended for future rainfall prediction in WA. However, the developed models have some limitations as they were not able to predict all extreme cases. Investigating the existing nonlinear relationship between rainfall and climate drivers can provide a better understanding of the trend, and associated variabilities that the developed models failed to address. Possible analysis approaches can be considered for developing nonlinear models and hybrid models. Since rainfall is a complex mechanism, any linear or nonlinear model by itself, might not be able to predict or capture all the extreme cases. For such instances, a hybrid model could offer a better solution. The ARIMAX model residuals can be used to explain the nonlinear relationships, where the combined output of both ARIMAX and non-linear models can be used for improved forecasting.
Author Contributions: F.I. was involved in methodology development, data curation, formal analysis, investigation, software use, validation, writing-original draft preparation. M.A.I. was involved in conceptualization, resources allocation, supervision, project administration, and writing-review and editing. All authors have read and agreed to the published version of the manuscript.