Detection and Attribution of Climate Change Through Econometric Methods

Detection and attribution of climate change through econometric methods Abstract An ever growing body of evidence regarding observed changes in the climate system has been gathered over the last three decades and large modeling efforts have been carried to explore how climate may evolve during the present century (IPCC, 2001; 2007a). The impacts from both observed weather and climate endured during the 20th century and the magnitude of the potential future impacts of climate change have made this phenomenon of high interest for policy-makers and the society at large (IPCC, 2007b). Two fundamental questions arise for understanding the nature of this problem and the appropriate strategies to address it: is there a long-term warming signal in the observed climate, or is it the product of natural variability alone? if so, how much of this warming signal can be attributed to anthropogenic activities? As discussed in this review, these questions are intrinsically related to the study of the time-series properties of climate and radiative forcing variables and of the existence of common features such as secular co-movements. This paper presents a brief summary of how detection and attribution studies have evolved in the climate change literature and an overview of the time series and econometric methods that have been applied for these purposes. 1 A brief review of time series based detection and attribution of climate change. We first review work related to the detection and attribution of climate change based on time series techniques. We start with early studies related to the statistical properties of temperatures. We then review the work using cointegration techniques assuming temperatures and forcing variables to be unit root processes. Finally, we cover recent work arguing that the series of interest are better described as stationary fluctuations around a trend function with changes in slope. In such cases, attribution consists in establishing the existence of a common non-linear trend between temperatures and radiative forcings. All technical concepts are presented in subsequent sections. 1.1 Early studies on the detection and attribution of climate change. Since the late 1970 different research groups have published different estimates of global and hemispheric temperatures based on the available observational data. Large improvements in information recovery and processing/analysis techniques have occurred in the last decades, leading to longer temperature records based on more accurate estimates and to a better assessment of their uncertainties (e.g., Hansen et al., 2010; Morice et al., 2012; IPCC, 2007a). These efforts have permitted the analysis of observed global warming through econometric methods that is reviewed here. Two main types of stochastic processes have been proposed to represent global temperature series: differencestationary and trend-stationary. In earlier studies of the time-series properties of observed global temperatures, the presence of a deterministic trend was interpreted as evidence for a long-term, human-induced, global warming process. In contrast, under a stochastic trend, the recent warming trend was interpreted as a manifestation of the low-frequency variability and high persistence of climate, indicating that it is temporal and should not be expected to continue in the long-term. In this case, no human interference with the climate system had occurred and no mitigation strategies to avoid further changes are justified. Some recent publications continue to propose this interpretation (Mills, 2010a,b). For instance, based on standard Augmented Dickey-Fuller (ADF) test, Galbraith and Green (1992) found evidence in favor of global temperature series being a trendstationary process with a positive slope. Their main conclusions supported the idea of a lasting global warming process. In the same vein, Bloomfield (1992) concluded that statistical models consisting of a trend plus serially correlated noise may be fitted to temperature data and estimated a constant growth rate ranging from .2 to .8 degree Celsius per century, providing evidence for the presence of global warming during the last century. Zheng and Basher (1999) reported evidence reinforcing the view that global temperatures are affected by a long-term warming that is not of natural origin or at least that the existence of a linear trend cannot be ruled out. 1 In contrast, some argued that temperature data are better described by ARIMA (integrated autoregressive model average) models, e.g., Woodward and Gray (1993). Their conclusions gave support to the hypothesis that the upward trend in temperature was the product of random fluctuations and that no global warming was in effect. Moreover, Woodward and Gray (1995) fitted an ARIMA model to three different reconstructions of global temperature series, concluding that the series are better modeled as integrated processes, although a deterministic trend could also be present. Further evidence was provided by Kärner (1996) who argued that a unit root process could explain the trends and the variability of temperatures over the last century. Short-term trends are inevitable and may have little in common with increasing CO2 concentrations. Such findings increased the support in favor of the unpredictability and randomness of temperature patterns. To a lesser extent, fractionally integration has also been advanced as a possible representation of the data generating process of global and hemispheric temperature series (Gil-Alana, 2008a,b; Rea et al., 2011; Mann, 2011; Qu, 2011). Nevertheless, it is important to note that the use of fractionally integrated models have been limited to the detection of climate change and the notion of fractional cointegration has not yet been applied to investigate the attribution of climate change. 1.2 Cointegration approach: testing for a common stochastic trend. In light of such results, the Intergovernmental Panel on Climate Change (IPCC; Folland et al., 1992) concluded that present statistical tests were not able to resolve the question of statistically significant relationships given the differences in the time-series properties of possibly trend-stationary temperature series and the concentrations of greenhouse gases series assumed to be integrated processes. This gave rise to a second stage in the time series analysis of global and hemispheric temperatures, in which efforts concentrated on the issue of the attribution of climate change. The application of cointegration techniques to global and hemispheric temperature series and to radiative forcing variables provided a breakthrough on this issue (Richards, 1993; Tol and de Vos, 1993; Tol, 1994; Stern and Kaufmann 1997a,b). These papers changed the view of how a stochastic trend in temperature series should be interpreted. The presence of a unit root was no longer seen as evidence against anthropogenic global warming, but rather as confirmation that mean global temperatures were driven by anthropogenic forcing. According to Stern and Kaufmann (1997b) the evidence of cointegration provided statistically rigorous and direct attribution of climate change to anthropogenic activity. Furthermore, it provided an alternative method for estimating the climate sensitivity (i.e., the amount of warming expected as a result of a doubling of the atmospheric carbon dioxide concentrations; see Skinner, 2012) from observed data. The use of cointegration techniques seemed to have ended the debate on the data generating 2 process of global and hemispheric temperatures and on the statistical approach to estimate their relationship with forcing variables. Therefore, this representation became widely accepted in the climate change literature and has been used also to model the long-run relationships of different climate variables, e.g., global temperatures and sea level (Schmidth et al., 2012). Nevertheless, the problem of identifying the order of integration of temperature and radiative forcing variables was not solved. This lead to publications arguing different orders of integration without following adequate procedures: 1) according to unit root tests, temperatures were trend stationary processes [I(0)] but it was concluded that they must be integrated of order one [I(1)] because a cointegration relationship between them and radiative forcing (identified as I(1)) could be found (Stern and Kauffman, 1997a); 2) temperatures are I(1) based on evidence of a standard ADF unit root test applied to a longer sample which includes more recent records of global temperature values (Kaufmann et al., 2006a,b); 3) temperatures are integrated of order two [I(2)], meaning that they are characterized by having two independent stochastic trends (Stern and Kauffman, 1997b). Some authors have applied the concept of polynomial cointegration 1 , leading to opposing conclusions regarding the role of anthropogenic forcing in the observed warming (Beenstock et al., 2012; Liu and Rodriguez, 2005) and questioning the validity of previous results found within the cointegration framework assuming I(1) processes. Not being able to identify the order of integration of a time series has strong implications for cointegration analysis. In this case, cointegration techniques are not reliable (e.g., Gonzalo and Lee, 1998; Elliot, 1998; Leybourne and Newbold, 2003). It has been shown in the econometric literature that cointegration methods tend to find spurious cointegration, with probability approaching one asymptotically, when the correct order of integration is uncertain and when there are structural breaks in the trend function (Gonzalo and Lee, 1998). As discussed in Estrada et al. (2010), the lack of a proper analysis of the time-series properties of climate variables is a fundamental problem that bedevils the application of cointegration techniques for the attribution of climate change studies, as also is the lack of physical mechanisms that would support the existence of a unit root process in temperature series. Although cointegration offered a formal way to study the attribution of climate change, results reported in the literature pointed to the need of revising the assumptions regarding the univariate time-series properties of temperature and radiative forcing variables. This is as a necessary first step for the multivariate modeling that could be employed to investigate the attribution of observed climate change. 1 Multicointegration extends the cointegration methodology to integrated variables of higher order than one and to variables integrated of different order. 3 1.3 Trend stationarity and co-trending: testing for common nonlinear secular trends. The presence of breaks in global and hemispheric temperatures has been discussed extensively (e.g., Seidel and Lanzante, 2004; IPCC, 2007a). The rate of warming of the 20th century can hardly be considered constant and therefore a linear trend provides an inadequate representation to describe the secular movement of global and hemispheric temperatures. Consequently, research on the properties of temperature series extended the deterministic linear trend plus stationary noise to more realistic specifications of the trend function that allow the presence of nonlinearities. Using recent time-series techniques and unit root tests that allow for the presence of a structural break (Kim and Perron, 2009; Perron, 1997), Gay et al. (2009) showed that global and hemispheric annual temperature series are better represented as trend stationary processes with an “exogenous and permanent” change in the rate of growth that cannot be interpreted as part of the natural variability occurring during the 20th century. Their main results can be summarized as follows: 1) there is strong statistical evidence against the existence of a stochastic trend in global and hemispheric temperatures: their data-generating process can adequately be described as trend stationary with a single change in slope occurring at different dates and with different magnitudes for the various series. In general, the results show evidence for a “two stages” warming trend, the first a moderate one and the other rapid and of much larger magnitude; 2) the time-series properties of annual global and hemispheric temperatures suggest that, at least for these spatial scales and data frequency, climate change has manifested itself as a “change-in-the-mean phenomenon”, variability has not changed with climate change; 3) in terms of Article 2 of the Framework Convention on Climate Change significant anthropogenic interference with the climate system has already occurred and; 4) inference methods, such as cointegration, which assume that temperatures are integrated processes should be revised. This meant that a new approach to relate temperature and radiative forcing series was needed to investigate the anthropogenic contribution to the warming. Other authors also offered evidence supporting trend stationarity around nonlinear trends as a better representation of the data generating process of temperature series. Harvey and Mills (2001, 2002) showed that global and hemispheric temperatures can be better represented as stationary processes around one or two smooth transitions in the linear trend. Holt and Teräsvirta (2012) provided further evidence in favor of trend stationarity by applying a unit root test that uses Fourier series to approximate any unspecified shifts in the trend (Becker et al., 2006). Based on the work of Perron (1989), Gay et al. (2009) provided an explanation for the inability of standard unit root tests to determine the order of integration of temperature series. Using proper tests for structural breaks (Perron and Yabu, 2009a), they provided strong evidence for the existence of a break in the slope of the trend function irrespective of the order of integration of 4 temperature series. The recursive ADF tests shown in Gay et al. (2009; Fig. 1) suggest that the existence of a break can explain why standard unit root tests find temperatures either integrated or trend-stationary depending on whether the sample includes observations from the last two to three decades of the 20th century. The effects of breaks in the trend function on the behavior of unit root tests was established by Perron (1989) who showed that the sum of the autoregressive coefficients is highly biased towards unity if there is a shift in the trend function. In this case, the unit root null is hardly rejected even if the series is composed of white noise disturbances around the trend. If the break occurs in the slope of the trend, unit root tests are not consistent, i.e., the null hypothesis of a unit root cannot be rejected even asymptotically. Gay et al. (2009) generated a new phase in the debate of identifying the type of data generating process for temperature series. Shortly after its publication, Kaufmann et al. (2010) argued in favor of the cointegration approach mainly based on the potential usefulness of this technique for attribution studies as, if this representation was indeed adequate, it would offer “the possibility of greater insights regarding the potential causes of climate change and efforts to slow its progression”. By means of structural time-series models, Mills (2010a,b) supports the unit root hypothesis, but reaches very different results from those in Stern and Kaufmann (1997a,b), Kaufmann and Stern (1997) and Kaufmann et al (2006a,b; 2010) regarding the detection and attribution of climate change. According to this author, global temperature is better represented as a simple random walk, where the lack of a drift indicates that there is no long-term warming trend, and that what was observed during the 20th century were realizations of natural variability. In the light of the results of Gay et al. (2009), Harvey and Mills (2001; 2002) and Holt and Teräsvirta (2012), the difference in the order of integration between temperatures (trend stationary) and radiative forcing variables, commonly assumed as integrated, seemed to preclude the possibility of establishing empirical evidence on attribution through current statistical methods. Estrada et al. (2013a) investigated the time-series properties of radiative forcing variables and found that like temperature series they are trend stationary processes with time ordered breaks in the slope of their trend function. A new approach to conduct detection and attribution studies was introduced and it was shown that global temperatures and radiative forcing share a common nonlinear secular movement. This approach is based on “common features” concepts such as co-trending and co-breaking (Engle and Kozicki, 1993) and uses the nonparametric nonlinear cotrending test of Bierens (2000). Cointegration analysis is only one possibility for relating trends of non-stationary variables: relationships between nonstationary variables can be established when linear combinations of different time series cancel out common features such as trends and cobreaks in trend. It was shown that climate models’ simulations and observed global temperature 5 series share a common nonlinear trend that is imparted mainly by the radiative forcing of wellmixed greenhouse gases. This paper provided an approach that is both adequate for the time-series properties of temperature and radiative forcing variables and consistent with the physics of climate. Holt and Teräsvirta (2012) found Northern and Southern Hemispheres temperature series to be trend stationary processes and they use a Shifting-Mean Vector Autoregression (SM-VAR) to provide evidence in favor of a deterministic co-shifting between Northern and Southern temperature series. This modeling allows the possibility that shifts in the underlying data occurred gradually and that one or more of these shifts are common to both series. They argue that for hemispheric temperatures the evidence of co-shifting is particularly strong since the early 1980s. Recent advances in the understanding of low-frequency variability modes have shown that the 60- to 90-years natural oscillations produced by the Multidecadal Atlantic Oscillation (AMO) can mask or exaggerate the warming trend depending its phase (Wu et al., 2011; Swanson et al., 2009). Estrada et al. (2013b) showed that after filtering out the effects of AMO, the secular movement of global, northern and southern hemisphere temperatures is characterized by a common warming trend with a break in its slope near 1955, marking the onset of sustained global warming. Using co-trending analysis they show that the nonlinear warming trend is of anthropogenic origin mostly driven by the radiative forcing of well-mixed greenhouse gases. Also, the analysis of the common warming trend imparted by the radiative forcing indicate that human interventions contributed to slowing down global warming in two occasions: the first is related to the large socioeconomic disruptions caused by the two World Wars and the Great Crash which contributed to the cooling in the mid 20th century; and the second is associated to the Montreal Protocol and the technological change in agricultural production in Asia. Both of these changes in greenhouse emissions are major drivers of the slowdown in the warming that has been experienced since the late 1990s. 2 Difference stationary, trend stationary and fractionally integrated processes. Consider a time series yt with the following decomposition: yt = τ t +zt , where τ t is the deterministic trend function and zt is the noise component. A process is said to be integrated of order d or I(d) if the dth difference of the noise zt , ∆d zt = (1 − L)d zt , is stationary. To illustrate this concept consider an ARIMA process such that A(L)zt = B(L)et with et ∼ i.i.d. (0, σ 2e ), where A(L) and B(L) are polynomials in the lag operator L (defined as Lzt = zt−1 ). It is assumed that the roots of B(L) = 0 are outside the unit circle so that the process is identified, stationary and invertible. If a time series is stationary around an appropriately defined trend τ t its order of integration is zero or I(0). This occurs when all the roots of A(L) = 0 are outside the unit circle. The process is said to be I(1) if the deviations from the trend have to be differenced once to achieve stationarity, in 6 which case one of the roots of A(L) is one. Similarly, it is said to be I(2) if the deviations have to be differenced twice to achieve stationarity, in which case A(L) has two unit roots. 2.1 Unit roots. Consider the first-order autoregressive model zt = αzt−1 + et , et ∼ i.i.d. 0, σ 2 . (1) An example of a unit root process is when α = 1. Then, ∆zt = zt − zt−1 = et . The first difference of the process is i.i.d. This model has the following implications. First, each shock et has a longterm effect on the level of zt . To see this, write (1) with α = 1 as (by recursive substitution) zt = z0 + t j=1 ej . Since each shock has a permanent effect on future levels of zt , a 1% unexpected increase in zt today increases our predicted value of future zt ’s by 1% for all future periods. In this simple example, with et ∼ i.i.d., zt is called a random walk, the best predictor of zt tomorrow being zt today, i.e., E [zt+1 |Ft ] = E [zt |Ft ] + E [et+1 |Ft ] = 0, where Ft is all information known at time t, given that zt ∈ Ft and et+1 is whie noise. The most important implication of a unit root is the permanent effect of shocks. Assuming this type of process as a representation of global and hemispheric temperatures implies that the secular movement of the series is determined by the sum of random shocks: all shocks have permanent effects on temperature series and even shocks in the distant past are as important as present variations to determine the current trend. The long-term forecast is always influenced by historical events, and temperature predictability is limited, even if forcing factors are held constant (Estrada et al., 2013a; Gay et al., 2009; 2007). The second implication is that the variance of zt increases with t since V ar (zt ) = V ar( t j=1 ej ) = tσ 2 if z0 is fixed. This is a nonstationary process since its second moment depends on t. Hence, a unit root process is non-stationary in variance. The process can cross any line within a long enough period. 2.2 Difference versus trend-stationary models. The random walk model is quite restrictive. Most of the time, allowing for additional short-run correlation is needed; i.e., having ∆zt = vt , where vt is a stationary process exhibiting some correlation (without a unit root itself). For example, if vt is a MA(1), a moving-average of order one, vt = et + θet−1 . Then, zt = zt−1 + et + θet−1 . Suppose that there is a unit shock et = 1 today and none thereafter, then at time t, zt increases by 1 and by θ at time t + 1 and stays like that thereafter. So, the overall effect in the long run is to increase zt by (1 + θ). In general, we can specify vt as a stationary and invertible ARMA model (autoregressive moving-average model) A (L) vt = B (L) et with et ∼ i.i.d.(0, σ 2 ). Since vt is a stationary and invertible process, it can 7 be represented as an MA(∞), i.e., vt = A (L)−1 B (L) et ≡ ψ (L) et . Then, in analogy with the MA(1) case, the long-run effect of a unit shock et on the level of zt is given by the sum of the MA coefficients, i.e., 2.3 ∞ i=0 ψ i = ψ (1). Random walk with drift or DS process. A wide variety of time series have a tendency to show secular movements over time. Hence, it is common to specify the trend function as τ t = µ + βt. Here, β is the slope of the trend function, sometimes called the drift. Then, ∆yt = β + vt with A (L) vt = B (L) et , and yt = µ + βt + t j=1 vj . The trend function of yt is then composed of two parts: i) a deterministic part given by the drift term, ii) a stochastic part given by the permanent effect of each shock et on the level of yt . Since shocks have a permanent effect, they change our long term forecast of the level of the series which is a plausible definition of the trend of a series. In this case, the trend function shifts every period and is, hence, stochastic. For this reason, such models are called stochastic trend models (also called difference stationary models, DS). The trend function can be obtained using the so-called Beveridge-Nelson decomposition. Consider the moving-average representation of ∆yt given by ∆yt = β + ψ (L) et . Now decompose ψ (L) as ψ (L) = ψ (1) + (1 − L)ψ ∗ (L), where ψ∗ (L) = − ∞ ∗ i i=0 ψ i L with ψ ∗i = ∞ j=i+1 ψ j . Then, ∆yt = β + [ψ (1) + (1 − L)ψ ∗ (L)]et or yt = µ + βt + ψ (1) t j=1 ej + wt (2) where wt = ψ ∗ (L) et is a stationary random variable. A plausible definition of a trend function is the long term conditional forecast, E(yt+k |ys ; s ≤ t) for k large. Using (2), we have E(yt+k |ys ; s ≤ t) = µ + β(t + k) + ψ (1) t j=1 ej . since wt is stationary. Hence, we can view the component µ + βt + ψ (1) µ + βt the deterministic part and ψ (1) 2.4 t j=1 et , (3) t j=1 et as the trend with a scaled random walk, the stochastic part. The trend stationary model (TS). An alternative to model variables that increase over time is with a purely deterministic trend, i.e., yt = c + βt + wt , where C (L) wt = D (L) et with et ∼ i.i.d. 0, σ2 and all the roots of C(L) and D(L) outside the unit circle. This is called the TS model because the deviations from the trend function are transitory, i.e., the shocks et have no permanent long run effects. Note that in this case ψ (1) in (3) is 0 so that the trend has no stochastic component. Distinguishing between a TS and DS process is the so-called unit root testing problem. 8 2.5 Trend function with breaks. A particular case of the trend stationary process that has been discussed in the climate change literature is when breaks in the trend function are present. In general, the trend parameters and their structural changes need not to be assumed deterministic (Perron, 1989; Perron and Wada, 2009; see Perron, 2006 for a survey). In order to illustrate the class of models that applies in such cases, consider the framework offered in Perron and Wada (2009): yt = µt + β t t + zt where µt = µt−1 + vt and β t = β t−1 + ut . The intercept and slope of the trend function are time varying stochastic processes. The noise components are modeled as mixtures of normal distributions where the realizations from the variables are drawn from one of two normal distributions, one with high and the other with small or zero variance. These processes can be described as: ut = λt γ 1t +(1−λt )γ 2t and vt = κt δ 1t +(1−κt )δ 2t , where γ it ∼ i.i.d. N(0, σ2γi ), δit ∼ i.i.d. N(0, σ 2δi ) while λt and κt are Bernoulli variables that take value one with probability αλ and ακ , and value zero with probability (1 − αλ ) and (1 − ακ ), respectively. One can then obtain a model with infrequent changes in the slope and intercept parameters when αλ and ακ are close to one and σ 2γ1 and σ2δ1 are zero. If σ2γ1 > 0, there will be occasional changes in the slope, and correspondingly if σ2δ2 > 0 there will be infrequent changes in the intercept. When only one break occurs, it becomes difficult to model the change with a stochastic structure. Hence, the common approach in the literature has been to consider the change as being ‘exogenous’ in the sense of intervention analysis (Box and Tiao, 1975) and they are not explicitly modeled via a parametric stochastic structure. Under this parameterization, there are only some shocks that can change the long-term behavior of the time series, as opposed to the case of a unit root process where all shocks have long-term changes. In the climate context, long-term changes are not frequent in the scale of the sample under analysis and are produced by important changes in key external forcing factors such as Earth orbit changes, solar irradiance, and greenhouse gases concentrations (Estrada et al., 2013a; Gay et al., 2009). 2.6 Fractionally integrated processes. Fractionally integrated processes occur when the order of difference d needed to achieve stationarity is not necessarily an integer but can be any real value (Granger and Joyeux, 1980). The most commonly used fractionally integrated models are ARFIMA processes, which are a generalization of ARIMA models when d is not an integer. In this case, the fractional differencing term can be 9 written as an infinite order moving-average (MA) process using the binomial expansion: (1 − L)d = 1 − dL + where Γ(g) = ∞ g−1 −y e dy 0 y d(d − 1) d(d − 1)(d − 2) + + ... = 2! 3! ∞ k=0 Γ(k − d)Lk Γ(−d)Γ(k + 1) is the gamma function. The autocorrelation function of this process is given by ρτ = [Γ(1 − d)Γ(τ + d)]/[Γ(d)Γ(τ + 1 − d)] ≈ cτ 2d−1 for some c. Hence, the autocorrelation function decays at a slow hyperbolic rate when 0 < d < 0.5. For such processes, yt is stationary and invertible if the roots of A(L) and B(L) lie outside the unit circle and |d| < 0.5. If the integration order is 0 < d < 0.5, the process displays long memory (much longer persistence than that of an I(0) process with similar ARMA parameters). When the integration order is 0.5 < d < 1, the process is non stationary. Although its variance does not have a finite sum, the effect of a unit shock on the level of the series does decay to zero implying that the memory of the process is not infinite, as in the case of a I(1) process, and that it exhibits mean-reversion. In the case d ≥ 1 the process is not mean reverting and nonstationary (see Maddala and Kim, 1998; Baillie, 1996). 2.7 Examples. Figure 1 shows the global mean surface temperature (G) with the effects of the AMO filtered out (see Estrada et al., 2013b), the total radiative forcing (TRF), and three simulated temperature series representing: 1) a TS process with a break in its slope (Temp_TSB), calibrated from unfiltered G from the NASA dataset, and generated by T emp_T SBt = −0.3200 + 0.0039t + 0.0129DTt∗ + ut with ut = 0.3987ut−1 + et and et ∼ i.i.d. N (0, 0.0083), where DTt∗ = t − Tb if t > Tb and Tb = 1978; 2) a DS process (Temp_DS) generated by T emp_DSt = 0.01 + T emp_DSt−1 + ut with ut ∼ i.i.d. N (0, 0.01) and; 3) a stationary process (Temp_S) generated by T emp_St = 0.7T emp_St−1 + ut with ut ∼ i.i.d. N (0, 0.01). The filtered version of G is used since, as shown in the literature (Wu et al., 2011; Swanson et al., 2009; Estrada et al., 2013b), the AMO can mask or exaggerate the warming trend depending its phase. The unfiltered G as well as TRF are available from http://data.giss.nasa.gov/, while AMO is available from http://www.esrl.noaa.gov/. With the exception of AMO and Temp_S all series show a trending behaviour and visual inspection suggests the existence of a break in the slope of their trend functions, including Temp_DS and Temp_S even though by construction they do not include one. This illustrates the potential pitafalls of relying on visual inspection or curve fitting to indentify breaks. As will be shown, formal statistical tests reliably povide the correct results. 10 3 Testing for a unit root. The most commonly used unit root test is the ADF test (Dickey and Fuller, 1979; Said and Dickey, 1984). Leaving aside the deterministic components for a moment, suppose that the data generating process is an AR(p) process of known order p, i.e., p i=1 ai yt−i yt = + et or A(L)yt = et with A(L) = 1 − a1 L − ... − ap Lp . If the process has an autoregressive unit root, then A(L) = (1 − L)A∗ (L) where A∗ (L) is a polynomial of order p − 1 with all roots outside the p i=1 ai unit circle. Hence, A(1) = 0, which implies that p−1 i=1 di yt = αyt−1 + with α = p i=1 ai and di = − p j=i+1 ai . = 1. Now reparameterize the model as: △ yt−i + et (4) The test is based on estimating (4) and using the t- statistic for testing the null hypothesis that α = 1. More generally, in the presence of deterministic components denoted now by xt , the regression is yt = βxt + αyt−1 + p−1 i=1 di △ yt−i + et . Which deterministic components to include depends on the series analyzed. If the data are not trending (e.g., stratospheric aerosols, North Atlantic Oscillation, etc.), only a constant is included. If the data are trending (e.g., temperatures, greenhouse gases, etc.), a constant and a time trend need to be included. The importance of the correct specification of the deterministic components is discussed in Campbell and Perron (1991) and Perron (1988). The limit distribution of the t-statistic is non-standard and depends on the nature of the deterministic components included (see Dickey and Fuller, 1979). The tests remain valid under more general processes such as ARMA or linear processes. The idea is to approximate the process by a AR(k) for k suitably chosen so that the residuals are approximately uncorrelated. The validity of the t-statistic was analyzed by Said and Dickey (1984) when the noise is an ARMA(p, q) process of unknown order with i.i.d. innovations. Their conditions were relaxed by Ng and Perron (1995) and Chang and Park (2002). Several extensions have been proposed to reduce size distortions and improve the power of standard unit root tests such at the ADF and Phillips and Perron (1988) tests. Some relevant variants include the development of efficient unit root tests that use quasi-differentiation and detrending via a Generalized Least Squares (GLS) method to improve power against persistent local alternatives (Elliot et al., 1996; Ng and Perron, 2001) and the modification of information criteria (e.g., AIC, BIC) used to determine the appropriate lag length to correct for serial correlation in the residuals 11 (Ng and Perron, 2001). Here, we shall discuss the widely used tests of Ng and Perron (2001), which adopt the GLS-detrending approach of Elliot et al. (1996). Suppose that yt = γxt + vt , where vt = αvt−1 + ut . If the errors {ut } are i.i.d. N(0, σ2u ), then the most powerful test of the null hypothesis α = 1 versus the alternative hypothesis that α = ᾱ, is given by the likelihood ratio test which reduces to L(ᾱ, 1) = S(ᾱ)−S(1), where S(α) = inf γ with ytα = yt − αyt−1 and xαt = xt − αxt−1 for t = 2, ..., T T α α 2 t=1 (yt −γxt ) = and y1α = y1 , xα1 T GLS α 2 α xt ) t=1 (yt −γ̂ = x1 . The estimate γ̂ GLS is the GLS estimate of γ from the model yt = γxt + vt assuming that vt is an AR(1) with parameter α. Hence, the likelihood ratio test compares the sum of squared residuals from two GLS regressions with autoregressive parameters 1 and ᾱ. Assuming that the true value α = ᾱ, and using the local to unity parameterization ᾱ = 1 + c̄/T the local asymptotic power function of the LR test gives a local asymptotic power envelop, i.e., the maximum that can be achieved under normal errors. In practice the true value of α is unknown. A sensible recommendation by Elliot et al. (1996) is to choose ᾱ such that the asymptotic power is 50% evaluated at c̄. This implies c̄ = −7.0 when xt = 1 and c̄ = −13.5 when xt = (1, t). In the presence of serial correlation, a feasible point optimal test with the same asymptotic properties is PT = [S(ᾱ) − ᾱS(1)]/s2AR , with s2AR defined below. Denote the GLS detrended data by yt = yt − γ̂ GLS xt . Consider the ADF test constructed with GLS detrended data, i.e., the t-statistic for testing that α = 1 in the regression yt = αyt−1 + k i=1 bi ∆yt−i + etk (5) Elliot et al. (1996) found that it has an asymptotic power function very close to the maximum power possible under normal errors. Ng and Perron (2001) considered modified versions of the Phillips and Perron (1988) unit root tests and others suggested by Stock (1999). The MZα test constructed with GLS detrended data is M ZαGLS = [T −1 yT2 − s2AR ]/[2T −2 where s2AR is an autoregressive spectral density estimator given by s2AR = σ̂ 2ek /(1 − T 2 t=1 yt−1 ] k 2 i=1 b̂i ) T 2 t=k+1 êtk with b̂i and êtk obtained from (5). They also considered the al2 /s2 )1/2 and MZ GLS = MSB GLS ∗ MZ GLS . These tests MSB GLS = (T −2 Tt=1 yt−1 t α AR GLS collectively referred to as the M tests. The limit distributions of the tests are non- where σ̂2ek = T −1 ternative tests are standard but critical values are presented in Ng and Perron (2001). To have tests with good properties they proposed selecting the lag order k using a Modified Information Criterion given by MIC(k) = ln(σ̂2ek ) + CT (τ T (k) + k)/T , where τ T (k) = (σ̂2ek )−1 b̂20 T −2 T 2 t=1 yt−1 . We have the MAIC when CT = 2 and the MBIC when CT = ln(T ). They recommend using the MAIC. Perron and Qu (2007) further refined the procedure by suggesting to select k using OLS-detrended data but constructing s2AR with GLS-detrended data once the order is selected. Many tests were also proposed to test the null hypothesis of stationarity versus the alternative 12 of a unit root (e.g., Kwiatkowski et al., 1992). But since the unit root tests described above are nearly optimal, they offer no added advantages to discriminate between I(0) and I(1) processes. The results of applying a battery of standard unit root tests indicate that the unit root hypothesis cannot be rejected for G, TRF, Temp_TSB and Temp_DS (Table 1). However, the existence of breaks in the trend functions of some series could lead to incorrectly classify them as DS processes. As discussed in the next sections, structural break tests and unit root tests that allow for a break in the trend function are required to investigate the data generating process of this type of series. 4 Unit root tests that allow for the presence of structural changes in the trend. The standard unit root tests for the existence of stochastic trends provided a first step to investigate the time-series properties of global and hemispheric temperatures and radiative forcing series. Nevertheless, it is important to note that they can be severely affected when the trend function is subject to changes in level and/or slope. Perron (1989) showed that the sum of the autoregressive coefficients is highly biased towards unity if there is a shift in the trend function. In this case, the unit root null hypothesis is hardly rejected even if the series is composed of white noise disturbances around the trend. Also, if the break occurs in the slope of the trend, standard unit root tests are not consistent, i.e., the null hypothesis of a unit root cannot be rejected even asymptotically. As a consequence, standard unit root tests are not adequate to investigate the data generating process of series which exhibit such behavior. Several tests have been devised to allow for structural breaks, we shall focus on the Kim-Perron test (2009) and briefly mention other methods that have been applied when investigating global and hemispheric temperatures. To motivate the problem addressed, it is useful to look at some basic properties of unit root and trend-stationary processes. As discussed earlier a DS process is such that the trend function changes every period (i.e., all shocks have a permanent effect on the future level of the series). On the other hand, for a TS process the trend never changes. Hence, one can view the unit root versus trend-stationary problem as addressing the following question: do the data support the view that the trend is changing every period or never? The question is then why restrict the comparison to ‘never’ or ‘always’ ? Would it not be preferable to make a comparison between ‘always’ and ‘sometimes’ ? Ideally, the proper question to ask would be ‘what is the frequency of permanent shocks?’. This is a question for which no satisfactory framework has been provided. The basic motivation for the work initiated by Perron (1989, 1990) is to take a stand on what is ‘sometimes’. The specific number chosen then becomes case-specific. His argument was that in many cases of interest, the relevant number of changes is relatively small, and often only one. These changes are then associated with important historical events. As far as statistical modelling 13 is concerned, the main conceptual issue is to view such changes as possibly stochastic but of a different nature than shocks that occur every period, i.e., drawn from a different distribution. However, since such large changes are infrequent, it is difficult to specify and estimate a probability distribution for them. The approach is then to model these infrequent large changes in the trend as structural changes. The question asked by unit root tests is then: ‘do the data favor a view that the trend is ‘always’ changing or is changing at most occasionally?’ or ‘if allowance is made for the possibility of some few large permanent changes in the trend function, is a unit root present in the stochastic component?’. Note two important qualifications. First, the setup allows but does not impose such large changes. Second, by “permanent” what should be understood is not that it will last forever but that, given a sample of data, the change is still in effect. Perron (1989) proposed an extension of the ADF test that allows for a one-time break in the trend function of a univariate time series. Three different model specifications were considered: the “crash” model that allows for a change in the level of the series; the “changing growth” model that permits a change in the rate of growth; and a third model that allows both changes. When analyzing global and hemispheric temperatures and radiative forcing variables the “changing growth” model is of interest, which can be briefly described as follows. In this case, the model is: yt = µ1 + β 1 t + (β 2 − β 1 ) DTt∗ + ut where DTt∗ = t − T1 if t > T1 and 0 otherwise. The null hypothesis is that ut is I(1) and the alternative hypothesis is that ut is I(0). The “changing growth” model takes an “additive outlier” approach in which the change is assumed to occur rapidly and the regression strategy consists in first detrending the series according the following regression: yt = µ + βt + γDTt∗ + yt (6) where γ = β 2 − β 1 and yt is accordingly defined as the detrended series. The test is based on the value of the t-statistic for testing that α = 1 in the following autoregression: yt = αyt−1 + k i=1 ai ∆yt−i + et (7) In the original Perron (1989) test the break is assumed to occur at a known date. Later, this test was generalized to allow the break to occur at an unknown date to be determined endogenously from the data (Zivot and Andrews, 1992; Perron 1997). The break date was originally proposed to be selected by 1) minimizing the t-statistic to test for a unit root; 2) minimizing/maximizing the t-statistic of the parameter associated with γ in regression (6) or; 3) maximizing the absolute value of the t-statistic for γ in regression (6). The resulting unit root test is then the t-statistic for 14 testing that α = 1 in regression (7) estimated by OLS using yt defined according to the estimate of the break date. Note that procedures (2) and (3) yield the same estimate of the break date and, hence, the same unit root test (it is also equivalent to select the estimate by minimizing the sum of squared residuals from regression (6)). However, (2) allows one to impose a prior on the sign of the change. This affects the limit distributions and leads to tests with higher power. The critical values of the tests have been tabulated in Perron (1997), see also Zivot and Andrews (1992). A problem with most procedures for testing for unit roots in the presence of a one-time break occurring at an unknown date is that the change in the trend is allowed only under the alternative hypothesis of a stationary noise component (Perron, 1997; Zivot and Andrews, 1992; Vogelsang and Perron, 1998). Consequently, it is possible that a rejection occurs when the noise is I(1) and there is a large change in the slope of the trend function. A method that avoids this problem is that of Kim and Perron (2009). Their procedure is based on a pre-test for a change in the trend function, namely the Perron and Yabu (2009a) test described below. If this pre-test rejects, the limit distribution of the unit root test is then the same as when the break date is known (Perron, 1989; Perron and Vogelsang, 1993; Kim and Perron, 2009). This is very advantageous since when a break is present the test then has much greater power. It was also shown in simulations to maintain good size in finite samples and that it offers improvements over other commonly used methods. As mentioned above, the presence of structural changes can have considerable implications when investigating time-series properties by means of unit root tests (Perron, 1989) and at the same time most tests for structural breaks require to correctly identify whether the data generating process is stationary or integrated. This creates a circular problem: depending on whether the process is I(0) or I(1) the limit distributions of these tests are different and if the process is misidentified the unit root tests will have poor properties. The Perron-Yabu test (Perron and Yabu, 2009a) was designed explicitly to address the problem of testing for structural changes in the trend function of a univariate time series without any prior knowledge as to whether the noise component is stationary or contains an autoregressive unit root. Their approach builds on previous work which analyzed the problem of hypothesis testing on the slope coefficient of a linear trend model when no information about the nature, I(0) or I(1), of the noise component is available (Perron and Yabu, 2009b). 4.1 Perron-Yabu testing procedure for structural changes in the trend function. We present the case of a model with a one-time structural break in the slope of the trend function with an autoregressive noise component of order one (AR(1)). A more detailed presentation of this case and of other structural change models and extensions can be found in Perron and Yabu 15 (2009a). Consider the following data generating process: yt = x′t Ψ + ut , ut = αut−1 + et , for t = 1, ..., T , et ∼ i.i.d. 0, σ2 , xt is a (r × 1) vector of deterministic components, and Ψ is a (r × 1) vector of unknown parameters which are model specific and described in the next paragraphs. The initial condition u0 is assumed to be bounded in probability. The autoregressive coefficient is such that |α| ≤ 1 and therefore, both integrated and stationary errors are allowed. The interest resides in testing the null hypothesis RΨ = γ where R is a (q × r) full rank matrix and γ is a (q × 1) vector, where q is the number of restrictions. The restrictions are used to test for the presence of a structural change in the trend function. For this purpose, Perron and Yabu (2009a) consider three models where a change in intercept and/or slope in the trend function occurs. In what follows, the break date is denoted T1 = [λT ] for some λ ∈ (0, 1), where [·] denotes the largest integer that is less than or equal to the argument and 1 (·) is the indicator function. The model to test for a one-time change in the slope of the trend function is specified with xt = (1, t, DTt∗ )′ and Ψ = (µ0 , β 0 , β 1 )′ where DTt∗ = 1 (t > T1 ) (t − T1 ) so that the trend function is joined at the time of the break. The hypothesis of interest is β 1 = 0. The testing procedure is based on a Quasi Feasible GLS approach with a superefficient estimate of α when α = 1. The estimate of α is the OLS estimate obtained from an autoregression applied to detrended data and is truncated to take a value 1 when the estimate is in a T −δ neighborhood of 1. This makes the estimate “super-efficient” when α = 1 and implies that in the case of a known break date, inference on the slope parameter can be performed using the standard Normal or Chi-square distribution whether α = 1 or |α| < 1. Theoretical arguments and simulation evidence show that δ = 1/2 is the appropriate choice. When the break date is unknown, the limit distribution is nearly the same in the I(0) and I(1) cases when considering the Exp functional of the Wald test across all permissible dates, see Andrews and Ploberger (1994). Hence, it is possible to have a test with nearly the same size in both cases. To improve the finite sample properties of the test, they also use a bias-corrected version of the OLS estimate of α as suggested by Roy and Fuller (2001). The testing procedure suggested is: 1) for any given break date, detrend the data by OLS to obtain the residuals ût ; 2) estimate an AR(1) model for ût yielding the estimate α̂; 3) use α̂ to get the Roy and Fuller (2001) biased corrected estimate α̂M ; 4) apply the truncation α̂MS = α̂M if |α̂M − 1| > T −1/2 and α̂MS = 1 if |α̂M − 1| ≤ T −1/2 ; 5) apply a Generalized Least Squares (GLS) procedure with α̂MS to obtain the estimates of the coefficients of the trend and the long-run variance of the residuals and construct the standard Wald-statistic WF MS ; 6) since the break date is assumed to be unknown, these 5 steps must be repeated for all permissible break dates to construct the Exp functional of the 16 Wald test denoted by Exp-WF S = log[T −1 Λ exp(WF MS (λ) /2)] where Λ = {λ; ǫ ≤ λ ≤ 1 − ǫ} for some ǫ > 0. A commonly used value is ǫ = 0.15. The results of applying the Perron-Yabu test to the series in Figure 1 confirm the existence of a break in the slope of the trend function for G, TRF and Temp_TSB, while no evidence for a break is found for AMO, Temp_DS and Temp_S (Table 2). The Perron-Yabu test correctly indicates that the breaks suggested by the visual inspection of Temp_DS and Temp_S are not part of their data generating process. The results in Tables 1 and 2 allow to conclude without further testing that Temp_DS has a unit root and that Temp_S and AMO are stationary, since when no break in trend is present the standard unit root tests reported in Table 1 are reliable. In contrast, the results strongly suggest that the analysis of the order of integration of G, TRF, Temp_TSB needs to be extended using unit root tests that allow for a break in their trend function. 4.2 Kim-Perron unit root tests with a one-time break in the trend function As mentioned, a problem with most tests for unit roots in the presence of a one-time break occurring at an unknown date is that the change in the trend function is allowed only under the alternative hypothesis of a stationary noise component. Hence, it is possible that a rejection occurs when the noise is I(1) and there is a large change in the slope of the trend. A procedure that avoids this problem is that of Kim and Perron (2009). It is based on a pre-test for a change in the trend function, namely the Perron and Yabu (2009a) test. If this pre-test rejects, the limit distribution of the unit root test is then the same as if the break date was known (Perron and Vogelsang, 1993). This is very advantageous since when a break is present the test has much greater power. The testing procedure under the additive outlier approach for the changing growth model consists in the following steps: 1) obtain an estimate of the break date T̂1 = λ̂T by minimizing the sum of squared residuals using regression (6). Then construct a window around that estimate defined by a lower bound Tl and an upper bound Th . A window of 6 observations is common. Note that, as shown by Kim and Perron (2009), the results are not sensitive to this choice; 2) create a new data set {ytn } by removing the data from to Tl + 1 to Th , and shifting down the data after the window by S (T ) = yTh − yTl ; hence, ytn = yt if t ≤ Tl and ytn = yt+th −tl − S (T ) if t > Tl ; 3) perform the unit root test using the break date Tl . This is the t-test statistic for testing that α = 1 in the following AO regression estimated by OLS, denoted by tα (λ̂tr ): ytn = αytn + k n i=1 ci ∆yt−i + et (8) where λ̂tr = Tl /Tr , Tr = T − (Th − Tl ) and ytn is the detrended value of ytn . The number of lags in (8) can be chosen by means of different information criteria for model selection such as the BIC and AIC. Extensions of this test are presented in Carrion-i-Silvestre et al. (2009). 17 Other unit root tests that allow for nonlinear trends that can accommodate structural changes have been applied to global and hemispheric temperatures. Becker, Enders and Lee (2006) proposed a stationarity test that allows for an unknown number of intercept shifts by using a selection of terms from a Fourier approximation. Harvey and Mills (2002) proposed a unit root test that allows for double smooth transitions to accommodate non-instantaneous breaks. Table 3 shows the results of applying the Kim-Perron test to the series in Figure 1 for which evidence of a break in the slope of their trend function was found (i.e., G, TRF, Temp_TSB). In all cases, once a better description of the deterministic trend is allowed, the unit root null is strongly rejected, reversing the results of the standard unit root tests in Table 1. 5 Perron-Zhu methodology for estimating a confidence interval for the break date Perron and Zhu (2005) analyzed the consistency, rate of convergence and limiting distributions of parameter estimates in models where the trend exhibits a slope change at some unknown date and the noise component can be either stationary or have an autoregressive unit root. Their results are of particular relevance when considering the problem of selecting the break date when testing for structural changes and for deriving the limiting distributions of unit root tests that allow for a one-time structural change that occurs at an unknown date, such as in the Kim and Perron (2009) test. Another important practical application of deriving the limiting distribution of the estimate of the break date is that it permits forming a confidence interval for the break date. Perron and Zhu (2005) considered a total of six models with deterministic and stochastic trends. The random component was assumed to be either stationary or to contain a unit root, while for the deterministic component three cases were considered: 1) a first-order linear trend with a one-time change in the slope such that the trend function is joined at the time of the break; 2) a local disjoint broken trend and; 3) a global disjoint broken trend. The case of interest here is the first specification with a stationary noise component. The interested reader is referred to Perron and Zhu (2005) for the specifications and limiting distributions for the other models considered. The deterministic part is specified as τ t = µ + βt + γDTt∗ , where DTt∗ = t − T1 if t > T1 and 0 otherwise with T1 = λT the break date and λ the break fraction. Note that at the time of the break, the slope coefficient changes from β to β + γ but that the trend function is continuous at T1 . This specification is therefore referred to as the “joint broken trend”. The estimation method is simply to select the break date that minimizes the sum of squared residuals from a regression of the series of interest yt (t = 1, ..., T ) on the regressors {1, t, DTt∗ }, i.e., applying OLS to the model yt = µ + βt + γDTt∗ + ut 18 (9) Denote the estimate by T̂1 and the associated estimate of the break fraction by λ̂ = T̂1 /T . They 2 showed the limit distribution to be T 3/2 (λ̂ − λ) →d N(0, 4σ2 /[λ0 (1 − λ0 ) γ 0 ]) where γ 0 is the true value of the change in the slope parameter and σ2 = limT →∞ E( T 2 t=1 ut ) is the so-called long-run variance of ut . Hence, it is easy to construct a confidence interval using the estimates of σ2 , γ 0 and λ0 . A common estimate of σ2 is based on a weighted sum of the autocovariance function of ut of the form Ŝw,T = R̂u (0) + 2 T −1 j=1 w(j, m)R̂u (j), where R̂u (j) = T −1 T t=j+1 ût ût−j with ût the OLS residuals from regression (9). Here, w(j, m) is some weight function. A popular choice is the Bartlett triangular weight with w(j, m) = 1 − j/(m + 1) if j ≤ m and 0 otherwise. Other choices are available such as the Parzen or Quadratic Spectral. The parameter m is a bandwidth or truncation parameter. A popular method to select this parameter is due to Andrews (1991). To illustrate the use of this procedure with observed climate variables, the 95% confidence intervals for the estimated break dates of G and TRF were constructed. The confidence interval for G is (1946, 1966) and for TRF is (1956, 1964), which overlap considerably suggesting that the break in the slope of their trend functions is a common feature shared by both variables. The long-run variance was estimated using the Bartlett kernel with Andrews’ bandwidth selection method. 6 Cointegration A set of integrated variables is said to be cointegrated if there exists a linear combination of them that is stationary (Engle and Granger, 1987). Cointegration implies that there is a long-run equilibrium relationship between two or more variables because they share a common stochastic trend (Stock and Watson, 1988). To illustrate this concept consider the following integrated processes: xt = µxt + εxt and zt = µzt + εzt where µxt and µzt are unit root processes which contain a stochastic trend and εxt , εzt are stationary noise processes. For these two variables to be cointegrated there must exist a linear combination α1 xt + α2 zt (with α1 and α2 different from zero) such that µxt − (α2 /α1 )µzt is stationary, indicating that the stochastic trends µxt and µzt are identical up to a scalar. Consequently, xt and zt share a common secular movement and a regression between these two variables will produce stationary residuals. In general, the n variables yt = (y1t , ..., ynt ) are said to be cointegrated of order d, b, or yt ∼ CI(d, b), if all of the variables are integrated of order d and there exists a cointegrating vector β such that the linear combination β ′ yt is integrated of order (d − b), with b > 0. Granger and Lee (1990) showed that multicointegration can also occur where an equilibrium relationships among sets of variables that have different orders of integration can occur. For example, if two of the variables in yt are I(2) and the others are I(1) an equilibrium relationship can hold if the two I(2) variables are CI(2, 1) and this linear combination cointegrates with the other I(1) variables. 19 The tests that have been devised to analyze the existence of cointegration relationships can be classified in three types: single equation residual based, single equation error correction and multiple equations based. These three types of tests have been used to investigate the existence of a common secular stochastic trend in global and hemispheric temperatures and radiative forcing data. The most widely used residual based test is the Engle and Granger (1987) two-steps procedure. Consider the OLS regression partitioning yt into y1t and yt∗ = (y2t , ..., ynt ): y1t = c + β ′ yt∗ + et (10) where (1, −β ′ ) is the cointegration vector normalized to have value of 1 for the variable y1t , which assumes that y1t is part of the set of variables that are cointegrated. Note also that usually only a constant is included, which implies that the cointegrating vector that eliminates the stochastic non-stationarity also eliminates the non-stationarity due to the deterministic trends in the variables (see Campbell and Perron, 1991). The second step is to test if et is I(0) by means of a unit root test such as the ADF or the Phillips-Perron tests applied to the OLS residuals from (10). The critical values of the cointegration test follow a non-standard distribution and are different from those of the original unit root tests. The critical values depend on the number of I(1) regressors included in (10) and on whether a constant and/or a linear trend is included and, if only a constant is included, on whether the data are trending. Note also that this test is valid only if the vector of variables yt∗ is not cointegrated (see, Phillips and Ouliaris, 1989, and Hansen, 1992). To illustrate other tests and features of cointegrated systems, it is useful to adopt a Vector Autoregressive (VAR) model yt = µ + p i=1 Ai yt−i + et , △yt = c + Πyt−1 + where Π = −(I − p i=1 Ai ) and Γi = − with et ∼ i.i.d. (0, Σ). It can be written as p−1 i=1 Γi p j=i+1 Aj . △ yt−i + et (11) The set of variables yt is cointegrated if rank(Π) = r for some 0 < r < n. We then have the decomposition Π = αβ ′ where α and β are n × r matrices and β is the matrix of cointegrating vectors. We can then write (11) as △yt = c + αβ ′ yt−1 + p−1 i=1 Γi △ yt−i + et = c + αzt−1 + p−1 i=1 Γi △ yt−i + et (12) where zt−1 = β ′ yt−1 is the disequilibrium error at time t − 1. The representation (12) is called an error correction model (ECM) and the matrix α consists of the vectors of adjustment coefficients. For simplicity, suppose that there is a single cointegrating vector so that r = 1 and α and β are vectors of dimension n. Now, suppose that α = (α1 , 0, ..., 0). Then, yt∗ = (y2t , ..., ynt ) is said to be weakly exogenous. While this appears to be a special case, it is highly relevant in the context of climate change. Suppose we have two series with y1t temperatures and y2t some forcing variable 20 (e.g., well-mixed greenhouse gases). Then, from standard physics of climate change, y2t is weakly exogenous. A test for cointegration can then be based on the first equation of the ECM (12): ∗ ∆y1t = c + α′1 ∆yt∗ + γ 0 ẑt−1 + γ ′1 yt−1 + k ′ ∗ i=1 (π i ∆yt−i + φi ∆y1t−i ) + ǫt . (13) ′ where ẑt−1 = β̂ yt−1 with β̂ the estimate of the cointegrating vector. Using the representation (13), there is cointegration if γ 0 = 0. Hence, a test of no-cointegration amounts to testing H0 : γ 0 = 0, using a standard t-statistic; see, Kremers et al. (1992), Banerjee et al. (1993), Banerjee et al. (1986), Banerjee et al. (1998) and Ericsson and Mackinnon (2002). It is also possible to do a joint test H0 : (γ 0 , γ ′1 )′ = 0 using a Wald test (see Boswijk, 1994). These tests were shown to have much higher power when the condition of weakly exogenous regressors is satisfied. The most popular system-based tests for cointegration were proposed by Johansen (1988) and Johansen and Juselius (1990) and are based on the VAR model (11). To describe them as well as the method to obtain the Maximum Likelihood Estimates (MLE) of the cointegrating vectors, define the following: R0t (resp., R1t ) as the residuals from a regression of △yt (resp., yt−1 ) on a constant T ′ t=1 Rit Rjt . Then the model (11) can be written as R0t = ΠR1t +et . known, then R0t = α(β ′ R1t ) + et and the MLE (assuming normal errors) of α on β are α̂(β) = S01 β(β ′ S00 β)−1 and Σ̂(β) = S00 − S01 β(β ′ S11 β)−1 β ′ S01 . The and the lags of △yt ; Sij = T −1 Suppose that β is and Σ conditional MLE of β is then obtained by minimizing |Σ̂(β)|. This is achieved by solving the eigenvalue system −1 S01 = 0 λS11 − S10 S00 (14) subject to the restrictions β ′ S11 β = Ir . Note that these impose r(r+1)/2 restrictions on β, which are needed since only the space spanned by the cointegrating vectors is identified. The system (14) has n eigenvalues, represented in decreasing orders as λ̂1 λ̂2 ... λ̂r λ̂r+1 ... λ̂n . The estimates of the cointegrating vectors are then the eigenvectors corresponding to the r largest eigenvalues. For testing, note that if the rank of Π is r then the population values of the n − r smallest eigenvalues are zero. Hence, this leads to two testing procedures. The first is H0 : rank(Π) = r versus H1 : rank(Π) > r and the test statistic, called the trace test, is T n i=r+1 ln(1 − λ̂i ). The second test had a different alternative hypothesis so that H0 : rank(Π) = r and H1 : rank(Π) = r + 1, and the test statistic, called the maximum eigenvalue test, is T ln(1 − λ̂r+1 ). This test is especially useful to devise a strategy to estimate the number of cointegrating vectors r. Both tests have non-standard distributions, which however depend only on n and r. In the application of cointegration techniques to climate variables the common stochastic trend has been interpreted as the fingerprint of anthropogenic activities in global and hemispheric temperatures. It should be noted, however, that a requirement for this concept is the existence of stochastic trends with data having a linear deterministic time trend. 21 Figure 2(a) illustrates the concept of cointegrated series. DS_1 and DS_2 are simulated DS series that share a common stochastic trend which determines their secular movement, and although DS_1 and DS_2 can show deviations from this trend, these are only temporary. These were generated by DS_1t = T emp_DSt +e1t with e1t ∼ i.i.d. N (0, 0.09) and DS_2t = T emp_DSt +e2t with e2t ∼ i.i.d. N (0, 0.0025) where T emp_DSt is as defined in Section 2.7. In contrast, Figure 2(b) shows two simulated DS series with independent stochastic trends (DS_3 and DS_4). These were generated by DS_3t = 0.01 + DS_3t−1 + e3t with e3t ∼ i.i.d. N (0, 0.01) and DS_4t = 0.01+DS_4t−1 +e4t with e4t ∼ i.i.d. N (0, 0.01) (e1t through e4t are mutually independent). These series deviate persistently from each other, precluding the existence of any linear combination that could be stationary. Two of the most widely used cointegration tests were applied to DS_1 and DS_2 and DS_3 and DS_4. The results (Table 3) correctly reflect the true data generating process used to construct the simulations: both the Johansen (trace and max-eigenvalue statistics) and the Engle-Granger tests indicate the existence of a cointegrating relationship between DS_1 and DS_2, while clearly rejecting it for DS_3 and DS_4. 7 Co-trending We now discuss various procedures to test for the presence of a common non-linear trend. We start with Bierens’ nonparametric nonlinear co-trending test. Nonlinear co-trending is special case of the more general “common features” concept described by Engle and Kozicki (1993). As mentioned, cointegration analysis is only one possibility for investigating secular co-movements of non-stationary variables. A long-run equilibrium relationship can be established when linear combinations of different time series cancel out “common features” such as trends and co-breaks in trend. Nonlinear cotrending aims to formally establish the existence of a secular co-movement between time series that are stationary processes around nonlinear deterministic trends. The advantage of the test proposed by Bierens (2000) is that the nonlinear trend does not have to be parameterized. The nonlinear trend stationarity model considered can be expressed by zt = g (t) + ut with g (t) = β 0 + β 1 t + f (t), where zt is a k-variate time series, ut is a kvariate zero-mean stationary process and f (t) is a deterministic k-variate general nonlinear trend function that allow, in particular, structural changes. Nonlinear co-trending occurs when there exists a non-zero vector θ such that θ′ f (t) = 0. Hence, the null hypothesis is that the multivariate time series zt is nonlinear co-trending, implying that there is one or more linear combinations of the time series that are stationary around a constant or a linear trend. Note that this test is a cointegration test in the case when it is applied to series that contain unit roots. The nonparametric test for nonlinear co-trending is based on the eigenvalues of the matrices M1 and M2 defined by 22 [T x] T ′ −1 t=1 (zt − β̂ 0 − β̂ 1 t) if x t=1 F̂ (t/T )F̂ ((t/T )) where F̂ (x) = T −1 0, T with β̂ 0 and β̂ 1 being the estimates of the vectors of intercepts M1 = T −1 ∈ T −1 , 1 , F̂ (x) = 0 if x ∈ and slope parameters in a regression of zt on a constant and a time trend; also M2 = T −1 T −1 t=m [m m−1 j=0 (zt−j − β̂ 0 − β̂ 1 (t − j))][m−1 m−1 j=0 (zt−j − β̂ 0 − β̂ 1 (t − j))]′ where m = T α with T is the number of observations and α = 0.5 as suggested by Bierens (2000). Solving |M̂1 − λM̂2 | = 0 and denoting the solution λ̂r , the test statistic is T 1−α λ̂r . The null hypothesis is r co-trending vectors against the alternative of r − 1 co-trending vectors. This test has a non-standard distribution and the critical values have been tabulated by Bierens (2000). The existence of r co-trending vectors in r + 1 series indicates the presence of r linear combinations of the series that are stationary and that these series share a single common nonlinear deterministic trend. Such a result indicate a strong secular co-movement in the r + 1 series. Another procedure used to investigate long term co-movements of hemispheric temperatures is that of co-shifting (Holt and Teräsvirta, 2012), based on the Shifting Mean Vector Autoregression, a generalization of the Shifting Mean Autoregressive model of González and Teräsvirta (2008). A simple approach can also be used to test for a common long-run path in the bivariate context. Assume two trend stationary variables of the form: yt = dt + ut and xt = wt + vt , where dt and wt are nonstationary components which may be composed of a wide class of linear and nonlinear trends, with changes in slopes or intercepts, and ut and vt are stationary noise components. The procedure is based on testing for remaining nonstationarities in the residuals of the following regression estimated by OLS: yt = α + βxt + εt . If the individual nonstationary components dt and wt are not present in εt , that is the residuals are found to be stationary, then it is said that yt and xt share the same long-run path. The existence of remaining non-common nonstationarities can be evaluated by applying standard unit root tests. This procedure is similar to the Engle and Granger (1987) two-steps cointegration test but it does not require the assumption of unit roots in radiative forcing and temperature variables. Also, since all series are trend-stationary, the relevant critical values are those tabulated for standard unit root tests with no deterministic terms included. G and TRF can be used to illustrate the concepts of co-trending and common long-run path described above. The results the Bierens nonparametric nonlinear co-trending test indicate the existence of one co-trending vector and the existence of a common secular trend between G and TRF, as expected by climate physics (Table 4). This finding is also supported by the results of the simple common long-run path test based on applying an ADF test to the residuals of a regression of G on TRF (test statistic value equal to -7.36, significant at the 1% level). 23 8 Conclusions As evidenced in this review, the time-series based study of the detection and attribution of climate change has been an active area of research and the application of econometric methods had a prominent role. With few exceptions, recent publications have shown compelling evidence based on both observations and climate models’ simulations supporting the existence of a warming trend. Furthermore, although strong methodological differences have characterized these studies and a wide variety of tests and models have been applied, the results are remarkably in agreement regarding the dominant role of human intervention with the climate system. Whatever the evolution of climate change during this century turns out to be, the time-series and econometric approach will be of importance in understanding the changes in the climate system as they occur as well as their linkage to human and natural drivers. Global warming will continue to provide a fertile field for the development and application of new statistical models, tests and methods. 24