The Application of Singular Spectrum Analysis in Modeling Earthquake Frequency in the Sumatra Subduction Zone

The island of Sumatra is located along the collision path of two tectonic plates, namely the India-Australia Plate and the Eurasian Plate. The boundary of the collision of these two plates is clearly visible in the form of a deep-sea trench in the west of Sumatra to the Andaman Islands, where the India-Australia Plate infiltrates under the island of Sumatra at a speed of 50-60 mm/year and a slope of about 12 degrees to the east [1]. The National Center for Earthquake Studies (PuSGeN) as one of the forums for researchers on seismicity in Indonesia has succeeded in identifying several 4 segments of large earthquake sources in the Sumatra megathrust zone, namely the Aceh-Andaman (AA), Nias-Simeuleu (NS), Mentawai-Pagai (MP), Mentawai-Siberut (MS), and Enggano (EO) segments. Visually, the limitations of the five segments of the source of the large earthquake can be seen in Figure 1 [2]. In the figure, information about the identity of the segment, the maximum potential earthquake is displayed, and the shear rate of tectonic plates in cm/year.

1.png

Figure 1. Source segment in the Sumatra megathrust zone

2.png

Figure 2. Research data preparation process

The data used in this study is mainshock earthquake. To obtain this data, a data preparation stage is required. The data preparation stage in this study uses the help of the ZMapp 7.1 package on MATLAB R2018a software and Microsoft Excel which follows the procedure as shown in Figure 2.

The first stage of research data preparation is to obtain data on earthquake events from January 1st, 1900 to December 31st, 2023. The data collection was carried out using the help of the ZMapp 7.1 package on the MATLAB R2018a software. In the process of collecting earthquake data using the ZMapp 7.1 package using a source from the USGS for the subduction zone of Sumatra Island with latitude (6° LU-11° LS) and longitude (94° BT-105° BT). The data to be taken is magnitude 2.5 until 10 Mw with a depth of ≤70 km. One of the reasons for choosing the depth of the earthquake is because an earthquake with this depth has the potential to cause great damage and can trigger a tsunami wave after an earthquake.

After knowing the value of M_c then the declustering process will be carried out. Declustering is the process of separating earthquake events between the main earthquake (mainshock) and the pioneer earthquake (foreshock), and aftershocks [3]. In this study, the declustering process uses the GardnerKnopoff method. GardnerKnopoff declustering is a method to separate pioneer or minor earthquakes that occur in groups or within the same time and location from major or large earthquakes in a short span of time [4]. The main goal of GardnerKnopoff declustering is to identify and remove shear earthquakes that occur due to major earthquake activity, thereby improving accuracy in statistical analysis and earthquake modeling. The GardnerKnopoff method is based on two assumptions, namely the assumption that major earthquakes can trigger pilot or minor earthquakes in a short span of time with the same location and the assumption that the time and location of pilot earthquakes can be considered as random variables [5]. The subsequent process is the uniformity of the magnitude type into one scale in the magnitude type M_w and sorting data with the provision data of value ≥M_c. 

2.1 SSA

SSA is one of the time series data analysis techniques used for forecasting. As we know, the SSA method is a non-parametric technique used to observe a time series data that is not linear, not stationary, has a variable nature at any time or temporarily, and more prone to success for short time-series [6, 7]. In general, SSA combines elements of classical time series analysis, multivariate statistics, multivariate geometry, dynamical systems, and signal processing. SSA can be very useful in solving various problems such as forecasting, imputation of missing values, decomposition of the original time series into a small number of interpretable components such as slowly changing trends, oscillatory components, and unstructured noise [8]. In SSA, the segmentation procedure produces smaller time series, which is favorable for SSA. If a time-series is non-stationary, transient oscillations can be evidenced locally in a given segment, whereas they could go unnoticed in the general computation of the whole time-series [9].

The main steps in forecasting using the SSA method are: Decomposition: The time series data is converted into a matrix (trajectories) using a window length parameter. This matrix is split using Singular Value Decomposition (SVD) into principal components that represent different patterns. Reconstruction: Relevant components are selected based on their contribution or significance (usually those with the largest singular values). These components are then recombined to reconstruct the time series into trends, seasonality, or noise as per the analysis. Forecasting: Based on the reconstructed components, a model is built to extend (forecast) the time series into the future while maintaining the identified patterns [6, 9]. The purpose of the first stage, decomposition, is to decompose the original series into a small number of independent and interpretable components such as slowly varying trends, oscillatory components, and unstructured noise. At this stage, capabilities for vibration analysis will also be applied. Then, in the second stage, reconstruction will be used to reconstruct the original signal using a selected number of components. For example, one can select only a certain number of components that contain important signal features, such as trend and oscillation components, and discard the (unstructured) noise components that cause discrepancies in the analysis results [10].

2.1.1 Embedding

Embedding is a stage where the initial time series data is converted into a trajectory X matrix, meaning converting the initial data in the form of one-dimensional data into multidimensional data. The trajectory X matrix has dimensions L×K. L shows the window length in the X matrix which is the number of rows, while K indicates the number of columns in the matrix. The selection of an appropriate value of the window length L is of great importance in SSA operation. This is because the efficiency of SSA in separating different signal components from the composite signal depends on the value of L [11].

In practice, the selection of window length (L) and clustering strategy are critical for optimal analysis results. The window length is usually chosen between 2 to N/2 (where N is the length of the data) and is tailored to the nature of the data: if the data has periodic patterns, the L should be close to a multiple of the period to capture seasonal patterns, while for data with dominant trends, a larger L is more suitable. The clustering strategy is performed after decomposition, where component clustering will be performed based on patterns, such as smooth components with large singular values for trends, oscillating components for seasonality, and components with random fluctuations and small singular values for noise.

There is no special method to determine the value of L precisely, so to determine the value of L is done by trial and error for the value to be in the range 2≤L≤2^-1 N and the time series data throughout the period N does not contain missing data and X={x_i}; i=1,2,…,N. Furthermore, in determining the value of K by K=N-L+1, where, N is the amount of data and L is the length of the windows [12]. The X trajectory matrix formed is a Hankel matrix where all its anti-diagonal elements have constant values. In the form of a matrix, it can be expressed as follows:

$\boldsymbol{X}=\left[x_1, x_2, \ldots, x_N\right]=\left[\begin{array}{ccccc}x_1 & x_2 & x_3 & \ldots & x_K \\ x_2 & x_3 & x_4 & \ldots & x_{K+1} \\ x_3 & x_4 & x_5 & \ldots & x_{K+2} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ x_L & x_{L+1} & x_{L+2} & \ldots & x_{L K}\end{array}\right]$        (1)

2.1.2 Singular Value Decomposition (SVD)

The SVD stage in its application has similarities with principal component analysis, which is to reduce components from the initial data and reduce dimensions [13]. The purpose of SVD is to obtain the separation of components in decomposition from time series data. Due to the high computational cost associated with SVD, its application to large datasets may not be suitable [14]. Note that, SVD begins with determining eigenvalue values $\left(\lambda_1, \lambda_2, \lambda_3, \ldots \lambda_L\right)$ from the matrix $\mathbf{S}=\mathbf{X} \mathbf{X}^{\mathrm{T}}$ where $\lambda_1 \geq \lambda_2 \geq \lambda_3 \geq \cdots \geq \lambda_L>0$, and eigenvector $\left(u_1, u_2, u_3, \ldots, u_L\right)$ from the $S$ matrix. The eigenvalue and eigenvector calculations of the $S$ matrix can be done using the following definitions:

Definition 1. If A is a matrix n×n, then the non-zero vector x is inside, Rⁿ. called an eigenvector (vector eigen) of A if Ax is a scalar multiple of x; that is

$\boldsymbol{A} x=\lambda x$

for a scalar λ. The scalar λ is called the eigenvalue of A and x is said to be the eigenvalue corresponding to λ [15].

To determine the eigenvalue of the $\boldsymbol{S}$ matrix measuring $L \times L$, so $\boldsymbol{A} x=\lambda x$ as $\boldsymbol{S} x=\lambda I x$ where in this case the matrix $\boldsymbol{S}$ is the matrix $\boldsymbol{A}$ referred to in Definition 1. The equation is equivalent to $(\lambda I-\boldsymbol{S}) x=0$. Thus, for the eigenvalue $\lambda$ equation $(\lambda I-\boldsymbol{S}) x=0$ has a non-zero if and only if:

$\operatorname{det}(\lambda I-\boldsymbol{S})=0$           (2)

Next, the eigenvalue $\lambda_i$ will be obtained by $i=1,2,3, \ldots, L$ and $U_i=\left(u_1, u_2, u_3, \ldots, u_L\right)$ is the corresponding eigenvector of eigenvalue.

Definition 2. Given $\boldsymbol{X}$ matrix $m \times n$ with rank $r$. The positive eigenvalue of $\sqrt{\boldsymbol{X} \boldsymbol{X}^{\boldsymbol{T}}}$ is called the singular value of $\boldsymbol{X}$ [16].

In other words, if $\sigma$ is the singular value of $X$ then $\sigma$ is the positive eigenvalue of $\sqrt{X X^T}$ or $\sigma^2$ is the eigenvalue of $X X^T$. From Definition 2, the relationship between the eigenvalue and the singular value is obtained. In other words, a singular value matrix is a matrix that is $\mathrm{L} \times \mathrm{K}$ in size with the main diagonal containing the positive root of the eigenvalue $\left(\sqrt{\lambda_i}\right)$ with the order $\lambda_1 \geq \lambda_2 \geq \lambda_3 \geq \ldots \geq \lambda_L>0$. It can be seen from the matrix below:

$\sqrt{\lambda_i}=\left[\begin{array}{ccc}\sqrt{\lambda_1} & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & \sqrt{\lambda_L}\end{array}\right]=\left[\begin{array}{ccc}\sigma_1 & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & \sigma_L\end{array}\right]=\sigma_i$   (3)

Furthermore, in determining the principal component matrix, can be searched using the relationship of singular values with the left singular vector and the right singular vector as follows:

Definition 3. Suppose the real matrix $\boldsymbol{X}$ is $m \times n$. Positive real numbers $\sigma$ called singular values of matrix $\boldsymbol{X}$ if there are zero vectors $\boldsymbol{U} \in \boldsymbol{R}^{\boldsymbol{m}}$ and $\boldsymbol{V} \in \boldsymbol{R}^{\boldsymbol{n}}$ so that $\boldsymbol{X} \boldsymbol{V}=\sigma \boldsymbol{U}$ and $\boldsymbol{X}^{\boldsymbol{T}} \boldsymbol{U}=\sigma \boldsymbol{V}$. The $U$ vector is called the left singular vector and $\boldsymbol{V}$ is called the right singular vector. Furthermore $(\sigma, U)$ is called the left singular pair of $\boldsymbol{X}$, and $(\sigma, V)$ is called the right singular pair of $\boldsymbol{X}$ [16].

From the above definition, the following equation is obtained:

$\boldsymbol{X} \boldsymbol{V}=\sigma_i \boldsymbol{U}$ and $\boldsymbol{X}^T \boldsymbol{U}=\sigma_i \boldsymbol{V}$

for vectors $\boldsymbol{V}_{\boldsymbol{i}}$ with $i=1,2,3, \ldots, L$ follows the following equation:

$\boldsymbol{V}_i=\frac{\boldsymbol{X}^T \boldsymbol{U}_i}{\sqrt{\lambda_i}}$    (4)

$\boldsymbol{V}_i^T$ is a transpose of the $\boldsymbol{V}_i$ matrix, so that the $\boldsymbol{V}_i^T$ matrix can be written in the form of:

$\boldsymbol{V}_i^T=\left[\begin{array}{cccc}v_{11} & v_{12} & \ldots & v_{1 L} \\ v_{21} & v_{22} & \ldots & v_{2 L} \\ \vdots & \vdots & \ddots & \vdots \\ v_{K 1} & v_{K 2} & \ldots & v_{K L}\end{array}\right]$       (5)

After obtaining the singular value $\left(\sigma_i\right)$, eigenvector $(\boldsymbol{U})$, and principal component $(\boldsymbol{V})$, the decomposition of the singular value of the trajectory matrix $\boldsymbol{X}$ is as follows:

Definition 4. For each $\boldsymbol{X} \in \boldsymbol{\Omega}^{m \times n}$ of rank $r$, there is an orthogonal matrix $\boldsymbol{U}_{m \times m}, V_{n \times n}$ and a diagonal matrix $\boldsymbol{D}_{r \times r}=$ $\operatorname{diag}\left(\sigma_1, \sigma_2, \ldots, \sigma_r\right)$ so that

$\boldsymbol{X}=\boldsymbol{U}\left(\begin{array}{ll}\boldsymbol{D} & \mathbf{0} \\ \mathbf{0} & \mathbf{0}\end{array}\right)_{m \times n} \boldsymbol{V}^T$            (6)

with $\sigma_1 \geq \sigma_2 \geq \cdots \geq \sigma_r>0$ [17].

The $\boldsymbol{V}^{\boldsymbol{T}}$ matrix is a transpose of the $\boldsymbol{V}$ and $\boldsymbol{U}$ matrices of $L \times L$ and $\boldsymbol{D}_{r \times r}$ is a matrix with $L \times K$ dimensions whose main diagonal element is the singular value of the $S$ matrix, so that the value of $X_i$ can be found by the following steps:

$\begin{gathered}\boldsymbol{X}_{\boldsymbol{i}}=\left[\begin{array}{cccc}u_{11} & u_{12} & \ldots & u_{1 L} \\ u_{21} & u_{22} & \ldots & u_{2 L} \\ \vdots & \vdots & \ddots & \vdots \\ u_{L 1} & u_{L 2} & \ldots & u_{L L}\end{array}\right]\left[\begin{array}{cccc}\sigma_1 & 0 & 0 & 0 \\ 0 & \sigma_2 & 0 & 0 \\ 0 & 0 & \ddots & 0 \\ 0 & 0 & 0 & \sigma_L\end{array}\right]\left[\begin{array}{cccc}v_{11} & v_{12} & \ldots & v_{1 k} \\ v_{21} & v_{22} & \ldots & v_{2 K} \\ \vdots & \vdots & \ddots & \vdots \\ v_{K 1} & v_{K 2} & \ldots & v_{K K}\end{array}\right]\\=\boldsymbol{U}_1 \sigma_1 \boldsymbol{V}_1^T+\boldsymbol{U}_2 \sigma_2 \boldsymbol{V}_2^T+\cdots+ \boldsymbol{U}_d \sigma_d \boldsymbol{V}_d^T=\sum_{i=1}^d \boldsymbol{U}_i \sigma_i \boldsymbol{V}_i^T\end{gathered}$         (7)

with $i=1,2,3, \ldots, d$ and $d=\max \{i\}$ and $\sigma_i>0$.

The three components in the $\boldsymbol{X}_{\boldsymbol{i}}$ matrix namely singular value $\left(\sigma_i\right)$, eigenvector $\left(\boldsymbol{U}_i\right)$, and principal component $\left(\boldsymbol{V}_i\right)$ are called the i'th eigentriple of Singular Value Decomposition. So the SVD for the $\boldsymbol{X}_{\boldsymbol{i}}$ matrix can be written as follows [6]:

$\boldsymbol{X}_i=\boldsymbol{X}_1+\boldsymbol{X}_2+\boldsymbol{X}_3+\cdots+\boldsymbol{X}_d$.           (8)

2.2 Reconstruction

Reconstruction is a stage where the data is reconstructed into new time series data based on the values obtained in the previous stage, namely the decomposition stage [6]. The newly obtained time series data will go through two processes, namely grouping and diagonal averaging.

2.2.1 Grouping

Grouping is the grouping of a matrix of $\boldsymbol{X}_i$ into several groups. The grouping of the $\boldsymbol{X}_{\boldsymbol{i}}$ matrix can be done by subjectively looking at eigentriple value patterns or can perform spectral analysis to see patterns and periods in eigentriples. The grouping process is carried out by grouping the set of indices $\{1,2,3, \ldots, d\}$ into $m$ subsets denoted by $I=$ $I_1, I_2, I_3, \ldots, I_m$ with $d=m$. Then $\boldsymbol{X}_I$ adjusted to the group $I=\left\{I_1, I_2, I_3, \ldots, I_m\right\}$. Then $\boldsymbol{X}_i=X_1, X_2, X_3, \ldots, \boldsymbol{X}_m$ can be expanded to:

$\boldsymbol{X}_I=\boldsymbol{X}_{I 1}+\boldsymbol{X}_{I 2}+\boldsymbol{X}_{I 3}+\cdots+\boldsymbol{X}_{I m}$.      (9)

The stage to select the set $\mathrm{I}=\left\{\mathrm{I}_1, \mathrm{I}_2, \mathrm{I}_3, \ldots, \mathrm{I}_{\mathrm{m}}\right\}$ is called eigentriple grouping which is carried out by trial and error. The trial-and-error method can be done subjectively, that is, trying all possibilities by looking at the correlation plot or relationship of each eigentriple.

The concept of forecasting in the SSA method is a separation that characterizes how well different components can separate from each other. The decomposition of SSA from F_N series will only be successful if the additive components of the series can be separated from each other. The natural measurement of similarity between two series X⁽¹⁾ and X⁽²⁾ uses the W–correlation equation which is formulated as follows [6]:

$\rho^{(w)}\left(\mathbb{X}^{(1)}, \mathbb{X}^{(2)}\right)=\frac{\left(\mathbb{X}^{(1)}, \mathbb{X}^{(2)}\right)_w}{\left\|\mathbb{X}^{(1)}\right\|_w\left\|\mathbb{X}^{(2)}\right\|_w}$          (10)

with $\left(\mathbb{X}^{(1)}, \mathbb{X}^{(2)}\right)_w=\sum_{i=0}^N w_i x_i^{(1)} x_i^{(2)}$ and $w_i$ searched by using:

$w_i=\left\{\begin{array}{c}i \quad \text { for } 0 \leq i<L^* \\ L^* \quad \text { for } L^* \leq i \leq K^* \\ N-i+1 \text { for } K^*<i \leq N\end{array}\right.$

with $L^*=\min (L, K)$ and $K^*=\max (L, K)$.

3.png

Figure 3. Example of W-correlation plot

If the absolute value of W-correlation is small, then the series almost meets the W-orthogonal. On the contrary, the sequence is far from the W-orthogonal and is very separate. So, if both reconstructed components have a W-correlation value of zero, it means that the two components can be separated. The W-correlation plot is used as a visualization of the magnitude of the correlation between eigentriples [6]. The older the color, the higher the correlation. An example of a W-correlation plot can be seen in Figure 3.

2.2.2 Diagonal averaging

After grouping, the next stage will be transformed from the grouping results, namely diagonal averaging. This step is done by converting each matrix $\boldsymbol{X}_I$ into a new time series with a length of $N$. Suppose the matrix $\boldsymbol{F}$ is an arbitrary matrix of $\mathrm{L} \times \mathrm{K}$ size with $f_{i j}$ elements with $1 \leq \mathrm{i} \leq \mathrm{L}$ and $1 \leq \mathrm{j} \leq \mathrm{K}, L^*=$ $\min (L, K), K^*=\max (L, K)$, and $N=L+K-1$.

$\boldsymbol{F}=\left[\begin{array}{ccccc}f_{11} & f_{12} & f_{13} & \ldots & f_K \\ f_{21} & f_{22} & f_{23} & \ldots & f_{K+1} \\ f_{31} & f_{32} & f_{33} & \ldots & f_{K+2} \\ \vdots & \vdots & \vdots & \ddots & f_{K+3} \\ f_L & f_{L+1} & f_{L+2} & \ldots & f_N\end{array}\right]$      (11)

The matrix $\boldsymbol{F}$ is transformed into a series of $g_1, g_2, g_3, \ldots, g_N$ through diagonal averaging with the following formula [6]:

$g_k=\left\{\begin{array}{l}\frac{1}{k+1} \sum_{m=1}^{k+1} y_{m, k-m+1}^*, \text { for } 0 \leq k<L^*-1 \\ \frac{1}{L^*} \sum_{m=1}^{L^*} y_{m, k-m+2}^*, \text { for } L^*-1 \leq k<K^* \\ \frac{1}{N-k} \sum_{m=k-K^*+2}^{N-K^*+1} y_{m, k-m+2}^*, \text { for } K^* \leq k<N\end{array}\right.$             (12)

Based on Eqs. (11) and (12), for example, in the matrix $\boldsymbol{F}$ is chosen for $k=1$ to give the result $g_1=f_{11}$, for $k=2$ gives $g_2=\frac{f_{12}+f_{21}}{2}$, for $k=3$ gives $g_3=\frac{f_{13}+f_{22}+f_{31}}{3}$ and so on. Meanwhile, based on Eq. (12), if applied to the $X_{I m}$ resultant matrix, a series $\tilde{g}^{(k)}=\left(\tilde{g}_1^{(k)}, \tilde{g}_2^{(k)}, \ldots, \tilde{g}_N^{(k)}\right)$ so that the original series is decomposed into:

$g_k=\left\{\begin{array}{l}\frac{1}{k+1} \sum_{m=1}^{k+1} y_{m, k-m+1}^*, \text { for } 0 \leq k<L^*-1 \\ \frac{1}{L^*} \sum_{m=1}^{L^*} y_{m, k-m+2}^*, \text { for } L^*-1 \leq k<K^* \\ \frac{1}{N-k} \sum_{m=k-K^*+2}^{N-K^*+1} y_{m, k-m+2}^*, \text { for } K^* \leq k<N\end{array}\right.$        (13)

2.3 Linear Recurrent Formula (LRF)

Forecasting in the SSA method uses LRF or known as R-Forecasting. SSA-LRF is a non-parametric method that is widely used for time series forecasting and does not accommodate random behavior (noise) in the analysis set because the goal is to separate the signal from the noise component [18]. The time series used is the value of the reconstruction series obtained from the averaging diagonal results. R-Forecasting forecasting is related to the estimation of the LRF, namely $a_1, a_2, a_3, \ldots, a_d$ which is an eigenvector obtained from the SVD stage [19]. The prediction model is as follows:

$g_i=\left\{\begin{array}{c}\tilde{y}_l, \quad i=0,1,2, \ldots, N \\ \sum_{j=1}^{L-1} a_j g_{i-j}, i=N+1, N+2, \ldots, N+M\end{array}\right.$           (14)

If for example $\boldsymbol{U}=\left(\boldsymbol{u}_1, \boldsymbol{u}_2, \boldsymbol{u}_3, \ldots, \boldsymbol{u}_{\boldsymbol{L}-1}, \boldsymbol{u}_L\right)^{\boldsymbol{T}}, \boldsymbol{U}^{\boldsymbol{\nabla}}=$ $\left(\boldsymbol{u}_1, \boldsymbol{u}_2, \boldsymbol{u}_3, \ldots, \boldsymbol{u}_{L-1}\right)^{\boldsymbol{T}}$, and $\boldsymbol{\pi}_{\boldsymbol{i}}$ is the last component of the eigenvector $\boldsymbol{U}$ or can be written $\pi_i=u_L$ then the LRF coefficient can be calculated [19]:

$\mathfrak{R}=U=\left(a_{L-1}, \ldots, a_2, a_1\right)^T=\frac{1}{1-v^2} \sum_{i=1}^r \pi_i \boldsymbol{U}_i^{\nabla}$    (15)

with $v^2=\sum_{i=1}^r \pi_i^2$.

Furthermore, based on Eq. (14), a series of forecasting results will be formed which are shown as follows: 

$G_{N+M}=\left(g_1, g_2, g_3, \ldots, g_{N+M}\right)$          (16)

2.4 Forecasting reliability

Point forecast accuracy measures vary based on (i) the type of error (e.g., absolute or squared) considered for their calculation, (ii) the type of scaling (e.g., mean of the series, first differences of the series, or relative performance against a benchmark method) potentially used for making them scale-independent, and (iii) the operator (e.g., mean, median, or geometric mean) utilized for aggregating the errors across different forecast horizons and series [20].

The degree of accuracy of forecasting methods can be explained by comparing the values obtained from forecasting with actual values. The accuracy of a forecasting method is determined by the smallest value of each method, the smaller the error value, the more accurate a method is in forecasting. To measure forecasting accuracy, two popular measures can be used, namely RMSE and MAAPE.

2.4.1 Root Mean Square Error (RMSE)

RMSE is the sum of squared errors or the difference between the actual value and the predicted value that has been determined. The equation for RMSE is as follows [21]:

$R M S E=\sqrt{\sum \frac{\left(g^{\prime}-g\right)^2}{n}}$        (17)

2.4.2 Mean Arctangent Absolute Percentage Error (MAAPE)

MAAPE is a new measure of accuracy developed by observing Mean Absolute Percentage Error (MAPE) from different perspectives. MAAPE measure is the slope as an angle, while MAPE is a slope as a ratio, considering that a triangle with adjacent and opposite sides is equal to the real value and the difference between the actual value and the prediction of each [22].

$M A A P E=\frac{1}{n} \sum_{i=1}^n \arctan \left|\frac{g_i-\hat{g}_i}{g_i}\right|$          (18)

In addition to using these two accuracy measures, in this study, tests will be carried out on residual white noise in the form of residual white noise test and normally distributed residual test which can be seen as follows [23]:

2.4.3 Residual white noise test

The residual white noise test is written as follows:

    a. Hypothesis

$\mathrm{H}_0: \rho_1=\rho_2=\cdots=\rho_k=0$ (no residual autocorrelation)

$\mathrm{H}_1$: There is at least one $\rho_{{j}} \neq 0$ (there is a residual autocorrelation)

    b. Test statistics, namely Ljung-Box or Box-Pierce Modified test statistics:

$Q^*=n(n+2) \sum_{k=1}^K \frac{\hat{\rho}_k^2}{(n-k)}$       (19)

    c. Decision-making criteria:

If the p-value is > alpha (the level of significance set), then H_0 is accepted and it can be concluded that there is no significant residual autocorrelation whereas if the p-value is ≤ alpha, then H_0 rejected and it can be concluded that there is a significant residual autocorrelation. If there is no significant residual autocorrelation, then the model is considered suitable for the data and is good at forecasting.

2.4.4 Normally distributed residual test

In data analysis, given that various statistical methods assume that the distribution of the population data is normal distribution, it is essential to check and test whether or not the data satisfies the normality requirement [24]. The Anderson-Darling test is a modification of the Kolmogorov-Smirnov test. The critical values in the Kolmogorv-Smirnov test do not depend on the specific distribution being tested, whereas the Anderson-Darling test makes use of a specific distribution. In the Anderson-Darling test, the hypothesis used is:

$\mathrm{H}_0$: Residual is normally distributed

$\mathrm{H}_1$: Residual is not normally distributed

The decision is that if the p-value is ≤α, then H_0 rejected and it can be concluded that the residual distribution of the model is not normal while if the p-value>α (the level of significance), then $\mathrm{H}_0$ accepted, that is, the residual distribution of the model is normal [25]. In this case, it can be considered that the model is suitable for the data modeled.

This study demonstrates that SSA is an effective method for modeling and predicting earthquake frequency in the Sumatra subduction zone, a region with significant tectonic activity. The findings highlight SSA’s ability to analyze complex, non-linear, and non-stationary time series data, offering a novel approach to earthquake forecasting. A key contribution of this research lies in its comprehensive evaluation of four different insample-outsample data splits, which revealed that a 90% insample and 10% outsample composition provided the best performance for insample predictions. Additionally, the study underscores the practical relevance of SSA in addressing disaster preparedness in one of the world's most hazardous seismic zones. However, the reliability of SSA depends heavily on the availability of high-quality historical data, and its ability to forecast extreme events, such as megathrust earthquakes, remains uncertain due to the complexity of tectonic dynamics.

Future research could enhance SSA’s predictive capabilities by integrating it with techniques such as Machine Learning or statistical models like ARIMA or GARCH, applying it to other seismic regions to test its generalizability, and incorporating additional geophysical data such as tectonic stress or surface deformation. These advancements could broaden the applicability of SSA in disaster mitigation and improve preparedness in earthquake-prone regions.