## 4 Model Implementation and Validation
### Arima Modeling of SIC  
The SICs along the given NSR obtained in Section 4.1 are plotted in Figure 10 as a time series. The time series contains a decreasing trend and a seasonal period of 12 months. To stabilize variances further, we preprocess SICs in two steps: (1) adding 1 to the original SICs and (2) performing a logarithmic transformation, i.e., \(Y_{t}^{\prime}=\text{Log}\left(\widehat{Y}_{t}+1\right)\). The \(Y_{t}^{\prime}\) is then used to establish the ARIMA model. Therefore, the original SICs can be obtained from the back-transformation, \(\widehat{Y}_{t}=\text{e}^{Y_{t}^{\prime}}-1\). Following the flowchart (Figure 3), the tests of randomness and stationarity are required next. From Table 3, the Ljung-Box test suggests that the preprocessed SIC series in sub-region No. 2 contains relationships at lags of 1, 6, and 12 months. In other words, the series is non-random. Next, the Augmented Dickey-Fuller Test (ADF) tests the null hypothesis that a series is not stationary. Before carrying out the test, a seasonal difference, \(\widehat{\
abla}_{12}^{1}\,Y_{t}^{\prime}\) (Equation (11)) is taken for the series. The results in Table 4 show that the p-value is smaller than 0.05, indicating that the seasonally differenced series are stationary. The statistic value is smaller than the critical values of 5% and 10%, and very close to 1%, which also supports the stationarity of the series. \begin{table}
\begin{tabular}{c c c c c c c c} \hline
**Region** & **Linear** & **Power** & **Gaussian** & **Spherical** & **Exponential** & **Hole-Effect** & **Statistical** \\ \hline
1 & 0.139 & 0.110 & 0.101 & 0.117 & 0.093 & 0.076 & 0.043 \\
2 & 0.178 & 0.143 & 0.182 & 0.199 & 0.199 & 0.157 & 0.060 \\
3 & 0.071 & 0.076 & 0.100 & 0.101 & 0.081 & 0.102 & 0.036 \\
4 & 0.221 & 0.202 & 0.206 & 0.196 & 0.201 & 0.134 & 0.035 \\
5 & 0.157 & 0.154 & 0.159 & 0.162 & 0.154 & 0.171 & 0.022 \\
6 & 0.178 & 0.181 & 0.189 & 0.159 & 0.151 & 0.135 & 0.040 \\
7 & 0.124 & 0.118 & 0.199 & 0.152 & 0.130 & 0.143 & 0.029 \\
8 & 0.307 & 0.222 & 0.262 & 0.276 & 0.239 & 0.135 & 0.048 \\ \hline \end{tabular}
\end{table} 
Table 2: The errors of statistical interpolation and kriging with different variogram models. Figure 9: The comparison between kriging and statistical interpolation on SIC in October, November, and December. In the ARIMA \((p,d,q)\times(P,D,Q,S)\) model, orders \(d=0\), \(D=1\), and \(S=12\) have been confirmed. The orders \(p\), \(q\), \(P\), and \(Q\) can be determined from the autocorrelation function (ACF) and the partial autocorrelation function (PACF) in Figure 11, the blue shaded regions represent the range of zero plus/minus two standard errors. When values fall within the shaded region, the ACF and PACF are considered to be zero. The lags of the red dots provide information to determine the orders of the ARIMA model. Within the first seasonal period, i.e., lag \(h\in[0,12]\), the ACF cuts off at lags 3, 4, or 5, whereas the PACF cuts off at lag 2. In this regard, \(p\) is equal to 2, and \(q\) can be either 3, 4, or 5. Furthermore, the ACF has a peak at lag 12 and decays very rapidly at other multiples of the seasonal period, e.g., \(h=12k\), for \(k>1\). Therefore, \(Q\) is equal to 1. Nevertheless, the PACF decays slowly and oscillates at lags \(h=12k\), for \(k=1,2,\cdots\), and exhibits \"tailing off\" instead of \"cutting off.\" In this regard, \(P\) should be equal to 0. Now only the value of \(q\) is absent in the ARIMA \((2,0,q)\times(0,1,1,12)\) model. To this end, the cross-validation is designed to find the best ARIMA model that gives the lowest Akaike information criterion (AIC, [29]). The total set of 492 months is split into a training set and a test set in five different ways. The size of the test size is always 60 months, while the training set sizes vary from 60 to 432 months (Figure 12), each of fold uses the same test set size but monotonically increasing training set sizes, the test sets cover the period from month 432 to month 492.