Nonlinear Cointegrating Regression of the Earth’s Surface Mean Temperature Anomalies on Total Radiative Forcing

Abstract

This study proposes a nonlinear cointegrating regression model based on the well-known energy balance climate model. Specifically, I investigate the nonlinear cointegrating regression of the mean of temperature anomaly distributions on total radiative forcing using estimated spatial distributions of temperature anomalies for the Globe, Northern Hemisphere, and Southern Hemisphere. Further, I provide two types of nonlinear response functions that map the total radiative forcing level to mean temperature anomalies. The proposed statistical model provides a climatological implication that spatially heterogenous warming effects play a significant role in identifying nonlinear climate sensitivity. Cointegration and specification tests are provided that support the existence of nonlinear effects of total radiative forcing.

1. Introduction

Over the past few decades, the observed global mean surface temperature has increased. With such an evident fact, global warming has received rapidly increasing attention during the past decades. Since all countries are involved in both the causes and consequences of this issue in a variety of complex ways, there have been worldwide debates on global warming among scientists and policymakers. To identify whether human activities are causing the recent rise in global mean temperature or whether their effects will have serious effects on the Earth, the detection and attribution of an anthropogenic influence on climate change has been studied extensively.

Climate sensitivity measures how much global warming will occur in response to a doubling of atmospheric concentration of carbon dioxide. Ultraviolet light from the Sun passes through greenhouse gases (GHG) such as carbon dioxide, water vapor, methane, nitrous oxide, and chlorofluorocarbons, and is absorbed by objects on the ground. Since GHG absorb the infrared radiation released by the objects and then reradiates it back to the surface of the Earth, global temperature is increasing as a result. We refer to this as the “Greenhouse Effect” (Hansen et al. 2011). It is well-known that GHG in the atmosphere have been consistently increasing due to human activity. Note that the atmospheric lifetime of carbon dioxide is currently estimated at 5–200 years (IPCC 2007). As a result of such accumulation, the stock of carbon dioxide in the atmosphere has increased by approximately 130 ppm over the last 270 years, from a range of between 275 and 285 ppm in the pre-industrial era to 410 ppm in 2018.

According to IPCC (2014), the global mean temperature on the surface of the Earth has increased by approximately 0.85 $°$C since 1880 and “most of the observed warming over the last 50 years is likely to have been due to the increase in GHG emissions.” Scientists have attempted to estimate the effect of GHG, according to which doubling carbon dioxide concentration in the atmosphere (a forcing of 4 W/m${}^2$) may increase the average global temperature by 1.5 to 4.5 $°$C. To stress the serious impact of increasing global temperature, Stern (2008) states that, “around 10,000–12,000 years ago, temperatures were approximately 5 $°$C lower than today, and ice sheets came down to latitudes just north of London and south of New York. As the ice melted and sea levels rose, England separated from the continent, rerouting much of the river flow. These magnitudes of temperature changes transform the planet”.

A time series analysis has been employed to test the anthropogenic global warming hypothesis. This hypothesis test has resulted in extensive controversy over the last two decades. The main argument relates to whether observing trends in temperature series and radiative forcing contain stochastic trends or deterministic trends with a structural break. During early years, Kaufmann et al. (2006a, 2006b, 2010, 2013) and Kaufmann and Stern (2002) had a breakthrough on the linear cointegration analysis between temperature series and radiative forcing variables by assuming they are integrated processes or difference stationary processes. Using Dynamic Ordinary Least Squares estimation, they concluded that the increase in global mean temperature can be associated with the change in radiative forcing variables. Such a linear cointegration analysis has also been investigated by Pretis (2020). He linked a two-component energy balance climate model of global mean temperature with a testable cointegrated Vector Autoregressive model.

Some econometricians cast doubt on their statistical rigor and challenged their empirical results (Gay-Garcia et al. 2009; Perron and Estrada 2012; Estrada et al. 2013a, 2013b; Estrada and Perron 2014). They first argued that temperatures and radiative forcing variables are described more effectively as trend stationary processes rather than difference stationary processes (or random walk with drift). By defining variables of interest as stationary processes fluctuating around a common breaking deterministic trend, they claimed that the conventional Least Squares (LS) method on the regression may negate a common feature as in cointegration analysis. Moreover, they argued that the residual-based ADF test (or formally nonparametric nonlinear co-trending test of Bierens 2000) may identify the existence of a long-run relationship.

Recently, Chang et al. (2020) analyzed the global warming issue under a novel time series framework. By using Global, Northern Hemisphere, and Southern Hemisphere temperature anomaly data from 1850 to 2012, they generated the distributions of temperature anomalies for each year. Instead of only analyzing mean temperature anomalies, they analyzed and tested the persistent features of distributions of temperature anomalies by regarding them as functional time series observations on the Hilbert space. More importantly, they distinguished unit-root nonstationarity from deterministic or explosive nonstationarity in their testing procedure Chang et al. (2016b, 2020). They concluded that the first few moments of temperature anomaly distribution indicate stochastic trends, rather than linear/exponential/quadratic trends or explosive roots.

In particular, they reasoned that the seemingly structural break in the global mean temperature anomaly trend, as argued by many authors, are more likely inherited from unit-root type persistency (stochastic trend), than from higher order persistency associated with deterministic trends. Based on their analysis, global temperature anomaly distribution and radiative forcing variables can share common stochastic trends and their linear combination can produce a stationary process. In this context, the nonlinear cointegration analysis, which allows for bidirectional causality and postulates a long-run equilibrium relationship between mean temperature anomalies and total radiative forcing and possibly with nonlinear moments of temperature anomaly distributions, seems the most reasonable approach.

More importantly, the linear regression model, which postulates the linear relationship between global mean temperature and total radiative forcing, fails to consider the climatological mechanism regarding how to change the global mean temperature. Since atmospheric carbon dioxide increases the global mean temperature by generating an imbalanced energy equation, a channel for the greenhouse effect must be considered. As shown in Section 2, the climate channel could be represented as temperature-dependent (or spatially heterogenous) net incoming absorbed radiation, implying that the relationship between total radiative forcing and the global mean temperature could be expressed as a temperature heterogeneous function. Evidently, ignoring the net incoming absorbed radiation term generates an endogeneity problem, invalidating the slope estimator of the linear model.

Throughout this study, I aim to develop a new statistical model to provide a more informative estimate for climate sensitivity. To do so, I propose a nonlinear cointegrating regression of the Earth’s surface mean temperature anomalies on total radiative forcing under technical backgrounds provided by Chang et al. (2020). Put differently, I analyze the cointegrating relationship between the time series of spatial distributions of temperature anomalies and the time series of total radiative forcing. In a climatological sense, the proposed statistical methodology provides two types of nonlinear climate sensitivity that map the total radiative forcing to mean temperature anomalies for the Globe, Northern Hemisphere, and Southern Hemisphere.

Specifically, I explicitly estimate the nonlinear effects by defining the nonlinear temperature term from the temperature anomaly distribution. I refer to the nonlinear response function as the temperature-dependent effect of total radiative forcing by assuming that net incoming absorbed radiation is hypothetically determined at some regional spaces that correspond to the temperature anomaly. I also define the misspecification error from the nonlinear cointegration model and identify the source of error in terms of temperature anomaly. Lastly, I conduct cointegration and specification tests that support the existence of nonlinear effects of total radiative forcing.

The statistical result of the study contributes an important insight to climate research. As is well-known, the popular term, “polar amplification” describes that the warming speed of the higher latitude area has been faster than that of the lower latitude area. The proposed nonlinear cointegration model provides a better understanding of such spatially heterogeneous warming effects by incorporating the higher-order moments of spatial distributions of temperature anomaly. Moreover, the notion that human forcing generates the probabilistic change of the temperature anomaly indicates that the true effect of human forcing on the mean temperature anomaly would be overestimated in the literature. In other words, ignoring the spatially heterogenous effects on the change in the global mean temperature anomaly would generate the upward bias when we attempt to identify the effect of human forcing.

The remainder of the paper is organized as follows. Section 2 provides the climatological background according to the global energy balance climate model. Section 3 presents the statistical model and methodology. Specifically, I employ the functional coefficient model for a nonlinear cointegrating regression model and apply it to the current climate model. In Section 4, the details of the data are discussed. The empirical results and interpretation of the results are presented in Section 5, while Section 6 concludes.

2. Climatological Analysis

2.1. Global Energy Balance Climate Model

North et al. (1983) developed the two-dimensional Energy Balance Climate Model (hereafter EBCM) and Brock et al. (2013) considered their model with human relating forcing activity. The EBCM model recently provided by Brock et al. (2013) is given by

$$\begin{array}{c}\hfill C\left(\widehat{r}\right)\frac{\partial s_{\widehat{r},t}}{\partial t}=\underset{\mathrm{Incoming}\mathrm{absorbed}\mathrm{radiation}}{\underbrace{QS\left(x\right)\alpha (x,s_{\widehat{r},t})}}+\underset{(-)\mathrm{Outgoing}\mathrm{radiation}}{\underbrace{h_{\widehat{r},t}-I_{\widehat{r},t}}}+\underset{\mathrm{Divergence}\mathrm{in}\mathrm{heat}\mathrm{flux}}{\underbrace{\nabla [D\left(x\right)\nabla s_{\widehat{r},t}]}}\end{array}$$

where $\widehat{r}=(\theta ,\varphi )$ is the point on the surface, $\theta \in [-\frac{\pi }{2},\frac{\pi }{2}]$ is the latitude and $\varphi \in [0,2\pi ]$ is the longitude, $C\left(\widehat{r}\right)$ is effective local heat capacity per unit area, and $s_{\widehat{r},t}$ is the Earth surface temperature at location $\widehat{r}$ and time t. By defining $x=sin\theta$, incoming solar radiation hitting the surface of the Earth is $QS\left(x\right)\alpha (x,s_{\widehat{r},t})$, where Q is the solar constant divided by 4, $S\left(x\right)$ is the mean annual meridional distribution of solar radiation, in which its integral from $-1$ to 1 must be unity, and $\alpha (x,s_{\widehat{r},t})$ is the absorption coefficient or co-albedo function, which is one minus the albedo of the Earth–atmosphere system. Outgoing radiation is $I_{\widehat{r},t}$, which is determined by surface temperature and clouds.

As in Brock et al. (2013) who include atmospheric CO${}_2$ concentration in the EBCM model, I included an external radiative forcing component, $h_{\widehat{r},t}$, at specific location $\widehat{r}$ and time t in the model. Following Hansen et al. (2017), however, I employed the sum of well-mixed GHG, ozone, surface albedo, tropospheric aerosols, and solar irradiance, in order to represent Total Radiative Forcing (hereafter, TRF). Clearly, the included TRF reflects human-related activities that influence the climate. The diffusion term for all different forms of heat transport, is $D\left(x\right)$. The heat transport (divergence in heat flux) term, $\nabla [D\left(x\right)\nabla s_{\widehat{r},t}]$, is due to incoming absorbed radiant heat not matched by net outgoing radiation.

As explained in the data section, global temperature data are expressed as anomalies from the zero-base period.1 Regarding data characteristics, I consider the discretized EBCM model at specific location $\widehat{r}$ and time t as given by

$$\begin{array}{c}\hfill C_{\widehat{r}}\left(s_{\widehat{r},t}-s_{\widehat{r},0}\right)=QS\left(x\right)\alpha (x,s_{\widehat{r},t})+h_{\widehat{r},t}-O_{\widehat{r}}{s}_{\widehat{r},t}+\nabla [D\left(x\right)\nabla s_{\widehat{r},t}]\end{array}$$

where $\widehat{r}=(x,\varphi )$ instead of $\widehat{r}=(\theta ,\varphi )$ for simplicity and $s_{\widehat{r},0}$ is the mean temperature at location $\widehat{r}$ during the zero-base period, $[\underline{t},\overline{t}]$. Note that the outgoing radiation term, $I_{\widehat{r},t}$, can be replaced with empirical formula, $O_{\widehat{r}}{s}_{\widehat{r},t}$ (Budyko 1969).2 Throughout this study, I consider the EBCM of the entire surface area of the Earth by integrating over $\widehat{r}$ (the Global Energy Balance Climate Model),

$$\begin{array}{cc}\hfill {\int }_{-1}^1{\int }_0^{2\pi }{C}_{\widehat{r}}\left(s_{\widehat{r},t}-s_{\widehat{r},0}\right)d\varphi dx=& \underset{\mathrm{Incoming}\mathrm{absorbed}\mathrm{radiation}}{\underbrace{{\int }_{-1}^1{\int }_0^{2\pi }QS\left(x\right)\alpha (x,s_{\widehat{r},t})d\varphi dx}}-\underset{\mathrm{Outgoing}\mathrm{radiation}}{\underbrace{{\int }_{-1}^1{\int }_0^{2\pi }{O}_{\widehat{r}}{s}_{\widehat{r},t}d\varphi dx}}\hfill \\ \hfill & +\underset{\mathrm{Total}\mathrm{radiative}\mathrm{forcing}}{\underbrace{{\int }_{-1}^1{\int }_0^{2\pi }{h}_{\widehat{r},t}d\varphi dx}}+\underset{\mathrm{Divergence}\mathrm{in}\mathrm{heat}\mathrm{flux}}{\underbrace{{\int }_{-1}^1{\int }_0^{2\pi }\nabla [D\left(x\right)\nabla s_{\widehat{r},t}]d\varphi dx}}.\hfill \end{array}$$

(1)

Note that the last term, ${\int }_{-1}^1{\int }_0^{2\pi }\nabla [D\left(x\right)\nabla s_{\widehat{r},t}]d\varphi dx$, is assumed to be zero because heat transports (typically from low to high latitudes) are negated across the entire surface area of the Earth. Moreover, dividing $4\pi$ on both sides of Equation (1) provides an interpretation that the global mean temperature anomaly is net incoming absorbed radiation across the globe plus TRF.

Although integrating divergence in heat flux across the globe is zero, there are spatial heterogeneities in net incoming absorbed radiation across the globe. Assuming a constant heat capacity (i.e., $C_{\widehat{r}}=C$ for all $\widehat{r}$)3 in Equation (1) provides

$$\begin{array}{cc}\hfill C\left(s_t-s_0\right)& =\frac{1}{4\pi }{\int }_{-1}^1{\int }_0^{2\pi }\left(QS\left(x\right)\alpha (x,s_{\widehat{r},t})-O_{\widehat{r}}{s}_{\widehat{r},t}+h_{\widehat{r},t}\right)d\varphi dx\hfill \\ \hfill & =\frac{1}{4\pi }{\int }_{-1}^1{\int }_0^{2\pi }\left(QS\left(x\right)\alpha (x,s_{\widehat{r},t})-O_{\widehat{r}}{s}_{\widehat{r},t}\right)d\varphi dx+h_t\hfill \\ \hfill & \equiv B_t\left(s_{\widehat{r},t}\right)+h_t\hfill \end{array}$$

(2)

where

$$B_t\left(s_{\widehat{r},t}\right)=\frac{1}{4\pi }{\int }_{-1}^1{\int }_0^{2\pi }\left(QS\left(x\right)\alpha (x,s_{\widehat{r},t})-O_{\widehat{r}}{s}_{\widehat{r},t}\right)d\varphi dx,h_t=\frac{1}{4\pi }{\int }_{-1}^1{\int }_0^{2\pi }{h}_{\widehat{r},t}d\varphi dx$$

and

$$\begin{array}{cc}\hfill s_t-s_0& =\frac{1}{4\pi }{\int }_{-1}^1{\int }_0^{2\pi }\left(s_{\widehat{r},t}-s_{\widehat{r},0}\right)d\varphi dx\hfill \end{array}$$

Note that $s_t$ is the global mean temperature at time t, $s_0$ is the global mean temperature during the zero-base period, $[\underline{t},\overline{t}]$, and therefore $s_t-s_0$ is the global mean temperature anomaly at time t. Further, note that integrating local temperatures across $x=sin\theta$ provides a way for temperatures in low latitudes to carry more weight in calculating the global mean temperature $s_t$, so that the increase in temperature at high latitudes due to divergence in heat flux from low to high latitudes, is roughly proportional to the temperature difference between high and low latitudes (Held and Suarez 1974 for one-dimensional EBCM model).

Equation (2) indicates that the spatially distributed TRF level along with different spatial characteristics over the Earth’s atmosphere, creates spatially heterogeneous warming effects. Consequently, the Global EBCM illustrates that spatially averaged net incoming absorbed radiation could be proportional to the global mean temperature anomaly we observed according to the data. More specifically, outgoing radiation to space is reduced by the TRF level at each location (greenhouse effect) and discrepancies between incoming radiation and outgoing radiation that is offset by exogenous forcing such as human and volcanic activity (without considering divergence in heat flux), generatetheirtorical global mean temperature anomaly data. In other words, the global mean temperature anomaly is determined by spatially heterogeneous net incoming absorbed radiation, which is a function of local temperatures, $s_{\widehat{r},t}$, as well as the TRF level in a complex way.

Notice that global climate change could be linked not only to ocean heat content, but also to the Top-of-Atmosphere (TOA) net energy imbalance. More specifically, the ocean absorbs more than 90 percent of excess energy inherited from human forcing (Loeb et al. 2012), and that clouds would enhance the anthropogenic influences by generating positive feedback (Fasullo and Trenberth 2012). Therefore, the Global EBCM should encompass the dynamics of both the change in TOA radiation and upper-ocean heating rate. Due to observational uncertainty and its complexity, however, I primarily focus on analyzing the Global EBCM (2) on the surface of the Earth.

2.2. Nonlinear Relationship between the Total Radiative Forcing Level and the Global Mean Temperature Anomaly

The precise estimation of climate sensitivity is of vital importance to climate scientists. The literature mainly considers three types of climate sensitivity. The first is Earth system sensitivity, which predicts resulting global warming associated with the doubled atmospheric CO${}_2$ concentration from the pre-industrial value. It attempts to reflect the long-run Earth system feedback such as the change in ice sheets and the Earth’s orbit for millennial timescales (Knutti et al. 2017). According to Vostok Ice-Core long-term temperature data for the last 420,000 years (Petit et al. 1999), more specifically, there has been a high correlation between atmospheric CO${}_2$ concentration and temperatures and sea level through the four glacial-interglacial cycles. The distribution of sunlight on Earth has experienced small changes in the Earth’s orbit over hundreds of thousands of years. Melting ice sheets due to more sunlight decreases the Earth’s surface albedo and increasingly warms the Earth. Consequently, the atmospheric CO${}_2$ concentration increased because the warming natural reservoirs such as the ocean release more carbon dioxide into the atmosphere (Hansen and Sato 2012).

The second one is equilibrium climate sensitivity, which measures the change in eventual (equilibrium) global surface mean temperature ($\Delta s_{eq}$) per unit of atmospheric CO${}_2$ concentration (h) or TRF (Hansen et al. 2013).4 The word “equilibrium” imposes a notion of the duration of time between two equilibrium climate states, which refer to the climate state in which the unit of TRF bumps the energy budget out of balance, and the climate state in which the Earth’s energy balance system is restored. Therefore, equilibrium climate sensitivity measures the difference in two of the Earth’s surface mean temperatures that compose the Earth’s energy balance body in a different time. Mathematically, equilibrium climate sensitivity is expressed as

$$\begin{array}{c}\hfill CS=\frac{\Delta s_{eq}}{h}.\end{array}$$

(3)

Note that the energy budget is determined by numerous feedbacks through which the temperature change affects the overall imbalance of the Earth’s climate system. For instance, carbon storage in the ocean can also be affected by ocean circulation. The decline in latitudinal overturning oceans (Ocean Conveyor Belt) as a result of global warming, slows down absorption of atmospheric CO${}_2$. Such nonlinear positive feedbacks enhance the increase in atmospheric CO${}_2$ concentration. However, the time scale of equilibrium climate sensitivity is moderate compared to Earth system sensitivity as it only considers the ocean heat uptake for decades.

The final one is the transient climate response. Similar to equilibrium climate sensitivity, it measures the magnitude of global warming per unit of atmospheric CO${}_2$ concentration or TRF. However, the difference is that it fails to consider the equilibrium climate states. Since the global surface temperature cannot respond to TRF or other climate forcing variables instantaneously due to thermal inertia of the climate system, and it is difficult to identify the equilibrium (steady) state from the observational dataset, the applied statistician would be more interested in measuring the transient climate response. Bindoff et al. (2013) provided an estimate of the change in temperature associated with a doubling of the concentration of carbon dioxide in the Earth’s atmosphere, which is likely to be in the range of 1 to 2.5 $°$C.

To estimate a constant climate sensitivity parameter, numerous authors estimate the slope of the following linear model: Kaufmann and Stern (2002) and Estrada et al. (2013b), inter alia.

$$\begin{array}{c}\hfill s_t-s_0=\mathrm{constant}+\mathrm{slope}\times h_t\end{array}$$

(4)

Note that the slope estimate indicates the change in the global surface mean temperature anomaly with respect to unit change in TRF and is expressed in units of $°$C/(W/m${}^2$), and therefore the estimation model (4) postulates that TRF effects are heterogeneous across the evaluated point, unlike the defined climate sensitivity in (3). More importantly, there is an endogeneity problem, inherited from the omitted term $B_t\left(s_{\widehat{r},t}\right)$ in estimating Equation (4), implying the invalidity of the slope estimator.

Specifically, the reduced form regression model for structural Equation (2) cannot be simply estimated, because two terms, $B_t\left(s_{\widehat{r},t}\right)$ and $h_t$, are intertwined physically. To produce physically sound estimates, there must be a channel for the TRF level, $h_t$, to affect the global mean temperature through net incoming absorbed radiation, $B_t\left(s_{\widehat{r},t}\right)$, (indirect effect) because TRF affects the global temperature by generating an imbalanced energy equation. In the absence of a net incoming absorbed radiation term, however, a constant TRF effect is evaluated at all temperature anomalies, assuming a linear dependence structure between the global mean temperature anomaly and TRF. In this light, considering the omitted term, $B_t\left(s_{\widehat{r},t}\right)$, in evaluating climate sensitivity is required.

Note that net incoming absorbed radiation cannot be entirely observable, and consequently the nonlinear climate sensitivity cannot be directly identified by unobservable information. Here, I consider that the TRF level affects the entire temperature anomaly distribution and that nonlinear temperature terms (or higher order moments of temperature anomaly distribution), which are defined as distributional temperature terms minus a linear temperature term, are regarded as information for net incoming absorbed radiation. Further, I assume that the way in which TRF affects net incoming absorbed radiation is temperature-dependent.

In the estimation procedure, I not only identify a constant (linear) TRF effect on the global mean temperature anomaly, but also estimate the remaining nonlinear TRF effect at each temperature anomaly. Climatologically, I may refer to the sum of a constant linear effect and a nonlinear effect of TRF at a temperature anomaly, $r_0$, the temperature-dependent TRF effect (or the TRF response function) by assuming that net incoming absorbed radiation is hypothetically and solely determined by a temperature anomaly, $r_0$. The positive gap between the TRF response function and a linear TRF effect at some temperature anomaly, therefore, implies that their relationship given the anomaly is underestimated. The reader is referred to Miller (2017) for the advanced climatological modeling.5

3. Statistical Model and Methodology

3.1. Reverse Functional Coefficient Model

In what follows, the statistical procedure is proposed to estimate the greenhouse effect. The reverse functional coefficient model is of particular interest of this study, as given by

$$\begin{array}{c}\hfill {\int }_{-\infty }^{\infty }{B}^{*}\left(r\right)f_t\left(r\right)dr=a_0+a_1{h}_t+{\epsilon }_t\end{array}$$

(5)

where $r_t$ is the random variable representing the raw temperature anomaly with density $f_t\left(r\right)$, and $h_t$ is the TRF level. The functional coefficient that is divided by some constant, $B^{*}\left(r\right)/a_1$, measures the TRF effect on the temperature anomaly, r. Therefore, the interpretation of Equation (5) is that the TRF level can affect the entire distribution of temperature anomaly in an arbitrary manner. The last term, ${\epsilon }_t$, is the mean-zero error. Note that I assume that there is no deterministic component in time series of spatial distributions of raw temperature anomaly, $r_t$ in Equation (5).

Here, it is assumed that $B^{*}\left(r\right)$ can be approximated by a Functional Fourier Flexible form as given by

$$B_{pq}^{*}\left(r\right)=\sum _{j=1}^p{b}_j^{*}{r}^j+\sum _{j=1}^q(b_{p+2j-1}^{*},b_{p+2j}^{*}){\varphi }_j\left(r\right),$$

where ${\varphi }_j\left(r\right)={(cos2\pi jr,sin2\pi jr)}^{\prime }$ as introduced by Gallant (1981) and extensively utilized by Park et al. (2010) and Chang et al. (2014). At first, I decompose the left-hand side into the linear part, $b_1^{*}r$, and the nonlinear part, $B_{pq}^{*}\left(r\right)$, where $B_{pq}^{*}\left(r\right)={\sum }_{j=2}^p{b}_j^{*}{r}^j+{\sum }_{j=1}^q(b_{p+2j-1}^{*},b_{p+2j}^{*}){\varphi }_j\left(r\right)$.

Specifically, the linear part contains a term according to which the functional coefficient is linear in temperature anomaly r and expressed as the mean temperature anomaly. On the other hand, the nonlinear part contains terms that refer to the functional coefficient as nonlinear in temperature anomaly r. As a result, I obtain

$$\begin{array}{c}\hfill b_1^{*}{\int }_{-\infty }^{\infty }r{f}_t\left(r\right)dr+{\int }_{-\infty }^{\infty }{B}_{pq}^{*}\left(r\right)f_t\left(r\right)dr=a_0+a_1{h}_t+{\epsilon }_t\end{array}$$

(6)

After rearranging Equation (6), I obtain

$$\begin{array}{c}\hfill {\int }_{-\infty }^{\infty }r{f}_t\left(r\right)dr=b_0+b_1{h}_t+{\int }_{-\infty }^{\infty }B\left(r\right)f_t\left(r\right)dr+u_t\end{array}$$

(7)

where $B\left(r\right)=-B_{pq}^{*}\left(r\right)/b_1^{*},b_0=a_0/b_1^{*},b_1=a_1/b_1^{*},$ and $u_t={\epsilon }_t/b_1^{*}$. Note that Equation (7) specifies the statistical relationship between the global mean temperature anomaly (the left-hand side), the TRF level (second term of the right-hand side), and the sum of probability-weighted net incoming absorbed radiation (third term of the right-hand side). Specifically, it is implicitly assumed that net incoming absorbed radiation, $B_t\left(s_{\widehat{r},t}\right)$, in the Global EBCM (2) is decomposed into linearly additive functions at each local temperature anomaly, $b\left(s\right)$, with associated probability weight, $f_t\left(s\right)$. Mathematically, it is expressed as

$$\begin{array}{c}\hfill B_t\left(s_{\widehat{r},t}\right)={\int }_{s}b\left(s\right)f_t\left(s\right)ds\equiv \mathbb{E}\left[b\left(s_t\right)\right]\end{array}$$

(8)

which evaluates the effect of net incoming absorbed radiation of the regions that are represented by the temperature anomaly s at time t. Therefore, the third term, $\int B\left(r\right)f_t\left(r\right)dr$, which is also expressed as $\mathbb{E}\left[B\left(r_t\right)\right]$ in Equation (7), provides a channel for net incoming absorbed radiation, $B_t\left(s_{\widehat{r},t}\right)$, in the Global EBCM (2). This term affects the global mean temperature anomaly through changed temperature probabilities. Consequently, the reduced form model (7) corresponds to the Global EBCM in Equation (2).6

Statistically, Equations (5) and (7) can be analyzed as a cointegrating relationship given by

$$\begin{array}{c}\hfill h_t-{\beta }_{pq}^{\prime }{x}_{pqt}\sim I\left(0\right)\end{array}$$

(9)

where the multivariate time series vector, $x_{pqt}={(\int r{f}_t\left(r\right)dr,\int r^2{f}_t\left(r\right)dr,\cdots ,\int {\varphi }_j\left(r\right)f_t\left(r\right)dr)}^{\prime }$, with cointegrating vector, ${\beta }_{pq}={(1,c_1,c_2,\cdots ,c_{p+2q})}^{\prime }$. Note that the mean-level regression implemented in climate research assumes that the TRF level only affects the aggregated temperature anomaly distribution towards the mean process. In this study, however, I consider the cointegrating relationship between the temperature anomaly distribution and the TRF level. It turns out that the relationship between the two is the cointegrating relationship between the global mean temperature anomaly, the TRF level, and nonlinear temperature terms.

Given the fact that GHG are well-mixed in the atmosphere, furthermore, the nonlinear cointegration analysis in this study becomes the distributional cointegration analysis between the temperature anomaly distribution and the TRF distribution (Chang et al. 2020). In other words, since TRF would be distributed homogeneously across the Earth’s atmosphere, it would be uninformative to consider the higher-order moments of the TRF distribution for the distributional cointegration analysis.

In particular, two types of estimators are provided from a nonlinear cointegration model. The first estimator, $D^1$, provides the slope estimator that considers the TRF effect on the global mean temperature anomaly through nonlinear temperature terms. On the other hand, the second estimator, $D^2\left(r\right)$, provides the temperature-dependent slope estimator, which is illustrated by the TRF response function. Statistically, the derivative $\partial (\int r{f}_t\left(r\right)dr)/\partial h_t$—which is equal to $b_1$ in the linear case—estimates climate sensitivity.

I first suppose that the spatial temperature anomaly distributions $f_t\left(r\right)$ across the globe for all times t are observed over the domain, $D_r=[{\lambda }_1,{\lambda }_2]$, where ${\lambda }_1$ and ${\lambda }_2$ are the lowest and highest temperature anomalies observed for the entire sample period, respectively. Moreover, it is assumed that the spatial temperature anomaly distributions themselves are regarded as random variables under predetermined temperature anomalies rs. Under fixed temperature anomalies rs, the total derivatives of $f_t(r,h_t)$ with respect to $h_t$ become the partial derivatives of $f_t(r,h_t)$ with respect to $h_t$. Thus, we have

$$\begin{array}{c}\hfill D^1\equiv \frac{\partial {\int }_{D_r}r{f}_t\left(r\right)dr}{\partial h_t}={\int }_{D_r}r\frac{\partial f_t\left(r\right)}{\partial h_t}dr=b_1+{\int }_{D_r}B\left(r\right)\frac{\partial f_t\left(r\right)}{\partial h_t}dr=b_1+\frac{\partial {\int }_{D_r}B\left(r\right)f_t\left(r\right)dr}{\partial h_t}\end{array}$$

(10)

Note that total differentiating Equation (7) and replacing it with partial derivatives provide Equation (10). As stated above, the last equality of Equation (10) holds under fixed temperature anomalies (i.e., $\int \frac{\partial B\left(r\right)}{\partial h_t}{f}_t\left(r\right)dr=0$), which means that the temperature anomalies rs are invariant with respect to $h_t$. The probabilities at fixed temperature anomalies, however, are affected by $h_t$.

More importantly, model (7) is reduced to the standard linear model in climate literature when $B\left(r\right)=0$ (Estrada et al. 2013a, 2013b). If $B\left(r\right)\ne 0$ and the term $\int B\left(r\right)f_t\left(r\right)dr$ is omitted, however, the estimator of climate sensitivity in the regression Equation (7) would be biased, and the conventional test statistics on estimator ${\widehat{b}}_1$ would be also invalid. By including a second term in Equation (10), temperature-dependent climate sensitivity can be identified and therefore the heterogeneous TRF effects on the global mean temperature anomaly can be identified at each anomaly. When net incoming absorbed radiation is hypothetically determined at $r_0\in [{\lambda }_1,{\lambda }_2]$, more specifically, the TRF effect on the global mean temperature anomaly is estimated by

$$\begin{array}{c}\hfill D^2\left(r_0\right)=b_1+B\left(r_0\right)d_T\left(r_0\right)\end{array}$$

(11)

where

$$d_T\left(r_0\right)=\frac{\partial f_t\left(r\right)}{\partial h_t}{|}_{r=r_0}$$

Note that $d_T\left(r\right)$ is estimated by the linear regression of the probability series at each temperature anomaly on TRF with a constant using the least squares method. I refer to the second terms in Equations (10) and (11) as the nonlinear effect of TRF on the global mean temperature anomaly. Subsequently, I define the misspecification error of the linear model, which represents error from ignoring the nonlinear effects, as

$$\mathrm{Error}\equiv b_1^{Res}-D^1=b_1^{Res}-\left(b_1^{Unres}+\frac{\partial {\int }_{D_r}B\left(r\right)f_t\left(r\right)dr}{\partial h_t}\right)$$

where $b_1^{Res}$ and $b_1^{Unres}$ are the slope estimators calculated by the linear and nonlinear regressions of the global mean temperature anomaly on TRF, respectively.

Further, note that the sum of nonlinear effects of estimator $D^2\left(r\right)$ across temperature anomalies rs is equal to the nonlinear effect of estimator $D^1$. In this sense, error calculation from estimator $D^2\left(r\right)$ enables us to identify the contribution of misspecification error at each temperature anomaly, $r_0\in [{\lambda }_1,{\lambda }_2]$, which is represented as

$$-B\left(r_0\right)\frac{\partial f_t\left(r\right)}{\partial h_t}{|}_{r=r_0}$$

Notably, nonlinear cointegration models (5) and (9) yield the implication that the linear model imposes an implausible restriction on the relationship between the temperature anomaly distribution and the TRF level. To illustrate this, consider the linear model as given by

$$\begin{array}{c}\hfill {\int }_{-\infty }^{\infty }{b}_1^{*}r{f}_t\left(r\right)dr=a_0+a_1{h}_t+{\epsilon }_t\end{array}$$

(12)

Note that only the linear part $b_1^{*}r$ of the functional coefficient $B_{pq}^{*}\left(r\right)$ is preserved in Equation (12). Then, the climate sensitivity of the linear model can be expressed as the temperature-weighted linear sum of the derivatives of the global temperature anomaly distribution with respect to the TRF level. Mathematically, climate sensitivity

$$b_1=a_1/b_1^{*}={\int }_{-\infty }^{\infty }r\frac{\partial f_t\left(r\right)}{\partial h_t}dr$$

assumes that the derivative at the high (low) temperature should be weighted more (less) than the low (high) temperature.

As this study provides a structural modeling for temperature anomaly on a yearly global scale, it is worth comparing the nonlinear cointegration model with the temperature model provided in the climate economics literature. In particular, Campbell and Diebold (2005) provide a daily average temperature model of the U.S. cities. Their econometric model consists of the conditional mean and variance components. The conditional mean part captures the deterministic trend using time trends, seasonality using a fourier series, and high-frequency cycle using autoregressive lags. In the meantime, the conditional variance part captures the seasonal volatility dynamics using a fourier series and the cyclical volatility dynamics using a generalized autoregressive conditional heteroscedasticity (GARCH) process.

Similarly, Morana and Sbrana (2019) provides a monthly average temperature model in hemispheric and regional scales. However, the econometric model of Morana and Sbrana (2019) differs from that of Campbell and Diebold (2005) in several aspects. First, the former captures the stochastic trend in the temperature dynamics using a radiative forcing variable. Second, the former does not allow the seasonal volatility dynamics, but impose the conditional covariance structure between the average temperatures in different regions. Third, the former allows a structural break in the conditional mean equation, while the latter does not.

Notice that the well-known temperature regression models such as the linear cointegration model of Kaufmann and Stern (2002) and the cotrending model of Estrada et al. (2013b) are different from the reduced-form regression models of Campbell and Diebold (2005) and Morana and Sbrana (2019). Specifically, the global temperature models employed by Kaufmann and Stern (2002) and Estrada et al. (2013b) assume that the radiative forcing is cointegrated, or shares a structural break with the mean temperature, and the conditional variance of mean temperature is time-invariant. More importantly, their global temperature models do not attempt to explain the short-run and medium-run variations in the temperature dynamics using autoregressive lags and a fourier series, explaining only the long-run component using the radiative forcing. In this light, the nonlinear cointegration model of this study would be regarded as a long-run version of Morana and Sbrana (2019) with extension for the spatial heterogeneity. Put differently, this study incorporates the spatial heterogeneity of the effects of radiative forcing in the existing long-run temperature model.

3.2. Estimation and Inference

To estimate the derivative of interest $D^1$ and $D^2\left(r\right)$ from Equations (10) and (11), $b_1$ and $B\left(r\right)$ must be estimated. Prior to displaying the estimation procedure, the time-varying domain of the temperature anomaly distribution is discussed. The Functional Fourier Flexible form has been exploited to generate smoothed nonlinearity as a basis function. The polynomial and trigonometric basis functions turn out to be useful in describing the nonlinear relationship between temperature and the economic variable (Chang et al. 2016a). It is worth noting that the literature restricts the support to unit interval, in order to ensure the linear independence between basis functions on the domain. Indeed, the $n$th-order trigonometric function is periodic with $\frac{1}{n}$ period, making the different order trigonometric functions linearly independent. Likewise, the different order polynomial functions are linearly independent.

Nevertheless, when we integrate the polynomial basis functions with the temperature density on the unit interval, the multicollinearity problem would be generated in the sense that the polynomial function on the unit interval could be interpreted as a weight function. Put differently, since the low-order polynomial functions on the unit interval, which are integrated with identical time series of temperature density, cannot generate a notable difference, the generated regressors would have the high collinearity. Further, note that Chang et al. (2020) generated the global temperature anomaly distribution on the bounded common support, and that common support is calculated using 95 percent of the total probability mass over the entire time span. By definition, however, the nonstationarity of the distributional time series is not compatible with bounded common support, and that restricting the support of the entire time series to the common domain could also generate the multicollinearity issue, implying that the nonstationary domain can be the source of variability of the time series.

Given that the global mean temperature anomaly has a stochastically increasing trend, more specifically, 2.5 percent of the probability mass at the left end would eliminate the heterogeneity of the beginning of the distributional time series, and 2.5 percent of the probability mass at the right end would eliminate the heterogeneity of the end of the distributional time series. Indeed, the correlation coefficient between ${\int }_{s}s{f}_t\left(s\right)ds$ and ${\int }_s{s}^2{f}_t\left(s\right)ds$ is calculated as a value, 0.994, where $s\in [0,1]$. Consequently, the regression analysis with such a high collinearity could be misdirected.

To overcome this, the variability inherited from the polynomial basis function on the time-varying support must be preserved. For this, 99 percent of the total probability mass of the temperature anomaly distribution at each time t is calculated, making the support of spatial distribution of temperature anomaly vary over time. The time-varying support of the raw temperature anomaly, $r_t$, is denoted as $D_t^r=[{\lambda }_t^1,{\lambda }_t^2]$. Moreover, the unit interval is only applied to the trigonometric basis function by normalizing the raw temperature anomaly (i.e., $s=\frac{r_t-{\lambda }_t^1}{{\lambda }_t^2-{\lambda }_t^1}\in [0,1]$). In the meantime, the support of the raw temperature anomaly is directly applied to the polynomial basis function, implying that sample moments of the temperature anomaly distribution are included in the regression equation.

Subsequently, the mathematical procedure to estimate the derivative of interest is provided,

$$D^1\equiv \frac{\partial {\int }_{D_t^r}{r}_t{f}_t\left(r\right)dr}{\partial h_t}=\frac{\partial {\int }_{D_s}\{{\lambda }_t^1+({\lambda }_t^2-{\lambda }_t^1)s\}{f}_t\left(s\right)ds}{\partial h_t}\mathrm{where}{r}_t={\lambda }_t^1+({\lambda }_t^2-{\lambda }_t^1)s\in D_t^r\equiv [{\lambda }_t^1,{\lambda }_t^2]D_s.$$

Note that we have

$$\begin{array}{cc}\hfill {\int }_{D_r}B\left(r\right)f_t\left(r\right)dr& ={\int }_{D_r}\left(B^p\left(r\right)+B^t\left(r\right)\right)f_t\left(r\right)dr\hfill \\ \hfill & ={\int }_{D_r}{B}^p\left(r\right)f_t\left(r\right)dr+{\int }_{D_s}{B}^t\left(s\right)f_t\left(s\right)ds\hfill \end{array}$$

due to the normalization-invariance property. Estimator ${\widehat{D}}^1$ can be estimated by

$$\begin{array}{c}\hfill {\int }_{D_t^r}{r}_t{f}_t\left(r\right)dr=b_0+b_1{h}_t+\sum _{i=2}^p{b}_i{\int }_{D_t^r}{r}_t^i{f}_t\left(r\right)dr+\sum _{j=1}^q(b_{p+2j-1}^s,b_{p+2j}^s){\int }_{D_s}{\varphi }_j\left(s\right)f_t\left(s\right)ds+(u_t-e_{pqt})\end{array}$$

(13)

where s denotes the normalized temperature, which is defined as $\frac{r_t-{\lambda }_t^1}{{\lambda }_t^2-{\lambda }_t^1}$ and bounded in unit interval (i.e., $D_s=[0,1]$). Define

$$\begin{array}{cc}\hfill {\beta }_{p0}& ={(b_2,b_3,\cdots ,b_p)}^{\prime },{\beta }_{0q}^s={(b_{p+1}^s,b_{p+2}^s,\cdots ,b_{p+2q-1}^s,b_{p+2q}^s)}^{\prime },\hfill \\ \hfill {\pi }_{p0}\left(r\right)& ={\left(r^2,r^3,\cdots ,r^p\right)}^{\prime },{\pi }_{0q}\left(s\right)={\left(cos\left(2\pi s\right),sin\left(2\pi s\right),\cdots ,cos\left(2q\pi s\right),sin\left(2q\pi s\right)\right)}^{\prime }\hfill \end{array}$$

Then, functional coefficients $B^p\left(r\right)$ and $B^t\left(s\right)$ are estimated by

$${\widehat{B}}^p\left(r\right)={\widehat{\beta }}_{p0}^{{}^{\prime }}{\pi }_{p0}\left(r\right),{\widehat{B}}^t\left(s\right)={\widehat{\beta }}_{0q}^{s^{\prime }}{\pi }_{0q}\left(s\right)$$

Subsequently, we calculate the TRF response function as

$$\widehat{B}\left(r\right)={\widehat{B}}^p\left(r\right)+{\widehat{B}}^t({\lambda }_1+({\lambda }_2-{\lambda }_1)s)$$

Finally, we estimate

$$\begin{array}{c}\hfill \frac{\partial {\int }_{D_r}\widehat{B}\left(r\right)f_t\left(r\right)dr}{\partial h_t}\end{array}$$

(14)

by implementing a linear regression of ${\int }_{D_t^r}\widehat{B}\left(r\right)f_t\left(r\right)dr$ on $h_t$ with a constant.

It is widely known that the LS method on cointegration analysis provides a super-consistent estimator. However, LS estimators may be inefficient and asymptotically biased. Moreover, the hypothesis testing based on LS estimator is invalid, due to the presence of nuisance parameters. Throughout this study, the canonical cointegrating regression developed by Park (1992) is employed.

In a matrix form, Equation (13) is re-written as

$$\begin{array}{c}\hfill y_t={\alpha }_{pq}^{\prime }{z}_{pqt}+u_{pqt}\end{array}$$

(15)

where

$$\begin{array}{cc}\hfill y_t& ={\int }_{D_t^r}{r}_t{f}_t\left(r\right)dr\hfill \\ \hfill z_{pqt}& ={\left(h_t,{\int }_{D_t^r}{r}_t^2{f}_t\left(r\right)dr,\cdots ,{\int }_{D_s}cos\left(2q\pi s\right)f_t\left(s\right)ds,{\int }_{D_s}sin\left(2q\pi s\right)f_t\left(s\right)ds\right)}^{\prime }\hfill \\ \hfill u_{pqt}& =e_{pqt}+u_t\hfill \end{array}$$

For the convenience of theoretical development, I assume both $y_t$ and $z_{pqt}$ are mean zero processes by taking demeaning or detrending procedure from Equation (15).

Assume $\left\{w_{pqt}\right\}$ defined by

$$w_{pqt}={(u_{pqt},\Delta z_{pqt})}^{\prime }$$

satisfies the Invariance Principle (IP). Define ${\varphi }_{pq}\left(k\right)=\mathbb{E}{w}_{pqt}{w}_{pq,t-k}^{\prime }$ as the autocovariance function of $\left(w_{pqt}\right)$. The long-run variance matrix $\mathsf{\Omega }$ of $\left\{w_{pqt}\right\}$ is then given by

$$\mathsf{\Omega }=\sum _{k=-\infty }^{\infty }\varphi \left(k\right)$$

Henceforth, the subscript $pq$ is suppressed at the variance matrix for notational convenience. The contemporaneous variance $\mathsf{\Sigma }=\varphi \left(0\right)$ and the one-sided long-run variance $\mathsf{\Lambda }={\sum }_{k=1}^{\infty }\varphi \left(k\right)$ is also defined so that $\mathsf{\Omega }=\mathsf{\Sigma }+\mathsf{\Lambda }+{\mathsf{\Lambda }}^{\prime }$. Further, define

$$\mathsf{\Gamma }=\mathsf{\Sigma }+\mathsf{\Lambda }$$

which is $(p+2q+1)\times (p+2q+1)$ matrix. Consider the partitioned submatrices of $\mathsf{\Omega }$ and $\mathsf{\Gamma }$ by

$$\begin{array}{c}\hfill \mathsf{\Omega }=\left(\begin{array}{cc}{\omega }_{11}& {\omega }_{12}\\ {\omega }_{21}& {\mathsf{\Omega }}_{22}\end{array}\right),\mathsf{\Gamma }=\left(\begin{array}{cc}{\gamma }_{11}& {\gamma }_{12}\\ {\gamma }_{21}& {\mathsf{\Gamma }}_{22}\end{array}\right)=\left[{\gamma }_1{\mathsf{\Gamma }}_2\right]\end{array}$$

by assuming ${\mathsf{\Omega }}_{22}>0$. The Canonical Cointegrating Regression (CCR) method proposed by Park (1992) is based on the regression

$$\begin{array}{c}\hfill y_t^{*}={\alpha }_{pq}^{\prime }{z}_{pqt}^{*}+u_{pqt}^{*}\end{array}$$

(16)

Here, $\left\{y_t^{*}\right\}$ and $\left\{z_{pqt}^{*}\right\}$ are stationary transformations of $\left\{y_t\right\}$ and $\left\{z_{pqt}\right\}$, which are given by

$$\begin{array}{cc}\hfill y_t^{*}& =y_t-{\alpha }_{pq}^{\prime }{\mathsf{\Gamma }}_2^{\prime }{\mathsf{\Sigma }}^{-1}{w}_{pqt}-w_{12}{\mathsf{\Omega }}_{22}^{-1}\Delta z_{pqt}\hfill \\ \hfill z_{pqt}^{*}& =z_{pqt}-{\mathsf{\Gamma }}_2^{\prime }{\mathsf{\Sigma }}^{-1}{w}_{pqt}\hfill \end{array}$$

(17)

The CCR estimator ${\widehat{\alpha }}_{pq}^{CCR}$ of ${\alpha }_{pq}$ is the LS estimator of ${\alpha }_{pq}$ in the transformed regression (16). Note that the long-run relationship between $y_t$ and $z_{pqt}$ in Equation (15) remains in the relationship between $y_t^{*}$ and $z_{pqt}^{*}$ in Equation (16) because the stationary transformation does not affect their long-run relationship. In this sense, the CCR estimate ${\widehat{\gamma }}^{CCR}$, must be the same for the LS estimate, ${\widehat{\gamma }}^{LS}$.

For the feasible CCR estimation, Park (1992) and Park et al. (2010) suggest estimating the long-run variance in Equation (17) with LS estimator ${\widehat{\alpha }}_{pq}^{LS}$ from Equation (15). Note that LS standard errors are not efficient under the presence of endogeneity of error term. Meanwhile, the CCR estimator could be problematic if the long-run variance is not consistently estimated due to nonstationarity in errors, inherited from misspecified orders p and q. In this case, we cannot guarantee that the CCR estimator is more effective than the LS estimator, implying that both CCR and LS methods may fail to provide correct standard errors under misspecified orders p and q. Moreover, there has been no clear consensus that the CCR estimator is more effective than the LS estimator in the nonlinear model. Note that the CCR estimator’s role as per Park et al. (2010) is to fix only the linear part of the regression model. In light of this, finding an exact order of p and q is a critical issue in this study.

Practically, it is assumed that the exact orders of p and q are known for the consistent long-run variance estimation with the Akaike Information Criterion (AIC) and the Bayesian information criterion (BIC), so that nonlinearity or nonstationarity in errors is effectively removed. Along the lines of Park et al. (2010), the consistent kernel estimator is considered for the long-run variance $\mathsf{\Omega }$ of ${\widehat{w}}_{pqt}={({\widehat{u}}_{pqt},\Delta z_{pqt})}^{\prime }$ that is given by

$${\widehat{\mathsf{\Omega }}}_n=\frac{1}{n}\sum _{\left|i\right|<l_n}\tau \left(\frac{i}{l_n}\right)\sum _t{\widehat{w}}_{pqt}{\widehat{w}}_{pq,t-i}$$

with lag window $\tau$ and lag truncation number $l_n$. The LS residual from Equation (15), is $\left({\widehat{u}}_{pqt}\right)$. Note that the Parzen window is applied to $\tau$, and the data dependent rule of Newey and West (1994) to $l_n$.

To test the cointegration relationship between the temperature anomaly distribution and the TRF level, the variable addition test (VAT) developed by Park (1990), is employed. To implement the test, two types of instruments are considered to the regression Equation (16) as

$$\begin{array}{c}\hfill y_t^{*}=z_{pqt}^{{*}^{\prime }}{\alpha }_{pq}+s_t^{i^{\prime }}\zeta +u_{pqt}^{*},\mathrm{where}i=1,2.\end{array}$$

(18)

where

$$\begin{array}{cc}\hfill s_t^1& =(t^2,t^3)\hfill \\ \hfill s_t^2& =sin\left(13.75\left(\frac{t}{T}\right)+0.21\right),\mathrm{where}\frac{t}{T}\in [0,1]\hfill \end{array}$$

The first instrument, $s_t^1$, is the two-dimensional superfluous regressors, quadratic and cubic trends, and the second instrument, $s_t^2$, is a single sine function to capture the Oceanic Multidecadal Oscillation (Miller and Nam 2020).7

The null hypothesis $H_0:\zeta =0$ cannot be rejected if there is a cointegration between the variables in Equation (16) in the sense that the added superfluous regressors are irrelevant variables. Under the alternative of spurious regression, superfluous regressors would be significant in the sense that the error contains nonstationary trends. Therefore, the null hypothesis would be rejected if there is no cointegration in Equation (16) (i.e., spurious regression). The reason for consideration of the second instrument, is that the 65–80 year oceanic cycle does not strictly follow the covariance stationary process. That is, the test is implemented for the case that the time-varying mean violates the stationarity of the error, and therefore the cointegrating relationship. Note that Estrada et al. (2013b) presented their empirical analysis with the Atlantic Multidecadal Oscillation-filtered Global and Northern Hemisphere mean temperature anomalies.

The main benefit of the VAT test is that it is easy to implement in addition to the robustness of many possible specification errors. Moreover, the VAT test prevents the nonlinear cointegration model from including too many nonlinear temperature terms. The test statistics of the VAT is given by

$$\mathrm{VAT}=\frac{RSS-RS{S}^k}{{\widehat{\omega }}_{nk}^{*2}}{⟶}^d{\chi }_k^2\mathrm{as}n\to \infty$$

where RSS and RSS${}^k$ are the sum of squared residuals from Equations (16) and (18), respectively, and ${\widehat{\omega }}_{nk}^{*2}$ is the consistent long-run variance estimate for the CCR errors from Equation (18). Note that the VAT statistics is the cointegration test for the linear regression model and it is feasible in the presented model because the nonlinear cointegration model can be reduced to the linear model of Equation (15).

Further, the model specification test is implemented using the Wald test based on the CCR estimator. Note that Wald statistics follow the asymptotically chi-square distribution because the limiting distribution of the CCR estimator follows mixed-normal. To do so, the mean-level regression is considered (i.e., $p=1$ and $q=0$) as a restricted model. For an unrestricted model, the CCR regression model is considered at the optimal p and q order, which is selected according to AIC and BIC criteria. Then, the Wald statistics based on the CCR estimator is

$$\begin{array}{c}\hfill W=\frac{RS{S}^1-RS{S}^2}{{\widehat{\omega }}_{nk}^{*2}}{⟶}^d{\chi }_{p+2q-1}^2\mathrm{as}n\to \infty \end{array}$$

where $RS{S}^1$ and $RS{S}^2$ are the sum of squared residuals from a restricted model and a minimum IC-based model, respectively, and ${\widehat{\omega }}_{nk}^{*2}$ is the consistent long-run variance estimate for the CCR errors from a minimum IC-based model. Notably, the rejection of the Wald test indicates the statistical significance of the specification error.

4. Data

The following data sources are employed for the Global, Northern Hemisphere, and Southern Hemisphere temperature anomalies and the TRF variable. As in Chang et al. (2020) who provided the technical background for the nonstationarity of global temperature anomaly distributions, the HadCRUT4 temperature anomaly data from 1850 to 2015 are employed for the Global, Northern Hemisphere, and Southern Hemisphere temperature anomaly data (Morice et al. 2012).8 Basically, temperature anomaly data from HadCRUT4 and the Goddard Institute for Space Studies (GISS) are generated from the same raw dataset. However, their treatments on the same raw dataset are different. Specifically, the GISS dataset uses an interpolated sea-surface temperature analysis by filling the empty grid boxes in the sea-surface area, while HadCRUT4 does not attempt to calculate the values for empty grid boxes. In light of this, HadCRUT4 would understate the effect of Arctic temperature anomalies, in which the warming has been significant over the past decades.9

Since the temperatures in land stations are measured at various elevations and since different countries employ different methods, the Global, Northern Hemisphere, and Southern Hemisphere temperature data are expressed as anomalies in degrees Celsius from the monthly temperature average from 1961 to 1990, which is known as the “zero-base” (i.e., climatological normal temperature) period. Note that since the number of stations and the methods of temperature measurement are different across grid boxes, calculating deviation from the zero-base may eliminate the heterogeneity across grid boxes over the entire space of the Earth. In this context, the temperature anomaly data are directly exploited for each grid box instead of recovering the actual temperature dataset.

Further, 99 percent of the total probability mass is exploited as the support of the temperature anomaly distribution at each time t because 0.5 percent of the probability mass at each end would be an adequate threshold to minimize the estimation errors induced by the boundary problem from the standard kernel density estimation technique. Figure 1 shows the time series of left and right ends of the chosen support of the global temperature anomaly distribution, indicating that non-common support would be necessary to effectively generate the distributions of temperature anomaly. Specifically, the right end of the support of the global temperature anomaly distribution reveals an increasing trend, implying the warming phenomena.

For the radiative forcing variable, TRF data from 1850 to 2015 (Hansen et al. 2017) are employed, which represent the sum of anthropogenic forcing and natural variability. Specifically, TRF is the sum of well-mixed GHG (CO${}_2$, CH${}_4$, N${}_2$O, and CFCs), ozone, surface albedo and tropospheric aerosols, and solar irradiance.10 Figure 2, Figure 3 and Figure 4 provide detailed information of the temperature anomaly distribution for the Globe, Northern Hemisphere, and Southern Hemisphere and TRF. Specifically, the top-left panels of Figure 2, Figure 3 and Figure 4 show the temperature anomaly distribution generated by the procedure of Chang et al. (2020), and the top-right panels of the figures illustrate the graphical comparison between the first moment estimated from the generated temperature anomaly distribution and the web-posted mean temperature anomalies from the GISS11 and HadCRUT4 (the median of the 100 ensemble member time series)12 websites for the Globe, Northern Hemisphere, and Southern Hemisphere.

Indeed, 99 percent of probability mass in estimating the temperature anomaly distribution provides a good approximation for estimating the first moment (i.e., mean temperature anomaly) at each time t, in the sense that estimated mean temperature anomalies are similar to mean temperature anomalies widely-used by climate scientists. Note that the GISS surface temperature data are expressed as an anomaly in degrees Celsius with the base period 1951–1980, which is available after the year 1880. As mentioned earlier, the differences in the mean temperature anomalies between the HadCRUT4 and GISS sources are greater for the Southern Hemisphere, in which the GISS team interpolates the sea-surface temperature data.

The bottom panels of Figure 2, Figure 3 and Figure 4 provide the first four central moments of the generated temperature anomaly distribution. While the variances of the estimated temperature anomaly distribution have decreased, roughly similar to Chang et al. (2020), the means, the skewness, and the kurtosis of the estimated temperature anomaly distribution have increased. In particular, the skewness appears to have increased from negative to positive. These statistical facts imply that the temperature anomaly distributions have been concentrated around their increasing means, and the probabilities of extremely positive temperature anomalies have increased. Not surprisingly, the decreasing variances of the temperature anomaly distribution would be on the same lines with the movements of the other central moments. Lastly, the calculated mean temperature anomalies are compared with the TRF variable, clearly showing that they moved together for the last 165 years.

5. Empirical Analysis

Throughout this study, the statistical testing result of the unit-root type nonstationarity of Chang et al. (2020) was followed. Specifically, the estimated persistence of the global mean temperature anomaly was closer to a stochastic trend, but not high enough to a deterministic trend, implying that no linear deterministic trend has been detected. As some econometricians discover a broken deterministic trend from the global mean temperature anomaly (Gay-Garcia et al. 2009; Estrada et al. 2013b, inter alia), moreover, Gao and Hawthorne (2006) attempt to estimate a general deterministic trend by allowing the flexible nonlinearity in the deterministic trend component. However, Chang et al. (2020) argue that a deterministic trend with the excessive nonlinearity and variability could be better expressed as a stochastic trend. Subsequently, the failure of rejection of the cointegration test implies that the TRF variable shares a stochastic trend with the global temperature anomaly distribution.

In the literature, climate sensitivity for the Globe is estimated as the value 0.43 ${}^{°}$C/(W/m${}^2$) and 0.35 ${}^{°}$C/(W/m${}^2$) with AMO-unfiltered HadCRUT4 and NASA dataset, respectively (Estrada et al. 2013b).13 As expected, linear climate sensitivity is estimated as the value 0.435 ${}^{°}$C/(W/m${}^2$), which is a similar value to the Global case.

Table 1. Estimation results.

Global Temperature Anomaly
p	q	${\mathit{D}}^{\mathbf{1}}$	Error(r)	${\mathit{b}}_{\mathbf{1}}^{\mathit{CCR}}$	T-Stat	AIC	BIC	${\mathit{VAT}}^{\mathbf{1}}$Stat	${\mathit{VAT}}^{\mathbf{2}}$Stat	WALD Stat
1	0	0.435	0.000	0.435	20.191	−217.336	−211.112	1.303	42.614	0.000
2	1	0.380	0.055	0.398	18.085	−244.910	−229.350	2.998	22.737	22.962
Northern Hemisphere Temperature Anomaly
p	q	${\mathit{D}}^{\mathbf{1}}$	Error(r)	${\mathit{b}}_{\mathbf{1}}^{\mathit{CCR}}$	T-Stat	AIC	BIC	${\mathit{VAT}}^{\mathbf{1}}$Stat	${\mathit{VAT}}^{\mathbf{2}}$Stat	WALD Stat
1	0	0.503	0.000	0.503	16.427	−95.034	−88.810	1.229	53.397	0.000
2	1	0.395	0.108	0.420	13.289	−125.168	−109.608	0.081	21.760	28.836
Southern Hemisphere Temperature Anomaly
p	q	${\mathit{D}}^{\mathbf{1}}$	Error(r)	${\mathit{b}}_{\mathbf{1}}^{\mathit{CCR}}$	T-Stat	AIC	BIC	${\mathit{VAT}}^{\mathbf{1}}$Stat	${\mathit{VAT}}^{\mathbf{2}}$Stat	WALD Stat
1	0	0.367	0.000	0.367	19.058	−256.415	−250.191	3.846	8.214	0.000
2	0	0.327	0.040	0.336	15.553	−268.311	−258.975	0.592	12.040	7.708

provides the estimation result of Equation (13) with a derivative of interest. Based on AIC and BIC criteria, the optimal models for the Globe, Northern Hemisphere, and Southern Hemisphere are $(p=2,q=1),(p=2,q=1)$, and $(p=2,q=0)$, respectively. Moreover, the first VAT test statistics (VAT${}^1$) for the Globe, Northern Hemisphere, and Southern Hemisphere cases as well as for linear/optimally chosen nonlinear models, indicate that all considered models are authentic, supporting the cointegration technique using the CCR methodology. However, the second VAT test statistics (VAT${}^2$) indicate that all regression models could be spurious if we assume the strict stationarity in the residual. Notice that the second VAT statistics decreased by considering the nonlinear temperature term for the Globe and Northern Hemisphere cases, indicating that the nonlinear temperature term $s_t^2$ in Equation (18) plays a role in explaining the oceanic multidecadal oscillation.

The nonlinear estimator $D^1$ represents climate sensitivity that considers all nonlinear effects across the observed temperature anomalies. The resulting climate sensitivity, $D^1$, of the nonlinear cointegration model is provided as the value 0.380 for the Globe, indicating that the global mean temperature anomaly increases by $0.380{}^{°}$C when TRF increases by 1 W/m${}^2$. In the meantime, the misspecification error of the linear model is greatest for the Northern Hemisphere (0.1077 ${}^{°}$C/(W/m${}^2$)), and lowest for the Southern Hemisphere (0.0402 ${}^{°}$C/(W/m${}^2$)). The Wald test decisively rejects the null of no statistical significance of the nonlinear temperature term for the Globe, Northern Hemisphere, and Southern Hemisphere cases (p-values are less than 0.01).

As shown by Figure 2, Figure 3 and Figure 4, TRF has affected not only the mean temperature anomaly, but also variance, skewness, and kurtosis temperature anomalies for approximately 150 years. This implies that the Earth’s surface temperature has been affected by human and natural forcing in a spatially heterogenous manner. The nonlinear cointegration model enables us to estimate the true effect by including the spatial distributions of temperature anomalies in the model, showing that the change in the global mean temperature anomaly associated with TRF (0.380 ${}^{°}$C/(W/m${}^2$)) would be less than what we have observed (0.435 ${}^{°}$C/(W/m${}^2$)).

In the literature, the transient climate response is often calculated as the global mean temperature response in ${}^{°}$C to a doubling of atmospheric CO${}_2$ from pre-industrial level by an increase of 1 percent per year. Schwartz (2012) estimates the value 3.71 W/m${}^2$ of the TRF level at the time when the atmospheric CO${}_2$ is doubled from the pre-industrial level. Subsequently, the transient climate response is calculated as the value 1.410 ${}^{°}$C (=3.71 × 0.380 ${}^{°}$C/(W/m${}^2$)), which is 0.204 ${}^{°}$C lower than the estimated transient climate response from the linear model (1.614 ${}^{°}$C). Given that the global mean temperature anomaly was $-0.{40}^{°}$C in 1850, the predicted global mean temperature anomaly associated with a doubled atmospheric CO${}_2$ level is 1.01 ${}^{°}$C.

Climatologically, the TRF response function at temperature anomaly $r_0$, $D^2\left(r_0\right)$, indicates a change in the mean temperature anomaly, if net incoming absorbed radiation is solely determined by temperature anomaly $r_0$. As such, the nonlinear effect, in addition to the mean effect (or linear effect), would provide the temperature-dependent TRF effect on the mean temperature anomaly. The nonlinear effects of the estimated climate sensitivity, ${\widehat{D}}^1$, are calculated as values −0.018 ${}^{°}$C/(W/m${}^2$), −0.025 ${}^{°}$C/(W/m${}^2$), and −0.009 ${}^{°}$C/(W/m${}^2$) for the Globe, Northern Hemisphere, and Southern Hemisphere (i.e., the term, $\partial {\int }_{D_r}B\left(r\right)f_t\left(r\right)dr/\partial h_t$ in Equation (10)). Therefore, we may conclude that the nonlinearity of the relationship induced by the spatial heterogeneity was strongest for the Northern Hemisphere.

In Figure 5, both the net incoming absorbed radiation, $B\left(r\right)$, and the derivative of the temperature anomaly distribution with respect to TRF, $d_T\left(r\right)$, are presented for the Globe, Northern Hemisphere, and Southern Hemisphere. Note that the domain of the temperature anomaly is shortened below −1.0 ${}^{°}$C and above 2.0 ${}^{°}$C for a practical purpose. The left panels of Figure 5 illustrate that the nonlinear effect becomes more significant in the opposite direction as the temperature anomaly approaches the left and right ends. As mentioned in Section 3, more specifically, $B\left(r\right)$ estimates the effect of net incoming absorbed radiation of the regions that are represented by the temperature anomaly r.

Note that the popular term, “polar amplification,” states that the low (high) temperature anomaly region, which mainly represents higher (lower) latitude areas, is related to the mean temperature anomaly in a positive (negative) direction (Boer and Yu 2003). Consistent with this notion, a positive (negative) value of the net incoming absorbed radiation term, $B\left(r\right)$ implies that the TRF effect under linearity would be underestimated (overestimated). In other words, the TRF effect under linearity should be amplified for the negative anomalies (high latitude areas) for the Globe and Northern Hemisphere, and be attenuated for the positive anomalies (low latitude areas) for the Globe, Northern Hemisphere, and Southern Hemisphere.

By considering the nonlinear temperature term with the derivative of the temperature anomaly distribution with respect to TRF, $B\left(r\right)d_T\left(r\right)$, the nonlinear TRF effect would be obtained by amplifying or attenuating the linear TRF effect (see Equation (11)). Since the changes in probability with respect to the change in TRF (i.e., $d_T\left(r\right)$) are very small at extreme temperature anomalies, the nonlinear effects at extreme temperature anomalies could be much less than those at low or high temperature anomalies.

Figure 6 illustrates the TRF response function for the Globe, Northern Hemisphere, and Southern Hemisphere, which is the estimator $D^2\left(r\right)$ for $r\in D_r$ in Equation (11). The reference (dashed blue) line represents the linear TRF response function, postulating that the TRF effect on the mean temperature anomaly is constant across temperature anomalies. That is, it illustrates the TRF effect on the mean temperature anomaly without the nonlinear effect. For the Globe, the linear term provides an estimate of the constant climate sensitivity as the value 0.40 ${}^{°}$C/(W/m${}^2$), and the nonlinear TRF response function intersects with the linear TRF response function at two anomalies, −0.14 ${}^{°}$C and 0.18 ${}^{°}$C. This implies that the linear TRF response function underestimates the true effect between these two anomalies.

In Figure 6, more importantly, the nonlinear TRF response function (i.e., nonlinear climate sensitivity) provides an interpretation of the relative magnitude of the TRF effect on the mean temperature anomaly. In particular, the greatest temperature-dependent TRF effect on the global mean temperature anomaly is estimated as the value 0.398 ${}^{°}$C/(W/m${}^2$), if net incoming absorbed radiation is solely determined by a temperature anomaly, 0.015 ${}^{°}$C. That is, the spatial contribution to the change in the global mean temperature anomaly is strongest for the areas in which their temperature levels correspond to $0.015{}^{°}$C. Such spatially heterogeneous contributions could be better understood by incorporating the higher-order moments of spatial distributions of temperature anomaly in the nonlinear cointegration model.

Moreover, the smoothed global mean temperature anomaly using a polynomial function was approximately 0.01 ${}^{°}$C in 1974, implying that the TRF effect on the global mean temperature anomaly would be strongest in 1974. Note that the sum of nonlinear effects of estimator $D^2\left(r\right)$ across temperature anomaly in Equation (11) is equal to the nonlinear effect of estimator $D^1$ in Equation (10). In this light, the high degrees of nonlinearity could be estimated by adding all nonlinear effects, although the magnitude of nonlinear effect at each temperature anomaly is small. Further note that the nonlinear TRF effect for the Northern Hemisphere case shows a similar pattern to the Global case. However, its effect is greater than the Global case because the speed of global warming is faster in the Northern Hemisphere than in the Southern Hemisphere.

To elaborate on this discussion, Figure 7 provides the nonlinear TRF response function over time. Specifically, I postulate that the spatial regions where provide the biggest contribution on the change of mean temperature anomaly at time t is located at the areas where have a record of the mean temperature anomaly at time t. This is a reasonable situation, in the sense that the mean value of spatial distribution would be largely determined by the spatial areas where have a similar or closer mean value. Consequently, the value of the nonlinear TRF response function at time t indicates the relative magnitude of the TRF effect when the biggest spatial contribution is evaluated at the mean temperature anomaly at time t.14

The TRF response functions presented in Figure 7 imply that the warming speed of the Globe and Northern Hemisphere has decreased since 1980. However, the speed would start to increase when the mean temperature anomaly reaches approximately 0.8 ${}^{°}$C. Notably, the presented TRF response function for the Southern Hemisphere has decreased since 1910. Similar to the Globe and Northern Hemisphere cases, however, the warming speed of the Southern Hemisphere would start to increase when the mean temperature anomaly reaches approximately 0.7 ${}^{°}$C, which is much earlier compared to other cases. Based on the information provided by the HadCRUT4 website,15 the mean temperature anomalies for the Globe, Northern Hemisphere, and Southern Hemisphere in 2019 were estimated at 0.736 ${}^{°}$C, 0.972 ${}^{°}$C, and 0.502 ${}^{°}$C, respectively, indicating that the warming speed of the Northern Hemisphere has already started to increase.

6. Conclusions

In this paper, I proposed the nonlinear cointegration model based on the well-known EBCM. For this, the nonlinear cointegrating regression of the mean temperature anomaly for the Globe, Northern Hemisphere, and Southern Hemisphere, were estimated using the spatial distributions of temperature anomalies. Subsequently, the nonlinear TRF effect on the mean temperature anomaly was estimated, suggesting that the TRF effect on the mean temperature anomaly would be temperature-dependent for the Globe, Northern Hemisphere, and Southern Hemisphere. Graphically, the TRF response function has a flexible shape to represent the change in the mean temperature anomaly when net incoming absorbed radiation is hypothetically determined at some temperature anomaly.

Statistically, the linear model fails to take into account the net incoming absorbed radiation term, which would cause the slope estimator to be invalid. Considering the functional form of net incoming absorbed radiation, the proposed nonlinear cointegration model shows the reasonable nonlinear dependence structure between the mean temperature anomaly and TRF. Specifically, climate sensitivity is estimated to be temperature-dependent (or spatially heterogenous), providing that the estimated nonlinear climate sensitivity was highest in the mid-1970s for the Globe. Moreover, the TRF effect on the mean temperature anomaly, which considers all nonlinear effects, is less than the estimate of the linear climate sensitivity provided by the literature. Lastly, the statistical testing results indicate that the linear model possesses a significant misspecification error for the Globe, Northern Hemisphere, and Southern Hemisphere cases.

The next step regarding this research refers to climate variability. As emphasized by Brock et al. (2013), spatial complexity on Earth is a big challenge for analyzing temperature series at hemispheric scales. The complexity inherited from spatial diversity produces hardly explainable natural variability through the observed data. At least to date, the well-known inter-annual global variability or hemispheric climate variability is regarded as the El Niño/Southern Oscillation (ENSO), the North Atlantic Oscillation (NAO), and the Atlantic Multidecadal Oscillation (AMO). Among them, the most influential natural variation for the attribution study is the AMO, which could distort the long-term global warming trend. Specifically, the large ocean-atmosphere cycle over the North Atlantic, which is defined as approximately 60 to 90 years of low-frequency patterns of sea surface temperature variability, explains the larger variability in Northern Hemisphere temperatures and therefore, globally.

In particular, Estrada et al. (2013b) identified its difficulty when they conducted an attribution study. They argued that the detrended Global and Northern Hemisphere temperatures with forcing variables could be further explained by the AMO and therefore the difference between the dates of a structural break could be explained as well. In light of this, Estrada et al. (2013b) filtered the AMO information from the Global and Northern Hemisphere mean temperature anomalies to estimate constant climate sensitivity. To put these factors into perspective, it is worth estimating nonlinear climate sensitivity after extracting the major climate variabilities.

Moreover, the uncertainty of the estimator of the TRF response function needs to be evaluated. As the estimation is implemented via two-step approach, it is difficult to evaluate the bootstrapping confidence bands, in the sense that the estimator’s uncertainty comes from both steps. We may estimate the uncertainty by fixing the point estimate of the first step as a constant. Since the first step is the nonstationary nonlinear regression model, however, the uncertainty of the first step’s estimator would be larger than that of the second step’s estimator. I leave these tasks to future research.