A Hybrid Space–Time Modelling Approach for Forecasting Monthly Temperature

Spatio-temporal forecasting has various applications in climate, transportation, geo-statistics, sociology, economics and in many other fields of study. The modelling of temperature and its forecasting is a challenging task due to spatial dependency of time series data and nonlinear in nature. To address these challenges, in this study we proposed hybrid Space–Time Autoregressive Moving Average-Generalized Autoregressive Conditional Heteroscedasticity (STARMA-GARCH) model in order to describe and identify the behaviour of monthly maximum temperature and temperature range in Bihar. At the modelling process of STARMA, spatial characteristics are incorporated into the model using a weight matrix based on great circle distance between the regions. The residuals from the fitted STARMA model have been tested for checking the behaviour of nonlinearity. Autoregressive Conditional Heteroscedasticity-Lagrange Multiplier (ARCH-LM) test has been carried out for the ARCH effect. The test results revealed that presence of both nonlinearity and ARCH effect. Hence, GARCH modelling is necessary. Therefore, the hybrid STARMA-GARCH model is used to capture the dynamics of monthly maximum temperature and temperature range. The results of the proposed hybrid STARMA11,0,0-GARCH0,1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm{STARMA} \left({1}_{1}, 0, 0\right)-\mathrm{GARCH} \left(0, 1\right)$$\end{document} model have better modelling efficiency and forecasting precision over STARMA11,0,0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm{STARMA }\left({1}_{1}, 0, 0\right)$$\end{document} model.


Introduction
Agricultural sectors are highly influenced by weather and climatic factors. It is sensitive to short-term changes in weather and long-term variations in climate. Weather parameters' variability in small regions, a few kilometres in size, is important in many hydrological, agricultural and energy contexts. Agricultural yields are directly and indirectly affected by changes in climatic factors such as rainfall, temperature, moisture and solar radiation, and the severity of extreme events like floods, droughts and wind storms. Due to the effect of global warming, there is a significant increase in global greenhouse gases and as a result there is steady rise in temperature of the planet earth. Therefore, mitigating climate change is one of the biggest challenges before humankind.
Temperature is the most important environmental factors at the earth surface. Forecasting of air temperature has been a crucial environmental factor required in many different areas such as agriculture, industry, energy, environment and tourism. Therefore, there is a need of systematic and scientific prediction of temperature for better planning, policy making and disaster management in the developing countries like India. There are various statistical methods used in forecasting monthly temperature. One of the most popular statistical methods used for forecasting temperature time series is Box-Jenkins' Autoregressive Integrated Moving Average (ARIMA). This model is suitable for only univariate and linear time series data. If the time series data set has non-linear and spatio-temporal relationship, then ARIMA model is not suitable.
The extension of univariate ARIMA time series models into the spatio-temporal domain results in a general class of Space-Time Autoregressive Integrated Moving Average (STARIMA) models [6,7]. It is used for the computation of linear dependencies between the variables in both time and space.
Application of spatial statistics mainly found in geo-statistics but in recent years the application of spatial statistics has been extended to sociology, economics, environmental, ecological and agricultural sciences. Mention may be mode of some recent works which show interest in spatio-temporal models, using both spatial and temporal information, for example Kamarianakis and Prastacos [12] work on traffic flow data set, Kurta and Tunay [13] developed STARMA model using of regional bank deposits data set, Zou et al. [26] modelling on wind power data set, Rathod et al. [22] modelling and forecasting on temperature data set, and Saha et al. [23] developed forecasting model on rainfall data set. STARMA models are special cases of Vector Autoregressive Moving Average (VARMA) models but VARMA models involve more number of parameters compared to STARMA models (Brynjarsdottir and Berliner [5] and Gupta and Robinson [11]). However, one main disadvantage of using STARMA model is that it cannot deal with nonlinear data series. To analyse data sets with nonlinearity, there are many statistical techniques such as ARCH technique first use by Engle [9], GARCH technique (Bolerslev [1]) and Threshold Autoregressive (TAR) technique (Tong and Lim [24]).
In this paper, an attempt has been made to develop a hybrid STARMA-GARCH model in order to investigate and predict spatio-temporal behaviour of monthly maximum temperature and temperature range. There are various hybrid statistical techniques available in the literature for modelling the rainfall data (Yusof et al. [25]), forecasting the exchange rate data (Pahlavani and Roshan [16]) and forecasting oil price data (Dritsaki [8]) but they all deal with only temporal data series. In real applications, weather data (rainfall, temperature, humidity etc.) constitute of dynamics involving spatial, temporal, linear and nonlinear features. Very few literatures have been found on spatio-temporal time series data which deal with nonlinear phenomena. Hence, this paper has made an attempt to develop a hybrid STARMA-GARCH model that can capture the spatio-temporal behaviour with linear and nonlinear dynamics of the data series under consideration.

Study Data
The monthly maximum temperature and temperature range of Bihar have been used here. The state of Bihar is categorized into four comprehensive agro-climatic zones, viz., zone I (North-West Alluvial Plains), zone II (North-East Alluvial Plains) and zone III (South Alluvial Plains). Further, zone III is divided into two parts zone III-A and zone III-B. For the present study, monthly maximum temperature and temperature range of eight districts of zone III-B, viz., Patna (PAT), Nalanda (NAL), Bhojpur (BHO), Rohtas (ROH), Gaya (GAY), Jahanabad (JAH), Aurangabad (AUR) and Kaimur (KAM), have been considered, which are presented in Fig. 1. The data set covers the period from January 1981 to December 2020. The data set on temperature has been collected from the NASA prediction of worldwide energy resources (https:// power. larc. nasa. gov/ data-access-viewer/) [15]. The data set from January 1981 to December 2019 and data set from January 2020 to December 2020 have been used for the purpose of model building and model validation, respectively.

STARMA Model
It is to be noted that the space-time data series is a vector time series. The m-dimensional VARMA process of order (p, q) is given by: where Z t is the N × 1 vector of observations at time t at N locations and t is the random error vector with t ∼ N(0, 2 I N ) independently distributed at each time t, and 2 I N is an N × N positive definite variance matrix. Further, 1 , … , p are the autoregressive coefficients and 1 , … , q are the moving average coefficients at all N × N matrices with 0 ≡ I . VARMA is not capable of modelling the spatial relationships and no assumptions of spatial stationarity can be made. Thus, incorporation of spatial dependence into a general VARMA model turns VARMA model into STARMA model. A STARMA model is expressed as: where Z t is the N × 1 vector of observations at time t at the N locations. p is the order of autoregressive, q is the order of moving-average, is the spatial order of the kth autoregressive term,m is the spatial order of the kth moving-average term, kl is the autoregressive parameter at temporal lag k and spatial lag l (scalar), kl is the moving average parameter at temporal lag k and spatial lag l (scalar), W l is the N × N spatial weight matrix with spatial order l and t is the random disturbances that are assumed to be normally distributed with , and E[ Z t prime t+s ] = 0, for s > 0. The STARMA ( p 1, … p , q m 1 …m q ) becomes a STAR model when q = 0, which is expressed as: when p = 0, it becomes a STMA model, which is expressed as: (1) The STARMA model-building procedure is similar to univariate ARIMA (Box and Jenkins [2]) model building procedure, viz., identification of the model, estimation of the parameters and diagnostic checking of the appropriateness of model for the observed data (Pfeifer and Deutsch [17][18][19][20][21]).

Preliminary Stage
Similar to ARIMA model, the first step is always to check for stationarity of time series data. If the time series data is nonstationary, then differencing should be done.

Identification
After making stationary of data series, next step is to examination of sample space-time autocorrelation (ST-ACF) and space-time partial autocorrelation functions (ST-PACF) with incorporate the space weight matrices in the analysis. The theoretical space-time auto variances can be shown in Table 1.

Estimation
After identification of proper STAR and STMA orders, next step is to estimate the parameters ( kl , kl ) of the STARMA model. The parameters are estimated by conditional maximum likelihood estimates (Pfeifer and Deutsch [18]). To find the maximum likelihood estimators, the normality assumption of residuals should be satisfied, i.e. t ∼ N(0, 2 I N ) . Estimation of the parameters by maximum  where: for t = 1, 2,... T, while setting t and Z t equal to zero for t < 1.

Diagnostic Checking
The fitted model is verified by first checking if the residuals are white noise for the best fitted model, the residuals should be white noise, i.e. t ∼ N(0, 2 I N ) . Autocorrelation of residuals is tested by the multivariate Box-Pierce test.

Construction of Spatial Weight Matrix
Prior to STARMA modelling, there is a need to define suitable spatial weight matrix (W l ) . Choice of the spatial weight matrix is very important because different spatial weight matrices often lead to different results. There are different ways to construct spatial weight matrices, including shared border and distance between two locations. The simplest form of weights is the binary weight, where two areas shared a common boundary gets weight 1 and otherwise gets 0 (Griffith [10]). In a spatial weight matrix, row normalization is a common practice (each row total is 1) and weights reflect a hierarchical order of spatial neighbours. The first order neighbours are closer than the second order ones, and the second order neighbours are closer than the third order neighbours and so on. In addition to that, each space unit is defined as its zero spatial order neighbour, and the zero order spatial weight matrix ( W 0 ) can be defined as unit matrix ( Table 2). In this article, a spatial weight matrix ( W 1 ) has been formed by taking weights that are inverse the distance between the two regions.
The elements are scaled such that: for each i. The distance ( d ij ) between the two geographical regions i and j is computed by using great circle distance formula of spherical law of cosines which is expressed as: and R = radius of the earth (6371 km). For each location, great circle distance ( d ij ) between the locations has been calculated using longitude and latitude coordinate (Table 3) which is represented in Table 4. Finally, the inverse distance weight matrix ( W 1 ) has been created through row normalization processes as defined in Table 5.

BDS Test
The BDS test is a nonparametric method of testing for nonlinear patterns in time series and proposed by Brock et al. [3, . The correlation integral at the embedding dimension m as: where T m = T − m + 1 , and I Z m t , Z m s is an indicator function of the event which equal to 1 if the sup norm ‖Z m t − Z m s ‖ < and equals to zero otherwise. Thus, the correlation integral measures the probability that any two embedding dimension m points are within a distance of distance of each other, If the {Z t } are iid. The probability should be equal to the following in the limiting case: The BDS statistic is then defined for a time series of length T as: where T = sample size, = arbitrarily chosen proximity parameter and m,T ( ) = standard sample deviation. Under fairly moderate regularity conditions, BDS m,T ( ) converges to the standard normal distribution. Decision: If the value of | | BDS m,T ( ) | | > 1.96, then the null hypothesis will be rejected, otherwise accepted at 5% level of significance.

ARCH-LM Test
GARCH modelling starts by identifying whether or not the space-time series data contains heteroscedasticity. The residuals obtained from the data are checked by using ARCH-LM test to know the ARCH effect [9].
The steps for ARCH-LM test: I. Define the linear regression: where t denotes the error term and q is the pre-specified positive integer. II. The test statistic is: where: and T is total number of observation, x i = 2 t i . R 2 is the R-squared from the regression of 2 t in Eq. (13). Here R 2 follows a χ 2 (q) distribution. III. The null and alternative hypotheses are:

GARCH Model
If there is ARCH effect, then the next step involves in identifying the appropriate orders for the GARCH (p, q) model. The GARCH model is a generalized form of ARCH model. The GARCH (p, q) model in which conditional variance is modelled as:

Propose Hybrid STARMA-GARCH Model
There are two phases in the hybrid STARMA-GARCH modelling. In the first phase, a STARMA model is fitted on stationary and linear space-time data and in the second phase, residuals of the STARMA model are fed as an input to the GARCH model. This hybrid model, which combines STARMA and GARCH model containing non-linear residuals patterns, is applied to forecast the air temperature data.

Diagnostic Checking of the Hybrid STARMA-GARCH Model
The diagnostic tests of hybrid STARMA-GARCH models are based on residuals. Normality of the residuals is tested using Ljung and Box Q-statistics (Ljung and Box [14]) for the serial correlation test. Test hypotheses are given below: H 0 : no serial correlation in the residuals H 1 : there is serial correlation The Q-statistic is given as: where N is the sample size, L is the autocorrelation lags and ̂ 2 j is the square of the estimated autocorrelation at lag j. Decision: If the value of Q > 1.96, then the null hypothesis will be rejected at 5% level of significance.

Forecast Evaluation
Mean absolute percentage error (MAPE) has been used to evaluate the forecast efficiency for which the formula is given below: where n = total numbers of predicted value, Z t = actual value at time t and Ẑ t = corresponding predicted value.

Preliminary Stage
Before proceeding to use STARMA modelling, it is important to check the stationary of the time series data. The results of the augmented Dickey-Fuller test showed that spatio-temporal data series of all locations satisfy stationarity. The time series plots of data set (Figs. 2 and 3) indicate the presence of seasonality which needs to be removed before analysis. Therefore, the original time series data are seasonally adjusted by seasonal differencing with period 12. Finally, seasonally adjusted data series of each location become stationary and have been used for analysis.

Identification
The model identification is a very crucial step for spatiotemporal modelling. First, the presence of correlation in space-time data is checked. The values of chi-square of multivariate Box-Pierce test for maximum temperature and temperature range are 3346.149 and 2824.946 with p value < 0.001, respectively, which tell that reject the null hypothesis (H 0 : non-correlation) at 1% level of significance. After checking the presence of space-time correlation, the suitable candidate STARMA model has been selected with help of ST-ACF and ST-PACF plots. The ST-ACF and ST-PACF plots of original data series have shown the presence of seasonality. Therefore, seasonally adjusted data of ST-ACF and ST-PACF plots are used to determine STARMA model (Figs. 4 and 5). Here STAR 1 1 and STMA(0) have been selected for both the maximum temperature and temperature range. Therefore, the STARMA 1 1 , 0, 0 model has been identified. This model is also called STAR 1 1 model.

Estimation
Parameters of the fitted STARMA 1 1 , 0, 0 model are estimated by maximum likelihood method. The estimated parameters and its standard error are presented in Table 6. All the estimated parameters are statistically significant at 5% level of significance. The STARMA (1 1 , 0, 0) model can be expressed as: where Eqs. (18) and (19) represent the STARMA models for temperature range and maximum temperature, respectively.

Fig. 2 Time series plots for monthly maximum temperature
Values in the parenthesis indicates the standard error; * and ** are the significant at 5% and 1% level of significance, respectively.

Diagnostic Checking
After parameters estimation, diagnostic checking of residuals is applied by using multivariate Box-Pierce non-correlation test. The test applied to the residuals give the chi-square statistics values as 390.25 and 455.12 with p values < 0.001 for maximum temperature and temperature range, respectively, which tell that reject the null hypothesis (H 0 : non-correlation) at 1% level of significance. Presence of spatio-temporal correlation within residuals also be seen in the ST-ACF and ST-PACF plots of residuals (Figs. 6 and 7).

BDS and ARCH-LM Test for Residuals
The BSD test results show the presence of both the linear and nonlinear patterns in the correlated residuals of maximum temperature and temperature range. The ARCH-LM test of square residuals is suffering from serial correlation in the most of the locations. Hence, it is necessary to develop a best fitted model for analysis of monthly temperature. A GARCH modelling is necessary to handle heteroscedasticity in the residuals series.

Hybrid STARMA-GARCH Model
The STARMA-GARCH model is a combination of two models, first one is the STARMA model that take into account of mean behaviour of a space-time series and the second one is the GARCH model that used to model the variance of the series. In this hybrid model, residuals of the fitted STARMA model are used for suitable GARCH modelling. The results of minimum

Comparative Analysis
To evaluate the forecast efficiency, MAPE has been computed for the forecasting performance for all eight locations. The MAPE values under training data set for maximum temperature and temperature range are given in Table 9. According to the results of the MAPE, the hybrid model has lowest mean absolute percentage errors compared to STARMA model for

Conclusions
In this study, the problem is to forecast the monthly maximum temperature and temperature range as a spatio-temporal forecasting problem for III-B agro climatic zone of Bihar.
In the real-world phenomena, weather parameters work in a very complex manner such they can vary over small regions and for short periods of time which lead to linear, as well as nonlinear patterns. In these variations, there are various statistical techniques used in forecasting of weather parameters. One of the most popular statistical methods used for forecasting temperature time series is Box-Jenkins ARIMA but this model is suitable for only univariate and linear time series data. If the multivariate spatio-temporal time series data set shows non-linear pattern, then ARIMA model is not suitable for its modelling. In order to overcome the drawback of traditionally used models, hybrid model combining the STARMA and GARCH has been investigated. Based on the results obtained by the STARMA model and hybrid model for both the variables, viz., maximum temperature and temperature range, it can be seen that fitted hybrid model performed better as compared to STARMA model for the training data set. After validation of proposed hybrid model, the fitted hybrid model has lower forecast errors than the STARMA model. Hence, it can be concluded that the spatio-temporal hybrid model has better forecasting accuracy because of the inclusion of both spatial information and nonlinear patterns in the temperature data sets. The proposed hybrid model also be used to model variations is other weather parameters, which will be taken up in a forthcoming paper.
Author Contribution All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Ravi Ranjan Kumar, Kader Ali Sarkar, Debasis Bhattacharya and Digvijay Singh Dhakre. The first draft of the manuscript was written by Ravi Ranjan Kumar. Kader Ali Sarkar, Debasis Bhattacharya and Digvijay Singh Dhakre commented on its improvement. All authors read and approved the final manuscript.
Data Availability All data are available online as provided in [15].

Code Availability
The codes written for the current study are available on reasonable request from the corresponding author.

Declarations
Ethics Approval The authors confirm that this work is original and has not been published elsewhere nor is it currently under consideration for publication elsewhere.

Consent to Participate Not applicable.
Consent for Publication Not applicable.