Data-driven models for atmospheric air temperature forecasting at a continental climate region

RESEARCH ARTICLE Data-driven models for atmospheric air temperature forecasting at a continental climate region Mohamed Khalid AlomarID1, Faidhalrahman Khaleel1, Mustafa M. Aljumaily1, Adil Masood2, Siti Fatin Mohd Razali3, Mohammed Abdulhakim AlSaadi4, Nadhir Al- AnsariID5*, Mohammed Majeed HameedID1,3* 1 Department of Civil Engineering, Al-Maarif University College, Ramadi, Iraq, 2 Department of Civil Engineering, Jamia Millia Islamia, New Delhi, India, 3 Department of Civil Engineering, Faculty of Engineering and Built Environment, Universiti Kebangsaan Malaysia, Selangor, Malaysia, 4 Natural and Medical Sciences Research Center, University of Nizwa, Nizwa, Oman, 5 Civil Engineering Department, Environmental and Natural Resources Engineering, Lulea University of Technology, Lulea, Sweden * m.majeed@uoa.edu.iq, mohmmag1@gmail.com (MMH); nadhir.alansari@ltu.se (NAA) Abstract Atmospheric air temperature is the most crucial metrological parameter. Despite its influ- ence on multiple fields such as hydrology, the environment, irrigation, and agriculture, this parameter describes climate change and global warming quite well. Thus, accurate and timely air temperature forecasting is essential because it provides more important informa- tion that can be relied on for future planning. In this study, four Data-Driven Approaches, Support Vector Regression (SVR), Regression Tree (RT), Quantile Regression Tree (QRT), ARIMA, Random Forest (RF), and Gradient Boosting Regression (GBR), have been applied

PLOS ONE | https://doi.org/10.1371/journal.pone.0277079 November 3, 2022 1 / 31 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 OPEN ACCESS Citation: Alomar MK, Khaleel F, Aljumaily MM, Masood A, Razali SFM, AlSaadi MA, et al. (2022) Data-driven models for atmospheric air temperature forecasting at a continental climate region. PLoS ONE 17(11): e0277079. https://doi. org/10.1371/journal.pone.0277079 Editor: Yangyang Xu, Texas A&M University, UNITED STATES Received: April 10, 2022 Accepted: October 18, 2022 Published: November 3, 2022 Peer Review History: PLOS recognizes the benefits of transparency in the peer review process; therefore, we enable the publication of all of the content of peer review and author responses alongside final, published articles. The editorial history of this article is available here: https://doi.org/10.1371/journal.pone.0277079 Copyright: © 2022 Alomar et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability Statement: The data are within the paper. Funding: The authors received no specific funding for this work.

1. Introduction It is well-known that numerous meteorological and ecological events, human life, and crops in agricultural areas are significantly influenced by climate conditions as well as several factors related to the environment’s physical conditions. The natural resources that provide humans basic needs and opportunities for social and economic development are part of the physical environment, including land, air, and water. A clean and healthy environment is one of the essential principles that should be preserved and protected [1]. The temperature parameter is seen as the most influential parameter out of all meteorological parameters, which reflects the effect of climate change on earth and its surrounding atmosphere. Recently, climate change has caused extreme natural phenomena such as heat waves, severe winters, heavy snowfall, and droughts worldwide, leading to environmental and health crises [2–6]. Air temperature prediction helps meteorologists to know the likelihood of hurricanes and floods in an area [7]. Various meteorological parameters such as rainfall, humidity, atmospheric pressure, wind speed, solar energy, and soil temperature are significantly correlated to air temperature [8]. Moreover, air temperature is one of the most influential factors in evapotranspiration, which is vital for managing water resources and agricultural activities [9]. Accurate air temperature prediction is substantial in many decision-making sectors, such as energy, agriculture, trans- portation, and tourism [10]. Additionally, accurately predicting air temperature is the most

crucial aspect of environmental studies involving operational eco-environmental systems. From the industrial aspect, predicting air temperature is essential in energy management strat- egies to obtain comfortable indoor temperatures and eventually reduce the consumption of energy [11]. According to the literature, two main approaches have been used to predict air temperature: general circulation (GCM) and statistical models. The GCM is utilized to comprehend the dynamics behind climate system physical components, derive global temporal and spatial changes and make predictions based on the future forcing of greenhouse gasses and aerosols [12]. GCMs can be applied to the problem of attributing climate change from a season to a decade ahead. Conversely, statistical models attempt to determine whether climate change is externally driven by minimizing the utilization of complex climate models. They are generally more straightforward and less computationally intensive than the GCMs, and several studies have showed that the use of statistical models has produced results consistent with GCMs. Var- ious statistically based approaches have been proposed recently, several of which have been developed in the econometric literature. The statistical models can be categorized into two approaches: cointegration approaches which determine the relationship between non-station- ary and stationary times series [13], and regression approaches which evaluate the characteris- tics of time series for a given temperature data [14, 15]. However, since temperature prediction

involves high nonlinearity and dimensionality, the statistical models faced some drawbacks in capturing them [16]. Meanwhile, machine learning (ML) approaches have attracted much attention due to their superlative performance in dealing with high nonlinearity phenomena [17, 18] and solving complex problems such as drought [19–24], rainfall [25–29], evapotranspiration [30–34] and streamflow [35–38]. For example, a study was conducted in the Queensland area where ML models’ performances were compared with the Australian Predicted Ocean-Atmosphere Model (POAMA) for precipitation prediction. The POAMA model showed a significant improvement in the predictive performance of the ML modeling framework. It was reported that the performance of the neural network (NN) model was superior to POAMA in precipita- tion prediction over three regions in Queensland [39–41]. For temperature prediction, A. Sekertekin et al. [6] used the adaptive neuro-fuzzy inference system (ANFIS) and long-short PLOS ONE Air temperature forecasting: Mutli-time scale PLOS ONE | https://doi.org/10.1371/journal.pone.0277079 November 3, 2022 2 / 31 Competing interests: The authors declare that there are no conflicts of interest.

term memory (LSTM) network to predict temperature for both ultra short-term and short term period(hourly and one day ahead). The results showed that the LSTM model was able to efficiently predict the temperature for both the time scales. However, the LSTM has several dis- advantages, such as it requires longer time and more memory to train. Besides, its parameters are difficult to assign and implement and the outcomes are vulnerable to various random weight initializations. S. Salcedo-Sanz et al. [12] used the support vector regression (SVR) and multi-layer perceptron (MLP) models to predict the mean monthly air temperature. The data- set from the monitoring stations located in New Zealand and Australia was used for the model development. The results showed that the SVR model provided the best accuracy in tempera- ture prediction. Overall, very few studies based on daily temperature prediction have been conducted for regions with continental climatic conditions Therefore, the main objective of this study is to forecast air temperature over a continental climate case study which is in North America. Two-time scales are adopted in this study, daily and weekly. For fulfilling this task, four Data-driven models i.e., Support Vector Regression (SVR), Regression Tree (RT), Quan- tile Regression Tree (QRT), and Gradient Boosting Regression (GBR) have been applied. These models have been used to predict one-day and one-week temperature ahead depending on the past temperature values for both time scales (weekly and daily). Comprehensive com-

parisons supported by statistical measures and comparative figures have been applied to select the most efficient models. 2. Methodology 2.1 Case study North Dakota is located in the middle of North America and is subjected to extreme climate conditions, with hot summers and cold winters. Due to its inland location and proximity to both the North Pole and the Equator, which are almost equal, there are noticeable temperature fluctuations. Furthermore, it has been observed that the temperature varies extremely from season to season, which may be responsible for the changes in weather throughout the time [42]. Since North Dakota has a continental climate, forecasting the patterns of meteorological parameters is a challenging task. The difficulty in simulating weather parameters in such region may be due to the nature of the fluctuating climate during the seasons. Tables 1 and 2 show the statistical characteristics of the minimum, mean, average, standard deviation, and skewness of the daily and weekly air temperature values at the Crary Table 1. Statistical characteristics of Crary station: Daily scale. Month/statics Min Average Max St. D Skewness Jan -35.50 -13.018 4.12 7.905 -0.194 Feb -32.55 -12.847 5.70 7.695 0.015 Mar -28.19 -4.191 14.85 7.465 -0.349 Apr -11.80 4.542 18.38 5.919 -0.320 May -1.67 11.860 24.87 5.072 -0.008 Jun 5.15 17.944 30.56 3.626 -0.101 Jul 12.28 20.799 27.65 2.916 -0.340 Aug 9.19 19.408 28.28 3.248 0.147 Sep 0.66 14.598 27.89 4.555 -0.078 Oct -8.15 6.098 22.86 5.436 0.137

Nov -19.79 -2.479 12.81 6.430 -0.165 Dec -30.83 -10.901 5.00 7.147 -0.320 https://doi.org/10.1371/journal.pone.0277079.t001 PLOS ONE Air temperature forecasting: Mutli-time scale PLOS ONE | https://doi.org/10.1371/journal.pone.0277079 November 3, 2022 3 / 31

meteorological station from 2000 to 2021. According to reported data, the three hottest months are June (17.944 C˚), July (20.799 C˚), and August (19.408 C˚) while the coldest months in this case study are December (-10.901 C˚), January (-13.018 C˚), and February (-12.847 C˚). Furthermore, the recorded air temperatures have extreme values (far from the mean) in four months (i.e., December to March). In these months, the standard deviation val- ues of the data are very high compared to other months. Nevertheless, the dispersion of the data through June (St. D = 3.626), July (St. D = 2.916), and August (St. D = 3.248) is very little, which means that the data are more consistent with their normal rates. Notably, the utilized data in this study are collected from the open-source website of the Crary station [43]. Finally, the location of the studied region is shown in Fig 1. 2.2 Support vector regression Support vector regression (SVR) is considered a powerful and efficient tool based on the notion of statistical learning and was first introduced by Vapnik [44] to describe regression as a part of the support vector machine (SVM). Based on the principle of structural risk minimi- zation (SRM), SVR has been successfully implemented in real-world challenge modeling by overcoming classification and regression tasks [45]. The linear relationship between indepen- dent variables (x1, x2, x3, � � �, xr) and dependent variable (y) is given in the equation below. y ¼ f x ð Þ ¼ w; x ð Þ þ b ð1Þ Where wi and b are the weight and bias of the model, respectively. ;(x) is the higher dimen-

sional feature space converted from the independent vector (input). These parameters can be determined by minimizing kwk2 = (w.w) as follows min 1 2 k w2 k þC X r i¼1 xi þ xi � ð Þ ð2Þ under the constraints yi w; x ð Þ b � ε þ xi yi w; x ð Þ b � ε xi � 8i xi ; xi � � 0 8 > < > : 2 1; . . . n f g ð3Þ Table 2. Statistical characteristics of Crary station: Weekly scale. Month/statics Min Average Max St. D Skewness Jan -28.23 -13.112 -1.84 5.882 -0.070 Feb -27.35 -13.191 0.26 6.093 0.111 Mar -20.56 -5.287 11.32 6.311 -0.115 Apr -10.33 3.833 12.74 4.926 -0.383 May 3.22 11.185 21.55 3.879 0.101 Jun 10.28 17.461 23.22 2.743 -0.133 Jul 15.76 20.644 24.47 1.949 -0.384 Aug 13.98 19.650 25.50 2.292 -0.056 Sep 6.20 15.149 22.96 3.319 -0.229 Oct -3.89 7.212 18.09 4.414 0.036 Nov -14.03 -1.822 8.59 5.210 -0.133 Dec -22.77 -9.993 1.34 5.816 -0.380 https://doi.org/10.1371/journal.pone.0277079.t002 PLOS ONE Air temperature forecasting: Mutli-time scale PLOS ONE | https://doi.org/10.1371/journal.pone.0277079 November 3, 2022 4 / 31

Where, C is the regularization constant, ξi, ξi� are the slack variables and ε is the size of the tube, "denoting the accuracy of the function to be approximated" [46]. Based on Lagrange multipliers, the standard SVR can overcome the following optimization problem. min 1 2 Xr ijðai a� i Þ aj a� j � � xi; xj � � Xr i¼1 ai þ a� i � þ Xr i¼1 gi ai a� i � ¼ 0 ð4Þ Fig 1. The location of Crary station. Source (https://d-maps.com/index.php?lang=en; https://www.usgs.gov/). https://doi.org/10.1371/journal.pone.0277079.g001 PLOS ONE Air temperature forecasting: Mutli-time scale PLOS ONE | https://doi.org/10.1371/journal.pone.0277079 November 3, 2022 5 / 31

Which is subjected to: X r i¼1 ai a� i � ¼ 0 ð5Þ ai; a� i 2 0; r ½ � ð6Þ Where ρ are the cost factor, ai; a� i � 0 are the Lagrange multiplier factors. The linear SVR can be written as follows f x ð Þ ¼ X r i¼1 ai a� i � xi; x ð Þ þ b ð7Þ This equation may be considered inappropriate for solving many engineering problems because of its linear characteristic, while engineering problems often need a non-linear regres- sion analysis. Therefore, in order to switch the input data to a much higher-dimensional space, nonlinear kernel functions are utilized. In this regard, the radial bases kernel function (RBF) is used in this study and can be expressed mathematically in the equation below. K xi; x ð Þ ¼ exp b k x xi k 2 ð Þ ð8Þ Where K(xi, x) represents the kernel function and β is the bandwidth of K(xi, x). 2.3 Regression tree and quantile regression tree A decision tree (DT) is a supervised machine learning-based technique that uses labeled data (data with known target attributes) to carry out simulations with the help of classification and regression algorithms [47]. In general, DT’s consist of three types of nodes: decision Root nodes, internal nodes, and leaf nodes, where each node or leaf denotes a class label while the branches denote the outcome of the test performed [48]. The technique splits the input dataset on the basis of the most significant splitter or differentiator in the input variables. This process of data division and selection of the most significant attribute in the dataset is governed by the

classification and regression algorithms. The technique follows a top-down approach as the top portion holds all the observations at one spot, which splits into two or more branches that further split. This approach is also referred to as the greedy approach, as it only incorporates the current nodes without focusing on the future nodes [49]. The decision tree algorithm con- tinues to run until a stop criterion such as the minimum number of observations etc., is attained. Once this criterion is achieved and a decision tree is developed, many nodes are detected as outliers which may be addressed through the tree pruning method. This, in turn, improves the forecasting accuracy of the DT-based model. In the same method that regression minimizes cost function (i.e., squared-error loss) when forecasting a single point estimate, the quantile regression tree (QRT) minimizes the loss func- tion in forecasting a particular quantile. The median, or 50th percentile, is the most commonly used quantile, and the quantile loss is just the sum of absolute errors in this case. Additionally, quantiles can be used as endpoints of prediction intervals; for instance, the 10th and 90th per- centiles define an 80-percentile range in the middle. It appears that the quantile loss differs according to the evaluated quantile, with higher quantiles penalizing more for negative errors and lower quantiles penalizing more for positive errors. Accordingly, in this study, we used the median (the 50th percentile), which is the most well-known quantile.

In the fields of artificial intelligence and search algorithms, pruning is a data compression method used to minimize the size of decision trees by deleting parts of the tree that are deemed PLOS ONE Air temperature forecasting: Mutli-time scale PLOS ONE | https://doi.org/10.1371/journal.pone.0277079 November 3, 2022 6 / 31

non-critical and repetitive to the regression of instances. By assessing the predictive value of each node in a regression tree, regression tree pruning decreases the danger of overfitting. Nodes that do not increase the anticipated prediction efficiency on new data are substituted with leaves. 2.4 Gradient Boosting Regression (GBR) GBR is an ensemble machine learning approach that enhances the prediction performance of a classical decision tree by incorporating a sequential statistical process called boosting, of which the principle idea is to combine a set of weak predictive models to form a single and high accurate predictive model [50, 51]. The technique applies an iterative procedure, where the estimates of the new tree model (weak learner) are updated with the pseudo residuals (neg- ative gradient of the loss function) of the current learner [52]. This process is repeated until the loss function of the model is reduced to a minimum value, thus improving the forecasting performance of the model. The iterative training process of the GBRT with K decision trees can be briefly explained as follows: For a given training dataset D = {(x1, y1), (x2, y2),. . .. . ., (xn, yn)}, the loss function is com- puted as: L y; f x ð Þ ð Þ ¼ y f x ð Þ ð Þ2 ð9Þ Step 1. Initialize the new tree model (weak learner) with a constant value: f0 x ð Þ ¼ arg minc XN i¼1 L yi; c ð Þ ð10Þ Step 2. Assume the number of iterations m = 1, 2, 3. . .. . .., K (a) For i = 1,2,3. . .. . ..,N.The pseudo residuals of the ith training data is calculated as: rmi ¼ @L y; f xi ð Þ ð Þ @f xi ð Þ � �

f xð Þ¼fm 1 xð Þ ð11Þ (b) Fit a regression tree in terms of rmi, and deduce the leaf node area Rml of the mth tree. Predict the leaf node area of the decision tree to attain an approximate value of the fitting residual. (c) For l = 1,2,3. . .. . .., L. Adopt linear search to attain the value in the leaf node range and minimize the loss function with gradient descent. The best residual fitting value of each blade is as follows: cm ¼ arg minc XN i¼1 L yi; fm 1 xi ð Þ þ c ð Þ ð12Þ (d) Update the regression tree fm x ð Þ ¼ fm 1 x ð Þ þ XL l¼1 cmlI x 2 Rml ð Þ ð13Þ Step 3. Obtain the final model f x ð Þ ¼ fM x ð Þ ¼ XM m¼1 XL l¼1 cmlI x 2 Rml ð Þ ð14Þ PLOS ONE Air temperature forecasting: Mutli-time scale PLOS ONE | https://doi.org/10.1371/journal.pone.0277079 November 3, 2022 7 / 31

2.5 Autoregressive model The autoregressive integrated moving average (ARIMA) is a historical data-based model and is considered the most common time series modeling approach first introduced by Box and Jenkin in 1976 [53]. ARIMA model is considered a hybrid model in which the Autoregressive (AR) and moving average (MA) models generalized forms are combined for modeling non- stationary univariate time series data by approximating the time series using a mathematical model based on past and current values. The utilization of ARIMA model is done by setting the order of three terms: autoregressive, sequence difference, and moving average. The general successive difference equation for dth order can be mathematically expressed as follows. DdTt ¼ 1 B ð ÞdTt ð15Þ Where d is the order of sequence difference and B is the backshift operator. The general ARIMA equation can be briefly presented as follows [54]. ;p B ð Þwt ¼ yq B ð Þet ð16Þ Where ;p(B) is the autoregressive operator of order p, θq(B) is the moving average term of order q and wt = ΔTt. 2.6 Random forest Random forest is a supervised machine learning technique which is made up of large number of small decision trees, known as estimators, which generate their own predictions. ’Forest’ generated by the random forest algorithm is trained through bagging or bootstrap aggregating [55]. Bagging is an ensemble meta-algorithm that fine-tunes the prediction accuracy of machine learning algorithms. The (random forest) algorithm produces the output based on

the predictions of the decision trees. It predicts by taking the average or mean of the output from various trees. Increasing the number of trees increases the precision of the outcome. The various advantages of this technique over other machine learning approaches such as need of less computation time, ease of working with high-dimensional data, strong fault tolerance and parallel processing make it suitable even for very high-dimensional problems like air tempera- ture forecasting. 2.7 Model development In this work, four regression models i.e., SVR, GBR, QRT, and RT have been used to predict the daily and weekly air temperatures over the continental climatic region of North America. The time-series data was collected from the Crary meteorological station from 2000 to 2021. For selecting the best input lags, the autocorrelation function (ACF) and partial autocorrela- tion function (PACF) have been used to analyze the data. According to Fig 2, the ACF provides more information on the time series properties like stationary, trend pattern, seasonality, and randomness. The daily and weekly temperature patterns were examined to determine the most appropriate predictors utilizing correlation statistics such as ACF and PACF, respec- tively. The statistical techniques used the time-lagged data from temperature time series to esti- mate the daily and weekly intervals between the present T value and prior T value for any given observation (i.e., a time lag) [55]. Thus, selecting which lags have a significant correlation

and significant information may benefit. Besides, the lags confined between upper and lower bounds are neglected because they have lower correlations and represent the white noise in the time series, which cannot be predicted. Both ACF and PACF are provided in the following PLOS ONE Air temperature forecasting: Mutli-time scale PLOS ONE | https://doi.org/10.1371/journal.pone.0277079 November 3, 2022 8 / 31

Fig 2. ACF and PACF for input determination. https://doi.org/10.1371/journal.pone.0277079.g002 PLOS ONE Air temperature forecasting: Mutli-time scale PLOS ONE | https://doi.org/10.1371/journal.pone.0277079 November 3, 2022 9 / 31

equations: ACF ¼ PN k t¼1 Xt �x ð Þ Xtþk �X � PN k t¼1 Xt �X ð Þ ð17Þ PACF ¼ ACF Pk 1 j¼1 PACFk 1;jACFk 1 1 Pk 1 j¼1 PACFk 1;jACFk 1 ð18Þ The N is the total number of temperature records, Xt and X is the time series record at time t and the mean of the temperature records, and finally, the k is number of lags in the time series data. lower and upper limits (UP, LO) can be determined at 95% significant level by the following equation: LO; UP ¼ � 1:96 ffiffiffiffi N p ð19Þ The ACF declines more slowly concerning the daily scale, which means that the time series data is not stationary. It is challenging to select the most effective lags for daily scale using ACF because of seasonality. The PACF, similar to the ACF, shows the association between two rec- ords that the shorter delays between those observations do not describe. For instance, the par- tial autocorrelation coefficient for the third lag in the daily scale temperature is only a correlation that the previous short lags (lag two and lag one) have not explicitly explained. Therefore, the PACF is more suitable for selecting the input lags for predicting the short scale of time series than the ACF. Table 3 shows the input combinations used in this study. It is worth mentioning that 70% of the data is used to train the suggested models, and 30% of the data included the end of the time series of the data utilized for the testing phase and checking the models’ performances (see Fig 3). The following steps summarize the primary process of developing the models for forecasting short-and mid-term air temperature.

Table 3. Air temperature input design for Baker and Crary stations. Model Input groups Output Scale Training data records Testing data records M1 Tt−1 Tt Daily 5624 2411 M2 Tt−1, Tt−3 Tt Daily 5622 2411 M3 Tt−1, Tt−3, Tt−4 Tt Daily 5621 2411 M4 Tt−1, Tt−3, Tt−4, Tt−5 Tt Daily 5620 2411 M5 Tt−1, Tt−3, Tt−4, Tt−5, Tt−6 Tt Daily 5619 2411 M6 Tt−1, Tt−3, Tt−4, Tt−5, Tt−6, Tt−7 Tt Daily 5618 2411 M7 Tt−1, Tt−3, Tt−4, Tt−5, Tt−6, Tt−7, Tt−8 Tt Daily 5617 2411 M8 Tt−1, Tt−3, Tt−4, Tt−5, Tt−6, Tt−7, Tt−8, Tt−9 Tt Daily 5616 2411 M9 Tt−1, Tt−3, Tt−4, Tt−5, Tt−6, Tt−7, Tt−8, Tt−9, Tt−10 Tt Daily 5615 2411 M10 TWt−1 TWt Weekly 803 344 M11 TWt−1, TWt−2 TWt Weekly 802 344 M12 TWt−1, TWt−2, TWt−3 TWt Weekly 801 344 M13 TWt−1, TWt−2, TWt−3, TWt−4 TWt Weekly 800 344 M14 TWt−1, TWt−2, TWt−3, TWt−4, TWt−5 TWt Weekly 799 344 M15 TWt−1, TWt−2, TWt−3, TWt−4, TWt−5, TWt−6 TWt Weekly 798 344 M16 TWt−1, TWt−2, TWt−3, TWt−4, TWt−5, TWt−6, TWt−7 TWt Weekly 797 344 M17 TWt−1, TWt−2, TWt−3, TWt−4, TWt−5, TWt−6, TWt−7, TWt−8 TWt Weekly 796 344 https://doi.org/10.1371/journal.pone.0277079.t003 PLOS ONE Air temperature forecasting: Mutli-time scale PLOS ONE | https://doi.org/10.1371/journal.pone.0277079 November 3, 2022 10 / 31

1. Collecting the daily data for air temperature from a continental climate. In this study, Crary station is selected. 2. Converting the temperature values from Fahrenheit to Celsius using the following formula. TC ¼ 5 9 ðTF 32Þ ð20Þ TC and TF are temperatures measured in degrees Celsius and Fahrenheit. 3. Computing the mean weekly temperatures. 4. Selecting the best lags using ACF, and PACF for both scales (weekly and daily). Fig 3. Daily and hourly air temperature measured in Crary station over the period from 2000 to 2021. https://doi.org/10.1371/journal.pone.0277079.g003 PLOS ONE Air temperature forecasting: Mutli-time scale PLOS ONE | https://doi.org/10.1371/journal.pone.0277079 November 3, 2022 11 / 31

5. Data partition: The time-series data is relatively long (2000 to 2021), 70% of the data are used for model calibration (2000 to 2015), and the rest are used for testing. The number of data points was fixed in that case to ensure a fair evaluation of the proposed model through- out the most critical step (testing phase). Notably, this procedure does not affect the data partitions. For example, for daily scale 2411, which represents about 30% of the entire daily records (Table 3) 6. Normalizing the training and testing dataset based on the minimum and maximum tem- perature in the training data set using the following formula [56]: Tnormalized ¼ Ti Tmin Tmax Tmin ð21Þ Tnormalized is the normalized temperature for ith temperature record (Ti) while Tmin and Tmax are minimum and maximum and temperature data obtained from the training data set. 7. Assigning the hyperparameters of the applied models. The trial- and error- method is used for this process where each model was trained 100 times over the training dataset with dif- ferent parameters. When these models were trained several times, the best ones were selected according to the statistical criteria. According to several statistical metrics, the model which generates lowest forecasting error in the training step is selected. In addition, the performance of the model should be stable so that there is no significant difference between its performance in the training and testing phase. 8. De-normalizing the data based on the following formula: Ti ¼Tmin þ Tnoralizei Tmax Tmin ð Þ � � ð22Þ

9. Evaluating the accuracy of the models with the testing dataset. It is important to mention that all models are constructed using MATLAB software. The candidate parameters of the applied models can be illustrated below: • RF: The number of trees is selected between 20–100 while the leaf node ranges from 1–5. • GBR: the learning rate 2 [0 1] and the number of trees 2 [150,1]. In Bag fraction = 1. • SVR: Box Constraints "regularization parameter "is set between 0.7 to 1. mean the sigma ranges from 0.8452 to 0.7071. • Epsilon parameter 2 [0.6 1] and sigma. Finally, the kernel scale parameter ranges from 0.8 to 1. • DT: MaxNumSplits (maximum number of decision splits) 2 [1 8]. Tree depth controllers 2 [5 10]. • Bag fraction = 1 for assembling models that’s mean "roughly 2/3 of input data is selected for training for every tree and the remaining 1/3 is used as out-of-bag observations". In the second scenario of this work, we used the randomization method to divide the data into the training and testing phase. In this scenario, the effect of the climate would be more PLOS ONE Air temperature forecasting: Mutli-time scale PLOS ONE | https://doi.org/10.1371/journal.pone.0277079 November 3, 2022 12 / 31

obvious on the performance of the AI models. Accordingly, the models in this scenario will be trained using some recent records. 2.8 Statistical metrics Different statistical metrics assess the best model accuracy in daily and weekly temperature forecasting. Furthermore, it is vital to recognize the most efficient model with the least fore- casted error. Four statistical measures can be adopted to examine the forecasting accuracy of the suggested modeling approaches, such as root mean square error (RMSE), correlation coef- ficient (R), Thiels’ U-statistics (U), and mean absolute error (MAE). The mathematical expres- sions of these measures are presented below [57, 58]. RMSE ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 N X N i¼1 ðAi PiÞ2 s ð23Þ MAE ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 N X N i¼1 jAi Pij s ð24Þ U ¼ PN i¼1 ðAi PiÞ2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiPN i¼1 Ai2 q þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiPN i¼1 Pi2 q ð25Þ R ¼ PN i¼1 Ai A ð Þ Ai P ð Þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiPN i¼1 ðAi A Þ2 q ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiPN i¼1 ðPi P Þ2 q ð26Þ Where, Ai, and Pi are the actual and forecasting temperature for ith observation. While A−, and P− are the mean of actual and predicted value, and N is the total number of observations. The

stated statistical parameters have been used frequently in the literature for model comparison, and their estimation can be achieved directly from the observed and predicted values. Based on the results of these statistical measures, the model presenting the lowest value of forecasting error and the highest value of R (close to one) is selected as the best model for predicting the air temperature for short- and mid-term forecasting. 3. Result and discussion 3.1 First scenario This scenario investigates the capability of the AI models to predict the one-step ahead values for daily and weekly temperature. In this part of the work, we divide the data into two phases: training and testing. The classical method is used to separate the data into two steps; the first 70% of the recorded temperature is used for training, and the last 30% of the data is used for testing. In this scenario, the effect of the current temperature trends is not considered. In other words, the current records of temperatures are used in the testing phase. Thus, the models are tested and evaluated based on their ability to predict the current temperature values of the time series data. Furthermore, the input lags were determined by AFC and produced to the adopted models like RF, SVR, GBR, RT, and QRT. Different statistical parameters and compa- rable figures are used to assess the models’ performances. This part discusses the performance of the proposed models for predicting the daily and weekly temperature over the training and testing phases for different input lags. Tables 4 and 5

show the performance of the proposed models during the training phase for both daily and PLOS ONE Air temperature forecasting: Mutli-time scale PLOS ONE | https://doi.org/10.1371/journal.pone.0277079 November 3, 2022 13 / 31

weekly temperature prediction. For daily forecast, all models provide satisfactory results. Nev- ertheless, the RF provides the best performance, with RMSE ranging from 1.807 to 2.824, R ranging from 0.978 to 0.991, U ranging from 0.065 to 0.101, and MAE ranging from 1.389 to 2.157, followed by the QRT model with RMSE ranging from 2.0874 to 3.1830, R ranging from 0.9722 to 0.9882, U ranging from 0.0745 to 0.1134, and MAE ranging from 1.3492 to 2.24, and RT model with RMSE ranging from 2.3687 to 3.09, R ranging from 0.9738 to 0.9849, U rang- ing from 0.0849 to 0.1106, and MAE ranging from 1.7951 to 2.374. While GBR and SVR mod- els came last with RMSE ranging from 3.4233 to 3.6633 and from 3.6821 to 3.7887, R ranging from 0.963 to 0.9677 and from 0.9604 to 0.9626, U ranging from 0.1229 to 0.1328 and from 0.1321 to 0.1353, and MAE ranging from 2.6659 to 2.8746 and from 2.8369 to 2.8965, respec- tively. Furthermore, increasing the number of input lags increases the accuracy of the QRT, RT, and SVR models. Here the QRT model reaches the optimum accuracy when the input lag is nine (QRT − M9), the RT model with eight input lags (RT − M8), SVR model with ten input lags (SVR − M9) and the GBR model with five input lags (GBR − M5). On the contrary, the RF model only requires one input lag to reach its optimum performance RF–M1. For weekly temperature prediction, all models provide satisfactory predictions reaching the best performance with the QRT model with RMSE ranging from 2.5069 to 3.7157, R ranging

from 0.9589 to 0.9817, U ranging from 0.0932 to 0.11386, and MAE ranging from 1.6300 to 2.5880, followed by RF model with RMSE ranging from 2.799 to 3.647, R ranging from 0.960 to 0.977, U ranging from 0.105 to 0.136, and MAE ranging from 2.116 to 2.787. Moreover, increasing the input lags for weekly prediction improves the accuracy of QRT, RF, RT, and SVR models which allows them to reach their optimum accuracy (QRT − M16, RF − M16, Table 4. The performance of the proposed models for daily temperature prediction: Training phase. Models GBR QRT RT SVR RF RMSE R U MAE RMSE R U MAE RMSE R U MAE RMSE R U MAE RMSE R U MAE M1 3.663 0.963 0.132 2.829 3.183 0.972 0.113 2.24 3.09 0.974 0.111 2.374 3.789 0.96 0.135 2.897 1.807 0.991 0.065 1.389 M2 3.573 0.965 0.128 2.771 2.763 0.979 0.099 1.878 2.851 0.978 0.102 2.171 3.745 0.961 0.134 2.886 2.824 0.978 0.101 2.157 M3 3.553 0.965 0.128 2.759 2.621 0.981 0.094 1.768 2.75 0.979 0.099 2.123 3.734 0.961 0.134 2.873 2.756 0.979 0.099 2.117 M4 3.521 0.966 0.127 2.734 2.357 0.985 0.084 1.562 2.529 0.983 0.091 1.923 3.724 0.962 0.133 2.864 2.521 0.983 0.090 1.919 M5 3.423 0.968 0.123 2.666 2.281 0.986 0.081 1.51 2.486 0.983 0.089 1.896 3.716 0.962 0.133 2.861 2.488 0.983 0.089 1.895 M6 3.472 0.967 0.125 2.701 2.267 0.986 0.081 1.49 2.488 0.983 0.089 1.892 3.709 0.962 0.133 2.853 2.497 0.983 0.090 1.904 M7 3.488 0.967 0.126 2.72 2.109 0.988 0.075 1.37 2.4 0.984 0.086 1.817 3.698 0.962 0.133 2.851 2.366 0.985 0.085 1.792

M8 3.639 0.965 0.133 2.875 2.094 0.988 0.075 1.361 2.369 0.985 0.085 1.795 3.692 0.962 0.132 2.841 2.368 0.985 0.085 1.800 M9 3.503 0.966 0.126 2.734 2.087 0.988 0.074 1.349 2.387 0.985 0.086 1.812 3.682 0.963 0.132 2.837 2.369 0.985 0.085 1.803 https://doi.org/10.1371/journal.pone.0277079.t004 Table 5. The performance of the proposed models for weekly temperature prediction: Training phase. Models GBR QRT RT SVR RF RMSE R U MAE RMSE R U MAE RMSE R U MAE RMSE R U MAE RMSE R U MAE M10 3.898 0.954 0.146 3.02 3.716 0.959 0.139 2.588 3.663 0.96 0.137 2.792 4.648 0.935 0.174 3.607 3.647 0.960 0.136 2.787 M11 3.797 0.957 0.143 2.971 3.486 0.964 0.13 2.379 3.522 0.963 0.132 2.713 4.607 0.936 0.173 3.584 3.458 0.965 0.130 2.669 M12 3.837 0.956 0.144 2.992 3.256 0.969 0.121 2.215 3.578 0.962 0.135 2.72 4.598 0.936 0.173 3.572 3.455 0.965 0.130 2.668 M13 3.333 0.967 0.125 2.605 2.898 0.975 0.108 1.93 3.167 0.971 0.119 2.398 4.543 0.938 0.17 3.528 3.126 0.972 0.117 2.400 M14 3.125 0.971 0.117 2.428 2.756 0.978 0.103 1.827 3.036 0.973 0.114 2.332 4.473 0.94 0.167 3.459 3.075 0.973 0.115 2.341 M15 3.466 0.965 0.131 2.72 2.645 0.98 0.099 1.737 3.019 0.974 0.113 2.302 4.422 0.941 0.166 3.383 3.036 0.973 0.114 2.305 M16 3.392 0.967 0.128 2.651 2.507 0.982 0.093 1.63 2.816 0.977 0.105 2.142 4.308 0.944 0.162 3.303 2.799 0.977 0.105 2.116 M17 3.827 0.962 0.147 3.044 2.853 0.976 0.106 1.880 2.822 0.977 0.106 2.144 3.532 0.963 0.131 2.689 2.737 0.978 0.102 2.109

https://doi.org/10.1371/journal.pone.0277079.t005 PLOS ONE Air temperature forecasting: Mutli-time scale PLOS ONE | https://doi.org/10.1371/journal.pone.0277079 November 3, 2022 14 / 31

RT − M16, SVR − M16). It is observed that any increase in the number of lags beyond a value of seven imparts a negative effect on the model performance. At the same time, the GBR model requires five lags (same as the daily prediction) to reach its optimum accuracy (GBR − M14). Overall, QRT, RF, RT, and SVR models better predict daily temperature than weekly during the training phase. On the other hand, the GBR model better predicts the weekly temperature. Based on the training phase, all models perform very well in predicting the daily and weekly temperature. However, the assessment of the model performances based on the testing dataset is also crucial. For the training phase, the models are provided with complete data (input and targets), which may result in overfitting. Thus, excluding the models’ performances in the test- ing phase may provide users with misleading results. It is known that in the testing phase mod- els received only input features and thus the forecasting accuracy would be more reliable than in the training phase [46, 59]. During the testing phase, the performance of the proposed models was assessed firstly by comparing the performance with each other and secondly by comparing the performance with the ARIMA model as a benchmark model. Table 6 shows the performance of the proposed models during the testing phase for daily temperature prediction. For daily forecast, the RF model provides the best performance, with RMSE ranging from 1.776 to 3.765, R ranging from 0.960 to 0.991, U ranging from 0.063 to 0.133, and MAE ranging from 1.353 to 2.898 fol-

lowed by the SVR model with RMSE ranging from 3.5915 to 3.6599, R ranging from 0.9621 to 0.9635, U ranging from 0.1265 to 0.1288, and MAE ranging from 2.7451 to 2.7902. Moreover, the RF model requires only one input lag (RF − M1) to reach the best accuracy, while the SVR model requires five input lags. On the other hand, despite the RT and QRT models showing the best performance during the training phase, they came last during the testing phase as they have a tendency to overfit during the training phase. For further assessment, the ARIMA model was implemented for daily and weekly predic- tions using two different scenarios. The first one, the ARIMA used for the prediction of tem- perature using a raw data set. However, the second one, data preprocessing, is used to improve the capacity of ARIMA. At that stage, the differencing method is used to remove seasonality. It is possible to utilize that method to get rid of the temporal reliance, also known as the series dependence on time. The best prediction results are used as a benchmark to validate the AI models. Itis important to mention that the time series data became smoother after the applica- tion of differencing transformation technique (see Fig 4a and 4b). Considering the PACF Table 6. The performance of the proposed models for daily temperature prediction: Testing phase. Models GBR QRT RT SVR RF� RMSE R U MAE RMSE R U MAE RMSE R U MAE RMSE R U MAE RMSE R U MAE M1 3.697 0.961 0.131 2.84 4.015 0.954 0.141 3.07 3.978 0.955 0.14 3.055 3.66 0.962 0.129 2.79 1.776 0.991 0.063 1.353

M2 3.654 0.962 0.129 2.808 3.808 0.959 0.134 2.913 3.748 0.96 0.132 2.883 3.616 0.963 0.127 2.77 3.765 0.960 0.133 2.898 M3 3.639 0.963 0.129 2.807 3.735 0.96 0.132 2.866 3.725 0.961 0.132 2.872 3.606 0.963 0.127 2.758 3.740 0.960 0.132 2.884 M4 3.634 0.963 0.129 2.795 3.744 0.96 0.132 2.866 3.693 0.961 0.131 2.849 3.602 0.963 0.127 2.75 3.687 0.961 0.130 2.841 M5 3.623 0.963 0.128 2.794 3.701 0.961 0.13 2.839 3.67 0.962 0.13 2.831 3.592 0.964 0.127 2.745 3.678 0.962 0.130 2.836 M6 3.612 0.963 0.128 2.791 3.697 0.961 0.13 2.851 3.691 0.961 0.131 2.848 3.593 0.963 0.127 2.747 3.686 0.962 0.130 2.841 M7 3.617 0.963 0.128 2.791 3.714 0.961 0.131 2.842 3.68 0.962 0.13 2.855 3.593 0.963 0.127 2.755 3.696 0.961 0.131 2.850 M8 3.714 0.963 0.134 2.917 3.727 0.961 0.131 2.843 3.684 0.962 0.13 2.838 3.593 0.963 0.127 2.752 3.690 0.961 0.131 2.849 M9 3.633 0.963 0.129 2.813 3.691 0.961 0.13 2.835 3.692 0.961 0.131 2.857 3.595 0.963 0.127 2.758 3.720 0.961 0.132 2.872 � Symbol is the hyperparameter for the best model are number of trees = 4, and leaf node = 6. https://doi.org/10.1371/journal.pone.0277079.t006 PLOS ONE Air temperature forecasting: Mutli-time scale PLOS ONE | https://doi.org/10.1371/journal.pone.0277079 November 3, 2022 15 / 31

presented in Fig 4c for daily temperature prediction, three input lags (2,3, and 4 days) are con- sidered for the model development. As shown in Table 7, the ARIMA model performs on par with the QRT model and underperforms in comparison to the other models for daily scale. Furthermore, the ARIMA model requires only three input lags (ARIMA 3,1,2) to reach its best performance. On top of that, the application of the data transformation approach enhances the efficiency of the ARIMA prediction considerably. On the other hand, compared to the ARIMA model, the RF improves the prediction accuracy by 3% in terms of R and reduces the prediction error by 52.33%, 51.16%, and 52.1% in terms of RMSE, U, and MAE, respectively. For weekly temperature prediction, as presented in Table 8, the RF model demonstrates excellent performance in weekly temperature prediction by providing high prediction accu- racy with R ranging from 0.933 to 0.982 and fewer prediction errors with RMSE ranging from 2.478 to 4.665, U ranging from 0.091 to 0.173 and MAE ranging from 1.874 to 3.614 compared to the other models. Furthermore, the high performance of the RF model was achieved using only one input lag (RF − M10), while the other models require significantly higher input lags (seven lags) to reach their optimum performance (SVR − M16, RT − M16, GBR − M16, Fig 4. Input determination for ARIMA model, a) is the original data, b) after applying differencing method, c) PACF for daily scale. https://doi.org/10.1371/journal.pone.0277079.g004 PLOS ONE Air temperature forecasting: Mutli-time scale

PLOS ONE | https://doi.org/10.1371/journal.pone.0277079 November 3, 2022 16 / 31

QRT − M16) and increasing the lags beyond seven tends to reduce the models’ performance. Overall, the proposed models performed slightly better in predicting the daily temperatures than the weekly ones. On the other hand, the performance of the ARIMA model for weekly temperature predic- tion is presented in Table 9. Notably, the differencing method smoothens the time series data by removing a seasonal signal from a series (see Fig 5a and 5b). According to Fig 5c, based on PACF, three lags (1, 2, and 3 weeks) have been considered for the model development. As shown in Table 9, the performance of the ARIMA model is significantly lower than the RF model, and the latter has been able to increase the prediction accuracy by 5.36% in terms of R and reduce the prediction error by 48%, 47%, and 48% in terms of RMSE, U, and MAE. Table 7. The performance of the ARIMA model for daily temperature prediction: Testing phase. Without removing With removing Model RMSE R U MAE RMSE R U MAE ARIMA (2,1,1) 21.868 -0.053 0.596 17.418 3.727 0.961 0.13 2.828 ARIMA (2,2,1) 169.046 -0.034 0.903 142.136 205.266 0.115 0.949 156.224 ARIMA (2,1,2) 21.489 -0.053 0.594 17.034 3.738 0.961 0.129 2.825 ARIMA (3,1,1) 21.526 -0.053 0.594 17.073 3.727 0.961 0.129 2.826 ARIMA (3,1,2) 21.493 -0.053 0.594 17.038 3.726 0.961 0.129 2.825 ARIMA (4,1,2) 21.714 -0.053 0.595 17.261 3.786 0.961 0.13 2.859 ARIMA (4,1,1) 21.361 -0.053 0.593 16.91 3.786 0.961 0.13 2.859 https://doi.org/10.1371/journal.pone.0277079.t007

Table 8. The performance of the proposed models for weekly temperature prediction: Testing phase. Models GBR QRT RT SVR RF� RMSE R U MAE RMSE R U MAE RMSE R U MAE RMSE R U MAE RMSE R U MAE M10 4.693 0.932 0.172 3.62 4.922 0.925 0.179 3.795 4.796 0.929 0.175 3.721 4.54 0.936 0.167 3.465 2.478 0.982 0.091 1.874 M11 4.611 0.934 0.17 3.556 4.703 0.932 0.172 3.632 4.615 0.934 0.169 3.561 4.51 0.937 0.166 3.476 4.580 0.935 0.168 3.516 M12 4.575 0.935 0.168 3.475 4.73 0.931 0.174 3.628 4.654 0.933 0.172 3.588 4.5 0.938 0.166 3.463 4.665 0.933 0.173 3.614 M13 4.618 0.934 0.17 3.5 4.649 0.933 0.17 3.556 4.605 0.934 0.17 3.576 4.461 0.939 0.164 3.44 4.597 0.935 0.169 3.564 M14 4.557 0.936 0.167 3.513 4.653 0.933 0.171 3.564 4.585 0.935 0.17 3.525 4.416 0.94 0.162 3.404 4.562 0.936 0.169 3.525 M15 4.443 0.94 0.165 3.494 4.487 0.938 0.165 3.404 4.441 0.939 0.165 3.406 4.365 0.941 0.161 3.348 4.492 0.938 0.167 3.506 M16 4.423 0.941 0.165 3.472 4.423 0.94 0.163 3.382 4.247 0.945 0.157 3.291 4.286 0.944 0.158 3.283 4.358 0.941 0.161 3.349 M17 4.582 0.941 0.174 3.707 4.552 0.936 0.168 3.521 4.297 0.943 0.160 3.262 4.613 0.934 0.169 3.495 4.293 0.943 0.159 3.305 � Symbol is the hyperparameter for the best model are number of trees = 50, and leaf node = 5. https://doi.org/10.1371/journal.pone.0277079.t008 Table 9. The performance of the ARIMA model for weekly temperature prediction: Testing phase. Without removing With removing Model RMSE R U MAE RMSE R U MAE