Systematic Literature Review on An Integrated Generalized Space Time Autoregressive Integrated Moving Average (GSTARIMA) Model with Heteroscedastic Error and Kriging Method for Forecasting Climate

1. Introduction

Climate is a statistical description of the average variability of the relevant quantities over months to years, referred to as average weather [1]. In addition, it includes several interrelated elements, such as temperature, rainfall, humidity, atmospheric conditions, and wind patterns [2]. Climate change is a pressing global issue of paramount importance that demands comprehensive research. The Intergovernmental Panel on Climate Change (IPCC) is an organization founded by scientists worldwide to research the concept. The sixth assessment report of the IPCC explains that climate change affects ecosystem conditions, human activities, the global water cycle, infrastructure, health, and others [3].

The handling of climate change is in the world's spotlight, which is included in the pillars of SDGs [4,5]. Meanwhile, climate management is the 13^th goal in the SDGs, with the mission statement "Take urgent action to combat climate change and its impacts by regulating emissions and promoting developments in renewable energy” [6]. One of the climate changes very important to study is rainfall patterns [7,8]. Based on NASA's Global Precipitation Measurement, the Indonesian region has high rainfall around 1000-4000 mm per year. This is because Indonesia is located on the equator, and it is vulnerable to natural disasters such as floods and landslides. Based on data from the West Java Regional Disaster Management Agency (known as BPBD), which can be accessed on the website https://opendata.jabarprov.go.id/id/ (accessed on 01/04/2023), there were 1954 and 5,662 flood and landslide events in 2012-2021. The number of flood and landslide events for each district and city in West Java is presented in Figure 1.a and 1.b.

Natural disasters due to extreme rainfall have a significant impact and cause damage to community settlements. Based on data from the West Java BPBD report, settlement damage totaled 943,160 units from 2012-2021. Damage is categorized into destroyed, heavily damaged, slightly damaged, moderately damaged, threatened, and submerged/buried, as presented in Table 1.

Climate change is very detrimental regarding materials, infrastructure, and people's lives. Therefore, it is essential to forecast future climate conditions to take preventive, mitigation, and adaptation actions.

Climate forecasting can also be conducted using the Spatio-Temporal Model. The Spatio-Temporal Model combines location and time to model a phenomenon to understand the relationship between changes in spatial and time. The Generalized Space-Time Autoregressive (GSTAR) Model assumes heterogeneous characteristics between locations and stationary data. Furthermore, it has different autoregressive and Space-Time parameters for each location. The Model was studied using stationary data by Borovkova et al.[9]. Giacinto developed the Model into the Generalized Space-Time ARMA (GSTARMA) Model [10]. In addition, the GSTARMA Model adds the effect of an error element with a Moving Average (MA) and is applied to stationary data. The Model Estimation uses the Maximum Likelihood Estimation (MLE) method and is applied to forecast the unemployment rate in each region of Italy based on historical data. The Model of non-stationary data is called the GSTARIMA Model, developed by Min et al., 2010 [11].

The GSTAR Model has a non-constant error variance (heteroscedastic error) for climate data. The GSTAR-ARCH Model is an extension of the Spatio-Temporal Model, which considers the heteroscedasticity of variance that depends on previous information in an autoregressive manner by Nainggolan et al [12]. The Model is applied to stationary data as an extension of the GSTAR-ARCH model for non-stationary data by Bonar et al. [13]. The GSTARI-ARCH model also overcomes the non-constant error variance. Bonar et al. used the concept to model and forecast the Consumer Price Index (CPI) in North Sumatra. In forecasting Spatio Temporal Model, the respon variable is influenced by exogenous variables. For example, in climate data, rainfall is affected by humidity and temperature. The Spatio-Temporal model with the addition of exogenous variables is known as the GSTARI-X Model by Elfiyan et al. [14]. Ditago et al. (2016) used the GSTARX-GLS model, with the exogenous variable being calendar variation[15]. Monika et al. (2022) developed the GSTARI-X-ARCH model to forecast rainfall with exogenous variables in humidity [16].

In previous research, the Kriging Method was used to predict phenomena at unobserved locations. Kriging method is used for interpolation and forecasting Temperature in Mosul and Baghdad City [17]. Kriging method, land-use regression (LUR), and LightGBM (light gradient boosting machine) methods were combined to predict PM2.5 concentrations [18]. In Spatio-Temporal modeling, the GSTAR Model is integrated with the Kriging method to forecast rainfall at unobserved locations in West Java [19].

This study aims to summarize previous research on Spatio Temporal forecasting models with heteroskedastic errors and the Kriging method applied to climate data. This research attempts to cover several areas, such as Spatio Temporal models for stationary and non-stationary data, methods for parameter estimation in the models, forecasting at unsampled locations, and the potential to integrate Spatio Temporal models with Heteroskedastic errors and Kriging for climate forecasting. Ultimately, this review contributes to a broader understanding of integrated Spatio Temporal Models with Heteroskedastic errors and the Kriging Method for climate and highlights avenues for further research and innovation in this critical area. To facilitate the analysis process, we formulate the following research questions (RQs):

RQ1:

How to integrate GSTARIMA model with heteroskedastic errors using Kriging method?

RQ2:

How to forecast climate phenomena using the integration of GSTARIMA and Kriging models through a data analysis life cycle approach?

RQs were examined and explored by reviewing previous results carried out by searching literature on databases. The results were filtered and selected using the Preferred Reporting Items for Systematic Review and Meta-Analyses (PRISMA) method. Furthermore, relevant articles are presented in a state of the art to obtain research gaps. A bibliometric method was also used to show the linkage of keywords for each article. The review stage was performed to analyze search results and discuss new research. Potential new research was provided to be studied and developed on the GSTARIMA Model and its application.

2. Materials and Methods

2.1. Theoritical Background

2.1.1. The Generalized Space Time Autoregressive Integrated Moving Average (GSTARIMA)Model

In 1980, Pfeifer and Deutch introduced the Space-Time Autoregressive (STAR) Model, assuming each location has the same characteristics [20,21]. In 2002, Ruchjana developed the STAR model into the Generalized Space-Time Autoregressive (GSTAR) model. This is because the assumptions in the STAR model do not match the reality in the field, where there is a diversity of characteristics at each location. The GSTAR model introduced by Ruchjana assumes that the characteristics of each location are heterogeneous. The GSTAR$(p,{\lambda }_k)$ model has a time order of $p$ and a spatial order of ${\lambda }_k$ expressed in matrix form through equation (1)[9]

$$\mathbf{z}\left(t\right)=\sum _{k=1}^p\sum _{l=0}^{{\lambda }_k}\left[{\mathsf{\Phi }}_{kl}{\mathbf{W}}^{\left(l\right)}\mathbf{z}\left(t-k\right)\right]+\mathbf{e}\left(t\right)$$

(1)

where,

$Z\left(t\right)$  

: the value of the observation at time $t$,

$Z\left(t-1\right)$ 

: the value of the observation at time $t-1$,

  $\varphi$ 

 : a parameter that indicates the influence of the value of $Z(t-1)$ on the value of $Z\left(t\right),$

$e\left(t\right)$   

 : the value of error.

The GSTARMA model expands the GSTAR model by adding MA error elements. The GSTARMA model is applied to stationary data [10]. The GSTARMA model developed on nonstationary data is called the GSTARIMA model. Min et al. (2010) first introduced the GSTARIMA model with application to urban traffic network modeling and short-term traffic flow forecasting. The GSTARIMA model $(p_{{\lambda }_k},d,q_{v_k})$with $d$ being the differencing order is expressed in Equation (2) [11]

$$\mathbf{y}\left(t\right)=\sum _{k=1}^p\sum _{l=0}^{{\lambda }_k}\left[{\mathsf{\Phi }}_{kl}{\mathbf{W}}^{\left(l\right)}\mathbf{y}\left(t-k\right)\right]-\sum _{k=1}^q\sum _{l=0}^{v_k}\left[{\mathsf{\Theta }}_{kl}{\mathbf{W}}^{\left(l\right)}\mathbf{e}\left(t-k\right)\right]+\mathbf{e}\left(t\right),$$

(2)

where,

$$\mathbf{y}\left(t\right)=\mathbf{z}\left(t\right)-\mathbf{z}\left(t-1\right),\mathbf{y}\left(t-1\right)=\mathbf{z}\left(t-1\right)-\mathbf{z}\left(t-2\right),\dots ,\mathbf{y}\left(t-k\right)=\mathbf{z}\left(t-k\right)-\mathbf{z}\left(t-k-1\right)$$

(3)

$\mathbf{z}\left(t\right)$

  : a vector of variables of size $(N\times 1)$ at time $t$,

$\mathbf{z}\left(t-k\right)$

: vector of variables of size $(N\times 1)$ at time $(t-k),$

${\lambda }_k$

    : spatial order in the $k$^th autoregressive,

$v_k$

    : spatial order of the $k$^th moving average,

${\mathsf{\Phi }}_{kl}$     

: autoregressive and space time parameters at time order $k$ and spatial order $l$ of size $(N\times N)$ in the form of diagonal matrix $\left({\Phi }_{kl}^{\left(1\right)},{\Phi }_{kl}^{\left(2\right)},{\Phi }_{kl}^{\left(3\right)},\dots ,{\Phi }_{kl}^{\left(N\right)}\right)$,

${\mathsf{\Theta }}_{kl}$  

: MA parameters at time order $k$ and spatial order $l$ of size $(N\times N)$ in the form of diagonal matrix $\left({\mathsf{\Theta }}_{kl}^{\left(1\right)},{\mathsf{\Theta }}_{kl}^{\left(2\right)},{\mathsf{\Theta }}_{kl}^{\left(3\right)},\dots ,{\mathsf{\Theta }}_{kl}^{\left(N\right)}\right),$

${\mathbf{W}}^{\left(l\right)}$    

: weight matrix of size $\left(N\times N\right)$at spatial order $l,l=\mathrm{0,1},2,\dots ,{\lambda }_k$containing $w_{ii}=0$ and $\sum _{i\ne j}{w}_{ij}=1$,

$\mathbf{e}\left(t\right)$

 : error vector of size $(N\times 1)$ at time $t$, assuming ${\mathbf{e}\left(t\right)}_{~}^{iid}N\left(0,{\sigma }^2\mathbf{I}\right).$

2.1.2. Autoregressive Conditional Heteroscedasticity (ARCH) Model

Although the GSTARIMA model assumes constant error variance, applying climate data often shows non-constant error variance. The GSTARIMA model is integrated with the Autoregressive Conditional Heteroscedasticity (ARCH) Model to overcome this. This time series model detects variance heteroscedasticity using historical data [22]. Describing the ARCH($p$) model, researchers use the following expression [22]

$$h_t={\sigma }_t^2={\alpha }_0+{\sum }_{i=1}^p{\alpha }_i{\mathbf{e}}_{t-i}^2;i=\mathrm{1,2},3,\dots ,p,$$

(4)

In Equation (4), the variables represented include:

$h_t$

        : the conditional variance at time $t,$

${\alpha }_0$

       : the intercept or constant error,

${\alpha }_1,{\alpha }_2,\dots ,{\alpha }_p$

: ARCH model parameters,

${\alpha }_0>0$ dan ${\alpha }_i\ge 0$.

 

2.1.3. Kriging Method

The Kriging method is a geostatistical interpolation technique used to predict variable values at unobserved locations based on variable values observed at other locations. This method assumes that variable values have a spatial structure related to the distance and direction between observation locations. In the calculation of the Kriging Method, a Semivariogram is required. An experimental semivariogram is calculated based on measurement data collected from the field or observations at a particular location. The formula for calculating the experimental semivariogram is as follows [23]:

$$\hat{\psi }(h)=\frac{1}{2N\left(h\right)}\sum _{i=1}^{N\left(h\right)}\left[Z\right(x_i+h)-Z(x_i){]}^2$$

(5)

where,

$\hat{\psi }\left(h\right)$

   : semivariogram value at distance $h$,

$Z\left(x_i\right)$

     : observation value at location $x_{i,}$

$Z(x_i+h)$

    : observation value at location $x_i+h,$

$N\left(h\right)$

    : many pairs of data that have distance $h$,

$h$

      : distance between 2 locations.

Theoretical semivariograms can be divided into Spherical, Gaussian, and Exponential Models. The Spherical Model is a model that assumes that spatial dependence has a certain maximum distance or radius. This Model is used if the spatial dependence decreases with distance and reaches a threshold value at a specific radius, after which the semivariogram value becomes constant. The semivariogram function of the Spherical Model can be expressed as [23]:

$$\psi \left(h\right)=\left\{\begin{array}{c}c\left[\left(\frac{3h}{2a}\right)-{\left(\frac{h}{2a}\right)}^3\right],h\le a\\ c\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em},h>a\hfill \end{array}\right.$$

(6)

The Exponential Model is a model that assumes that spatial dependence decreases exponentially with distance between locations. The semivariogram function of the Exponential Model can be expressed as [23]:

$$\psi \left(h\right)=\left\{\begin{array}{c}c\left[1-exp\left(\frac{-h}{a}\right)\right],h\le a\\ c\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em},h>a\hfill \end{array}\right.$$

(7)

The Gaussian Model is a model that assumes that spatial dependence has a symmetric pattern and decreases exponentially with distance between locations. The semivariogram function of the Gaussian Model can be expressed as [23]:

$$\psi \left(h\right)=\left\{\begin{array}{c}c\left[1-exp{\left(\frac{-h}{a}\right)}^2\right],h\le a\\ c\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em},h>a\hfill \end{array}\right.$$

(8)

where,

$h$

: distance between sample locations,

$c$

: sill value,

$a$

: range.

The semivariogram also provides the weights used in interpolation. The Kriging method aims to determine the value of the Kriging weight ${\theta }_i$, which minimizes the estimator's variance so that a BLUE (Best Linear Unbiased Estimator) estimator is obtained. The Kriging estimator $\hat{Z}\left(x_0\right)$can be written as follows [23]:

$$\widehat{Z}\left(x_0\right)-\xi \left(x_0\right)=\sum _{i=1}^n{\theta }_i\left[Z\left(x_i\right)-\xi \left(x_i\right)\right],$$

(8)

where,

$\hat{Z}\left(x_0\right)$

   : Kriging estimator at unobserved location $x$,

$x_i$

   : the i^th data location adjacent to the unsampled location $x$,

$\xi \left(x_0\right)$

 : expectation value of $Z\left(x_0\right)$,

$\xi \left(x_i\right)$

 : expectation value of $Z\left(x_i\right)$

$n$

    : many sample data used for estimation,

${\theta }_i$

  : weight value at location $i$.

2.1.4. Data Analytics Life Cycle

Climate data has the Big Data criteria of volume, variety, and velocity. Big Data could be more efficient when analyzed using traditional methods. The Data Analytics Life Cycle methodology is specifically designed to handle Big Data problems and data science projects. The Data Analytics Life Cycle consists of six phases, including [24]:

• Discovery -> At this stage, researchers must study, search and investigate facts, identify problems, and develop context and understanding of the data sources needed to support research.
• Data Preparation -> Next, data is cleaned to identify missing values or noisy data. The results of data cleaning are transformed by aggregating daily data into monthly or according to the needs of the analysis. In this case, pre-processing data is obtained and ready for processing and analysis.
• Model Planning -> At this stage, the model planning that will be used for analysis is carried out.
• Model Building -> Researchers divide the results of data preparation into in-sample data (training) and out-sample data (testing) to do forecasting.
• Communicate Results -> Researchers present forecasting results using visualizations in the form of time series plots, choropleth maps, diagrams, and others.
• Operationalize -> The final stage is operationalized, and researchers provide final reports, recommendations, scripts, and technical documents. In addition, researchers can also apply the Model to the appropriate environment.

2.2. Collected Article

The PRISMA method is a widely used guide and methodological framework for conducting and presenting systematic reviews and meta-analyses [25]. The method provides the results of a systematic review, including completeness and clarity in reporting. The PRISMA method is supported by flowcharts in selecting articles [26,27].

The first stage in the PRISMA method is a literature search. Meanwhile, literature search through keywords was carried out in this research in four databases, namely Google Scholar, Dimensions, Science Direct, and Scopus. The keywords entered in the database consist of four codes connected with "OR" and "AND." The criteria selected in the collection of articles include:

• The publication type selected is article research and conference paper.
• Written in English
• The range of article publications is 2000-2023.
• The title, abstract, or keywords contain the words presented in Table 2.

The keywords provided in Table 2 are input into the database, followed by pressing the enter key to initiate a search. After displaying the search results, criteria 1-3, which pertain to the publication type, language selection, and publication year range, are configured to filter articles under the specified parameters. Subsequently, eligible articles are downloaded in .bib, .csv, and .ris formats. The number of article findings in each database is recorded for utilization as reference material in the subsequent stage.

The second stage involves the selection of articles, which is carried out through a manual process to ensure relevance. Specifically, the criteria for selecting relevant articles are those that explore the GSTAR model and its application. The articles included at this stage comprise both the ones obtained from the initial database search and the ones found manually through citation searching. The stages in article selection are explained as follows [28,29,30]:

(a)

Duplicate selection aims to remove duplicate articles found. Duplication can be found in databases or literature sources with almost the same or similar structure. Duplication selection stage can be conducted with special software such as Jabref and Mendeley reference managers to compare titles, abstracts, and content.

(b)

The relevance of the title and abstract is selected by assessing and ensuring that it matches the topic criteria. Titles and abstracts of selected articles are read in their entirety and irrelevant ones are excluded at this stage.

(c)

The full selection aims to determine whether the discussion and content in the article are relevant to the topic. All articles are accessed and read manually to ensure their appropriateness. Articles that fail to meet the established criteria or do not pertain to the subject matter under investigation are hereby excluded from the subsequent phases of the process.

The final stage in the PRISMA method is the articles review, explaining, and answering the RQs presented in Section 1.

3. Results

3.1. Results of Literature Search and Dataset Analysis

The results of the literature search are presented in Table 3, where code A produces 213,557 articles, code B produces 1,121,262, code C produces 7,525,693, and code D is searched by combining code A, B, and C the "AND" connector to produce 286 articles.

The manual selection stage of the article is carried out as follows:

(a)

At the initial stage, duplicate selection is conducted to identify 161 articles as duplicates and removed from the study.

(b)

The selection stage is based on the relevance of the title and abstract, where 35 articles are selected as relevant and considered for further research.

(c)

In the full paper accessibility selection stage, a total of 60 articles can be accessed and downloaded for further selection.

(d)

In the full paper relevance selection stage, the entire contents of the 18 articles are read and analyzed to determine their relevance. Relevant papers were also added from another method with citation search, resulting in 32 relevant articles. So that a total of 48 review articles are obtained that are relevant to the topic discussed.

These stages are presented visually in the PRISMA diagram in Figure 2 with three stages, namely identification, screening, and inclusion. Identification includes the duplication selection stage in stage (a). Screening consists of stages (b) and (c) for selecting title-abstract and full paper. Finally, inclusion explains the number of research articles relevant to the topic.

3.2 Bibliometric Analysis

The next stage describes the selected articles in bibliometric mapping used as a visualization method to analyze the pattern of relationships between scientific articles [31,32,33]. This paper uses bibliometric maps to visualize scientific networks involving keywords in 48 articles. The visualization results are in the form of circles and clusters distinguished by different colors. The circles on the bibliometric map represent the number of related publications by keyword. A circle with a large size indicates several keywords with similar relationships between scientific articles. Clusters in a bibliometric map show connected circles and represent scientific articles with similarities in context, such as topics [34]. Furthermore, bibliometric mapping keyword analysis is obtained using VOSviewer to understand the structure, patterns, and relationships between scientific articles [35]. VOSviewer analyzes keywords that frequently appear in articles and identifies the relevant ones. The results of the bibliometric mapping for keyword analysis with VOSviewer are presented in Figure 3.

Figure 3 was created using the VosViewer software, and a total of 48 relevant articles are saved in .ris format. Article files are inputted into VOSviewer, which is a mapping selected for co-occurrence words. The bibliometric mapping in Figure 3 shows that the co-occurrence of keywords consists of five clusters. These clusters indicate the link between "Spatio-Temporal Models" and "Climate." Forecasting climate is done chiefly with Spatio-Temporal Models and Time Series Models. In Figure 3, it can be seen that there are clusters that show climate variables that are often used by researchers, such as rainfall, Pacific Decadal Oscillation (PDO), atmospheric pressure, etc.

As revealed by an analysis of 48 relevant articles, the state-of-the-art in this field highlights significant progress in several key topics shown in Table 4. First, "GSTARIMA models" are emerging as a prominent approach to analyzing Spatio-Temporal data. This cutting-edge model combines the capabilities of time series analysis and spatial relationships, enabling a comprehensive understanding of complex interactions. Secondly, the exploration of "Heteroscedastic Error" in this study is in terms of overcoming the non-constant variance of errors in the GSTARIMA Model. By addressing these heteroscedastic errors, researchers aim to improve the accuracy and reliability of their forecasts, ultimately leading to more robust modeling results. In addition, "Kriging," a geostatistical interpolation technique, plays an essential role in spatial analysis. This method incorporates the estimation of unknown values based on observed values in the vicinity, incorporating spatial correlation. Collectively, these advances show the evolving research landscape in spatial-temporal analysis, featuring the integration of cutting-edge methodologies such as the GSTARIMA Model, the consideration of heteroscedastic errors, and the application of techniques such as Kriging to unravel complex spatial patterns and relationships.

Table 4 provides a comprehensive overview of the research developments related to the GSTARIMA/Spatio-Temporal model. Several studies have been conducted in Spatio-Temporal modeling while considering heteroscedastic errors. Kumar et al. (2022) used a STARMA-GARCH model to forecast monthly temperatures, resulting in minimal Mean Absolute Percentage Error (MAPE) values in their predictions [36]. Similarly, Monika et al. (2022) used the GSTARI-X-ARCH model to forecast rainfall influenced by humidity, showing favorable forecast accuracy [16]. In a different context, Akbar et al. (2020) introduced the GSTARMAX model to forecast air pollutants in Surabaya, achieving low Root Mean Square Error (RMSE) values [51]. Furthermore, the integration of Spatio-Temporal and Kriging models is seen in several articles. Dhaher et al. (2023) applied the Spatio-Temporal-Kriging model for temperature interpolation and prediction in Baghdad and Mosul cities [17]. Dai et al. (2022) used four methods, including LUR, LightGBM, ML, and Kriging, to forecast PM2.5 concentrations, which resulted in satisfactory R² accuracy [18]. Pramoedyo et al. (2020) adopted the GSTARX-SUR-Kriging model to forecast cocoa plant diseases affected by rainfall with reasonably accurate forecast results [54]. However, Abdullah et al. (2018) used the GSTAR-Kriging model to forecast rainfall in unobserved locations and produced fairly reliable prediction results [19]. Shu-qin et al. (2014) explored two different approaches, namely GWR and Kriging methods [69].

4. Discussion

4.1. GAP Analysis

Conducting a GAP analysis based on relevant articles illustrates the evolving research landscape in Spatio-Temporal modeling, heteroscedastic errors, and Kriging methodologies for forecasting climate and environmental data. These articles collectively represent vital insights and areas that need further exploration. The research that has been evaluated demonstrates a high propensity to use GSTARIMA models' capacity to forecast climate-related variables like temperature, precipitation, and air pollutants [19,37,40,49,50,51,64]. A common thread is the evaluation of model performance metrics, especially MAPE, RMSE, R², and MSE. However, the gap lies in comprehensively exploring complex parameter configurations in the GSTARIMA framework, especially in dynamic Spatio-Temporal systems. In addition, progress still needs to be made in validating these models using more sophisticated techniques, especially in handling larger and higher-dimensional data sets.

It is clear that heteroskedastic errors are critical to climate prediction, and special attention has been paid to using ARCH and GARCH models to address this issue [16,36,57]. Researchers concentrate on achieving higher prediction accuracy, indicated by lower RMSE and MSE values. However, there are differences in research in dealing with complex and non-linear forms of heteroscedasticity, which can arise from complex climate datasets. Identifying more flexible methods to handle this complexity could be an interesting subject of investigation. By incorporating Kriging into a spatiotemporal model, discernible trends can be identified, especially in interpolation and forecasting climate variables [17,18,19,47,58]. This study relies heavily on RMSE and MAE as tools to assess prediction accuracy. However, areas still need to be addressed in creating an adaptive Kriging method that can capture the temporal and spatial changes present in complex climate data. Due to these limitations, there is a possibility for more complex techniques that are adapted to changing trends and non-stationary data.

4.2. The Framework for model integration for climate forecast

4.2.1. The Integration of GSTARIMA Model with Heteroskedastic error and Kriging Method for forecasting

After the review of previous researchers and gap analysis, a conceptual integration model of GSTARIMA with Heteroskedastic error and the Kriging method is made to answer RQ1. The GSTARIMA model is processed following the Box-Jenkins method, including identification, parameter estimation, and diagnostic checking. The initial identification of the GSTARIMA model is to determine the stationarity of the data. If the stationary test results show that the data is not stationary, then a differencing process is carried out until stationary data is obtained. Next, check the order of the model univariately with the ARIMA Model. The model order is received from the results of the Autocorrelation Function (ACF) and Partial Autocorrelation Function (PACF) plots. Models with the same order are selected for further multivariate and Spatio-Temporal modeling. In terms of Spatio-Temporal modeling, a weight matrix is used that shows the diversity in locations. The order of the Spatio-Temporal Model is obtained based on the calculation of the Space-Time Autocorrelation Function (STACF) and Partial Space-Time Autocorrelation Function (STPACF).

Furthermore, parameter estimation for the GSTARI Model is carried out using the Ordinary Least Square (OLS) method. The error generated by the GSTARI Model is re-modeled to obtain the GSTIMA Model using the Maximum Likelihood (MLE) method. The GSTARI Model and GSTIMA Model are combined to produce the GSTARIMA Model. On the other hand, if exogenous variables influence the response variable, it becomes the GSTARIMA-X Model. Furthermore, predictions are made on the testing data for the GSTARIMA Model. The last stage of the model diagnostic check to determine the model error is white noise and homoscedasticity.

The GSTARIMA Model errors with heteroscedasticity errors are re-estimated following the ARCH/GARCH Model to overcome the non-constant variance of the errors. GSTARIMA Model errors are divided into mean equations and variance equations. The mean equation of the GSTARIMA Model error is estimated using the MLE method, and the variance equation is estimated using the GLS method. Integrating the GSTARIMA Model with the ARCH/GARCH Model can minimize the model error. This model is only able to do forecasting at locations that have observed values.

Regarding climate data, some areas do not have observation stations. The GSTARIMA and ARCH/GARCH models are then integrated with the Kriging method. The Kriging method is proven to forecast phenomena at unobserved locations. Estimated parameters in the GSTARIMA-ARCH Model are input to obtain parameters at unobserved locations. Furthermore, experimental and theoretical semivariogram calculations are carried out to obtain Kriging weights from unobserved locations. The estimated parameters for the unobserved locations are simulated to get the data at the unsampled locations. Finally, the data at unsampled locations are forecasted with the GSTARI-MA-ARCH Model. The integration of the GSTARIMA Model, ARCH/GARCH Model, and Kriging Method can forecast the phenomenon at unobserved locations in the future.

4.2.2. Data Analytics Life Cycle for Climate Forecasting

The conceptual Integrated Model of GSTARIMA, ARCH/GARCH, and Kriging Method is then used in forecasting climate that meets the criteria of Big Data. Regarding answering RQ2, the modeling flow follows the data analytics life cycle methodology presented in Figure 4. The initial stage begins with discovery, problem identification, determination of data sources to be processed, and hypotheses that are proven using theorems and mathematical formulas. The next step involves data preparation, inputting climate data into the process. Raw climate data is taken at a daily interval and cleaned to eliminate missing value data. Daily data is transformed by aggregating daily data into monthly data. In model planning, mathematical model integration is carried out. At this stage, the theorem that answers the research hypothesis is created. The GSTARIMA model is developed following the Box-Jenkins method. Integrating the GSTARIMA model, ARCH/GARCH, and Kriging method requires complex mathematical reasoning, especially in estimating model parameters. The integrated model is used in the Model Building stage with the input of data preparation results. Climate data is divided into training data and testing data. The results of forecasting are interpreted by the model obtained. Furthermore, visualization is carried out at the communication results stage, and recommendations are obtained. The last step is to operationalize the results of discoveries in Model development with theorems on mathematical modeling and dissemination.

5. Conclusions

In conclusion, a systematic literature review was conducted in developing the Integration GSTARIMA model with heteroscedastic error and the Kriging method for climate forecasting. A comprehensive search and analysis of the literature was performed to provide a clear understanding of the latest research. This research uses the PRISMA and bibliometric methods in analyzing the developments on this topic. In this paper, the results of the study in integrating the GSTARIMA Model with the ARCH/GARCH Model can overcome the problem of non-constant error variance. The GSTARIMA and ARCH models provide an overview of multivariate modeling affected by time, location, and non-stationary data. On the other hand, the GSTARIMA/Spatio-Temporal Model can only forecast at the observed location. Through the integration of the GSTARIMA Model with the Kriging Method, it has been discovered that the prediction of Spatio-Temporal phenomena becomes feasible for unobserved locations in the future. The development of the GSTARIMA Model, ARCH/GARCH, and Kriging Method allows the discovery of theorems in mathematical modeling. The application of this model to climate data uses the data analytics life cycle methodology for more detailed processing and more accurate information.