The global warming causes changes to the Earth’s climate, or long-term weather patterns that vary from place to place. While we think of “Global warming” and “Climate change” as synonyms, scientists use the term climate change when describing the complex shifts affecting our planet’s weather and climate systems in different parts, because some areas actually get cooler in the short term, while the others become warmer.
Climate change encompasses not only rising average temperatures but also extreme weather events, shifting wildlife populations and habitats, rising seas and a range of other impacts. All of those changes are emerging as humans continue to add heat-trapping greenhouse gases to the atmosphere, changing the rhythms of climate that all living things have come to rely on. It has become clear that humans have caused most of the past century’s warming by releasing heat-trapping gases called “greenhouse gases”. Their levels are higher now than at any time in the last 800,000 years and, as a result, glaciers are melting, sea levels are rising and cloud forests are dying.
2. Global Temperature and the Greenhouse Effect
The warming that happens when certain gases in Earth’s atmosphere trap heat is considered as the greenhouse effect. These gases let in light but keep heat from escaping, just like the glass walls of a greenhouse, hence the name Greenhouse. Scientists have known about the greenhouse effect since 1824, when Joseph Fourier calculated that the Earth would be much colder if it had no atmosphere (   ). This natural green house effect is what keeps the Earth’s climate livable; and without it, the Earth’s surface would be an average of about 60˚F (33˚C) cooler.
2.1. Global Average Temperature
The concept of “global average temperature” is convenient for detecting and tracking changes in planet’s energy budget that is how much sunlight Earth absorbs minus how much it radiates to space as heat over time. The concept of an average temperature for the entire globe may sometimes seem odd, as the highest and lowest temperatures on Earth are about more than 55˚C or 100˚F apart. In the Northern and Southern Hemispheres, temperatures vary from night to day and between seasonal extremes, means that some parts of Earth are quite cold while other parts are downright hot.
In order to calculate a global average temperature, scientists begin with temperature measurements taken at various locations around the globe. Because the goal is to track changes in temperature, these measurements are converted from absolute temperature readings to “temperature anomalies”. These are the differences between the observed temperature readings and the long-term average temperature for each location and time. Multiple independent research groups across the world performed their own analysis of the surface temperature data, and they all showed a similar trend in upward direction  .
2.2. Trends in Northern and Southern Hemisphere Temperature
From increasing greenhouse gas concentrations, different parts of the world respond in different ways to warming. For example, high-latitude regions including far north or south of the equator become warm faster than the global average due to positive feedbacks from the retreat of ice and snow, an increased transfer of heat from the tropics to the poles in a warmer world also enhances warming.
2.3. Warmest Years on the Earth
According to the American Meteorological Society’s State of the Climate in 2017, the year brought an end to new record temperatures that were set each year from 2014 to 2016. Depending on the data set used, 2017 came in second or third warmest, after 2016 (warmest) and 2015 (second or third warmest)  . The near-record temperatures occurred in the absence of “El Niño” event, which is usually a factor in extreme global warmth. For much of 2017, “El Niño-Southern Oscillation (ENSO)” conditions were neutral, and October 2017 brought the start of “La Niña”, which typically drops global temperatures. Despite this, 2017 readings were 0.38˚C - 0.48˚C (or 0.68 - 0.86˚F) above the average of 1981-2010. Hence 2017 was the warmest non-El Niño year in the instrumental record (   ).
Figure 1 depicts the Global surface temperature in 2017 compared to the average temperature during 1981-2010. From this, it is clear that temperatures across most of the planet had been warmer than average during 1981-2010 (red colors). The high latitudes of the Northern Hemisphere were especially warm. Based on NOAA data  , the 2017 average global temperature across both the land and ocean surface areas was 0.84˚C (1.51˚F) above the 1901-2000 average of 13.9˚C (57.0˚F). This is making 2017 as the third-warmest year on record behind 2016 (warmest) and 2015 (second warmest). Furthermore, it was the warmest non-El-Niño year in the record  . It is also noted that since the start of the
Figure 1. Global surface temperature in 2017 compared to the average temperature during 1981-2010. Source: NOAA Climate.gov map, based on data from NOAA NCEI  .
twenty-first century, the annual global temperature record has been broken five times, The top 10 warmest years on record have all occurred since 1998, and the four warmest years on record have all occurred since 2014.
3. Literature Review
In this section, we will review some existing literature on different models/methods used to measure the climate change.
 used recent advances in time series econometrics to estimate the relation among emissions of carbon dioxide and methane, the concentration of these gases, and global surface temperature. These models were estimated and specified to answer two questions; whether the human activity affects global surface temperature and whether the global surface temperature affects the atmospheric concentration of carbon dioxide and methane. In this study, regression results provided direct evidence for a statistically meaningful relation between radioactive forcing and global surface temperature. A simple model based on these results indicated that greenhouse gases and anthropogenic sulfur emissions were largely responsible for the change in temperature over the last 130 years.
 used statistical models consisting of a trend plus serially correlated noise fitted to observed climate data, for example global surface temperature, the trend and noise representing systematic change and other variations, respectively. When such a model was fitted, the estimated character of the noise determined the precision of the estimated trend. In this study, the results of fitting such models to global temperature implied that there was uncertainty in the amount of temperature change over the past century of up to 0.2˚C and that the change was significantly different from zero.
To characterize observed global and hemispheric temperatures, previous studies have proposed different types of data-generating processes (see e.g.     ). The most common among them are random walk and trend-stationary, however, these approaches offering contrasting views regarding how the climate system works.
 presented an analysis of the time series properties of global and hemispheric temperatures using modern econometric techniques. Their results showed that the temperature series can be better described as trend-stationary processes with a one-time permanent shock. They suggested that the climate change has affected the mean of the processes but not their variability. During the last century, it has manifested in global and Northern Hemisphere temperatures, while a second stage is yet possible in the Southern Hemisphere. They argued that significant anthropogenic interference with the climate system has already occurred.
In  , the authors provided evidence of anthropogenic influence over the warming of the 20th century is presented and the debate regarding the time-series properties of global temperatures is addressed in depth. The 20th century global temperature simulations produced for the Intergovernmental Panel on Climate Change’s Fourth Assessment Report and a set of the radiative forcing series used to drive them are analyzed using modern econometric techniques. Results show that both temperatures and radiative forcing series share similar time-series properties and a common nonlinear secular movement.
4. Data and Statistical Methodology
We obtained the Combined Land-Surface Air and Sea-Surface Water Temperature Anomalies (Land-Ocean Temperature Index, LOTI) from Goddard Institute for Space Studies (GISS), NASA https://data.giss.nasa.gov/gistemp/. The data are available for Global mean, Northern Hemisphere mean and Southern Hemisphere means (monthly, quarterly and annual) since 1880 to present, updated through the most recent month  .
Functional Time Series (FTS) and Sliced Functional Time Series (SFTS)
 first introduced the functional time series (FTS) models. Using these models, the interest lies in forecasting a series of functional data observed over time. The functional curves are observed (with error) at time t = 1, … n, and we wish to forecast the functions for times t = n + 1, … n + h. Let [ft(xj)] denote the observed data, where j = 1, … p. We assume that there are underlying L1 continuous and smooth functions [st(x)] such that:
where [ei.j] are independent and identically distributed variables with zero mean and unit variance, and allows for heteroskedasticity.
The technique in  uses non parametric smoothing on each curve ft(x) separately to obtain estimates of the smooth functions [st(x)]. Panelized regression splines are used for smoothing, and then a functional principal component approach  is used to decompose the time series of functional data into a number of principal components and their scores. The functional time series (FTS) model can be written as follows:
where Ψk(x) is the kth principal component, the set of coefficients [ ] are the corresponding scores, et(x) denote independent and identically distributed random functions with zero mean, and K is the number of principal components to be used.
To plot a functional time series,  proposed three new graphical methods. They include the rainbow plot, the “Functional Bagplot” and the functional highest density region “(HDR) Boxplot”. Their approach has a side benefit of identification of outliers, which may not be obvious from the plot of the original data. These outliers are two types, either 1) magnitude outliers (i.e. the curves lie outside the range of the vast majority of the data), or 2) they may be the shape outliers (the curves that are within the range of the rest of the data but they have different shape from other curves). It is also possible that the curves may exhibit a combination of these two features. The presence of the outliers may have serious effect on the modeling and forecasting series.
To detect the outliers from a functional time series, the first step is to obtain the functional curves and the data are transformed into sliced functional time series (SFTS). For this, the entire data are sliced for each year as a function of 12 months. These curves are plotted in rainbow order with red for the earlier years and violet for the most recent year. The functional curves are then projected into a finite dimensional subspace. The subspace R2 is chosen for simplicity. Each of the functional data point in R2 are ordered by 1) data depth and 2) data density, based on halfspace Bagplot in  and HDR Boxplot in  . Those curves with lowest depth and/or lowest density are considered to be the outliers (see  for details).
The functional bagplot uses halfspace location depths described in  which is based on the bivariate bagplot of  , applied to the first two principal component scores. The depth region Rk is the set of all θ, with r(θ, z) ≥ k. Since the depth regions form a series of convex hulls, we have for k2 > k1. The Tukey bivariate depth median is defined as the value of θ which minimizes r(θ, Z) if there is such a unique θ, otherwise it is defined as the center of gravity of the deepest region.
Functional HDR boxplot
The functional HDR boxplot is based on the bivarate HDR Boxplot  , which is applied to the first two principal component scores. The bivariate HDR boxplot is constructed using a bivariate kernel density estimate f(z), which is defined as
where Zi represents a set of bivariate points; Khi(×) = K(×/hi)/hi; K is the kernel function; and hi is the bandwidth for the ith dimension. The bandwidths were selected using smoothed cross validation. Using the kernel density estimates, a HDR is defined as
where fα is such that ∫Rαf(z)dz = 1 − α; that is, it is the region with probability of coverage 1 − α, where all points within the region have a higher density estimate than any of the points outside the region, hence the name highest density region.
5. Results of Statistical Analysis
Figure 2 represents the average monthly global temperature since 1880 till 2018, as compared to the long-term average of twentieth century. Though warming has not been uniform across the planet, the upward trend in the globally averaged temperature shows that more areas are warming than cooling, specifically after 1980. According to the international State of the Climate in 2017 report, it was observed that since 1901, the planet’s surface has warmed by 0.7˚ - 0.9˚ Celsius (1.3˚ - 1.6˚ Fahrenheit) per century, but the rate of warming has nearly doubled since 1975 to 1.5˚ - 1.8˚ Celsius (2.7˚ - 3.2˚ Fahrenheit) per century.
Figure 3 depicts the mean monthly temperature anomalies of Northern and Southern Hemispheres from 1880 to 2018. If we compare the two series, we can observe that the northern series is more volatile with increasing trend than the southern series. The trend is evident from 1980 in northern and from 1960 in the southern hemisphere. The larger values of temperature anomalies in the Northern hemisphere may be due to the fact that it comprises of more land areas (represented by green color), whereas the southern hemisphere has more ocean/sea areas (represented by blue color). The ocean temperatures increase more slowly
Figure 2. The graph shows average monthly global temperatures since 1880 compared to the long-term average (1901-2000). The zero line represents the long-term average temperature for the whole planet during the twentieth century.
Figure 3. Mean monthly temperature anomalies of northern hemisphere (green color) and southern hemisphere (blue color) during 1880-2018. The zero line represents the long-term average temperature during 1901-2000.
than land temperatures because the oceans lose more heat by evaporation and they have a larger heat capacity.
Next, the data are transformed into sliced functional time series. The first step is to obtain the functional curves. For this, the entire data are sliced for each year as a function of 12 months, as plotted in Figure 4 and Figure 5. These curves are plotted in rainbow order with red for the earlier years and violet for the most recent year. We use R package rainbow  to construct these plots.
Figure 4 shows the global temperature anomalies (1880-2018) as sliced functional time series. The corresponding series for northern and southern hemispheres are plotted in Figure 5. The curves are plotted in rainbow order, with earlier years as red and most recent years as violet. It confirms that the average temperature is continuously rising in recent years. Some of the anomalies in Northern hemisphere series are as high as 1.0 - 1.5 (considered to be as outliers).
Figure 4. Sliced Functional Time series of Global Temperature Annomalies 1880-2018.
Figure 5. (a) Sliced Functional Time series of Northern Hemisphere 1880-2018; (b) Sliced Functional Time series of Southern Hemisphere 1880-2018.
Variability in Northern series is higher than the variability in southern series due to more land areas in Northern Hemisphere.
Outlier Detection in Temperature Series
Next, the functional curves are projected into a finite dimensional subspace, the subspace R2 is chosen for simplicity. Based on halfspace bagplot  and HDR boxplot  , each of the functional data point in R2 are ordered by data depth and data density. Those curves that have either lowest depth or lowest density are considered to be the outliers.
1) The Functional Bagplots
2) The Functional HDR Plots
3) The Functional Bivariate plots
Functional bagplots for Global, Northern and Southern Hemisphere series are plotted in Figures 6-8. Their respective functional HDR plots are shown in Figures 9-11; whereas, the functional bivariate plots based on the first two principle
Figure 6. Functional Bagplot for Global temperature anomalies (1880-2018). Black represents the modal curve, whereas the outliers are represented by different colors. Inner and outer regions are plotted with dark grey and light grey colors respectively.
Figure 7. Functional bagplot for average monthly temperature anomalies for Northern Hemisphere. The Median curve is denoted by black color, along with its confidence interval (blue dotted lines). The outliers are represented by red, green, yellow, blue and purple colors. Inner and outer regions are plotted with dark grey and light grey colors respectively.
Figure 8. Functional Bagplot for average monthly temperature anomalies for Southern Hemisphere. The median curve is denoted by black color, along with its confidence interval (blue dotted lines). Inner and outer regions are plotted with dark grey and light grey colors respectively.
Figure 9. Functional HDR plot for global temperature anomalies (1880-2018). The black line represents the modal curve, whereas the outliers are represented by red, yellow, green, blue and purple colors. Inner and outer regions are plotted with dark grey and light grey colors respectively.
Figure 10. Functional HDR plot for Northern Hemisphere temperature anomalies (1880-2018). The black line represents the modal curve, whereas the outliers are represented by different colors. Inner and outer regions are plotted with dark grey and light grey colors respectively.
Figure 11. Functional HDR plot for Southern Hemisphere temperature anomalies (1880-2018). The black line represents the modal curve, whereas the outliers are represented by different colors. Inner and outer regions are plotted with dark grey and light grey colors respectively.
components are constructed in Figures 12-14 respectively. The outliers depicted from these plots are shown in Table 1. In global series, 1916, 2015 and 2016 are confirmed outliers from all three methods. The corresponding outlier years in the Northern hemisphere are 2015, 2016 and 2017, whereas 1997, 1998 are consistently appeared to be the outliers in Southern series.
Application of FTS Model
Next, we apply the functional time series model of  and obtain the forecasts for next twenty years (2019-2038). The various components of FTS model (the mean function, the first two bases functions and corresponding time series coefficients) for global and the two hemisphere aeries are plotted in Figures 15-17 respectively. Figures 18-20 show the forecasts of average temperature in the three series; again the years are plotted in rainbow order with earlier year in red and most recent year in violet color.
The forecasts clearly show the warming in the three series. For global series, the forecasts values are relatively lower for the months of January and February, highest in March and then they are expected to be lower for April-July, slightly increase for August and October and relatively lower for the other months (Figure 18). The Northern Series shows the similar pattern with highest values in the month of March and lowest in July and relatively smaller in the other months (Figure 19).
The Southern hemisphere series forecasts depict a different pattern. The forecasts curves show increase in the average temperature in the next twenty years, with the maximum temperature in May and August and minimum in the months of November and December (Figure 20).
Forecast Comparison with the other Models
Finally the forecasting performance of Sliced Functional Time Series (SFTS) model will be measured by Mean Error (ME), Mean Absolute Error (MAE), Root Mean Square Error (RMSE) and Mean Absolute Percentage Error (MAPE). These measures are described below:
Table 1. Outliers in global, Northern and Southern Hemisphere temperature series (1880-2018) using different methods of outlier detection.
*The years in bold are magnitude outliers with high values and beyond the outer region.
Figure 12. Functional Bivariate plot for global average temperature series with first two principal components are plotted. The Red asterisk is the sample median, whereas the inner and outer regions are plotted with dark grey and light grey colors respectively.
Figure 13. Functional bivariate plot for Southern Hemisphere data with first two principal components are plotted. The Red asterisk is the sample median, whereas the inner and outer regions are plotted with dark grey and light grey colors respectively.
Figure 14. Functional bivariate plot for Southern Hemisphere data with first two principal components are plotted. The Red asterisk is the sample median, whereas the inner and outer regions are plotted with dark grey and light grey colors respectively.
Figure 15. Different components of FTS models applied to the global average temperature during 1880-2018, along with 10-year forecasts and 80% prediction intervals of the time series coefficients.
Figure 16. Different components of FTS models applied to the Northern Hemisphere average temperature during 1880-2018, along with 10-year forecasts and 80% prediction intervals of the time series coefficients.
Figure 17. Different components of FTS models applied to Southern Hemisphere average temperature during 1880-2018, along with 10-year forecasts and 80% prediction intervals of the time series coefficients.
Figure 18. 20-year Forecasts (2019-2038) of Global Average Temperature Anomalies using Sliced Functional Time Series Model. The forecast years are plotted in rainbow order with red color for the initial year (2019) and violet for the latest year (2038).
Figure 19. Twenty-year Forecasts (2019-2038) of Northern Hemisphere Temperature Anomalies using Sliced Functional Time Series Model. The forecast years are plotted in rainbow order with red color for the initial year (2019) and violet for the latest year (2038).
Figure 20. Twenty-year Forecasts (2019-2038) of Southern Hemisphere Temperature Anomalies using Sliced Functional Time Series Model. The forecast years are plotted in rainbow order with red color for the initial year (2019) and violet for the latest year (2038).
where Yi denotes the observed value and Fi denotes the corresponding forecast value. These measures of forecast accuracy are also computed for ARIMA/SARIMA models of Box and Jenkins  , exponential smoothing state space models  and random walk with drift models. For this, the global temperature series was divided into two sets: the training set 1880-1988 (109 years) and the test set 1989-2018 (30-years). The out-of-sample forecasts accuracy is measured and the results are summarized in Table 2.
Forecasts from these models are plotted in Figures 21-23 respectively. From these figures and Table 2, it is clear that the forecasts obtained by SFTS models have smaller values of error measures with relatively narrow prediction intervals.
As part of the Paris Agreement on climate change  , the international community committed in 2015 to limit rising global temperatures to well below 2˚C by the end of the 21st century and to pursue efforts to limit the temperature increase even further to 1.5˚C. However, these global temperature targets mask a lot of regional variation that occurs as the Earth warms. For example, land warms faster than oceans, high-latitude areas faster than the tropics, and inland areas faster than coastal regions.
In this paper, the global temperature data are analyzed through sliced functional time series (SFTS) model, a relatively new method of forecasting, and the
Figure 21. 10-year Forecasts (2019-2028) of global average temperature using ARIMA model, along with 80% prediction intervals.
Figure 22. 10-year Forecasts (2019-2028) of global average temperature using ETS model, along with 80% prediction intervals.
Figure 23. 10-year Forecasts (2019-2028) of global average temperature using Random walk with drift model, along with 80% prediction intervals.
monthly forecasts for the next twenty years (2019-2038) are obtained along with 80% prediction intervals. These forecasts are also compared with the forecasts
Table 2. Out of sample forecasting performance of different models.
obtained from Autoregressive Integrated Moving Average (ARIMA), exponential smoothing state space (ETS) and random walk with drift (RWD) models. It is found that the Sliced Functional Time Series models performed better than standard ARIMA, ETS and RWD models and the forecasts obtained from SFTS models are not only more accurate and reliable, but also they have narrow prediction intervals as compared to other models.
By 2038, the SFTS model projects that the average global surface temperature is expected to be 1.05 degree Celsius warmer than 1901-2000 average in the month of March, 0.95 degrees warmer in the months of January and February and about 0.85 degrees warmer in other months (see Figure 18). This similarity in temperatures regardless of total emissions is a short-term phenomenon: it reflects the tremendous inertia of Earth’s vast oceans. The high heat capacity of water means that ocean temperature doesn’t react instantly to the increased heat being trapped by greenhouse gases.
Given the size and tremendous heat capacity of the global oceans, it takes a massive amount of accumulated heat energy to raise Earth’s average yearly surface temperature, even a small amount. Behind the seemingly small increase in global average surface temperature over the past century is a significant increase in accumulated heat. That extra heat is driving regional and seasonal temperature extremes, reducing snow cover and sea ice, intensifying heavy rainfall, and changing habitat ranges for plants and animals by expanding some and shrinking others.