Significant advances have been made in documenting how global and hemispheric temperatures have evolved and in learning about the causes of these changes. On the one hand, large efforts have been devoted to investigate the time series properties of temperature and radiative forcing variables (Gay-Garcia et al., 2009;Kaufmann et al., 2006;Mills, 2013;Tol and de Vos, 1993). In addition, a variety of methods were applied to detect and model the trends in climate variables, including features such as breaks and nonlinearities (Estrada, Perron and Martínez-López, 2013;Gallagher et al., 2013;Harvey and Mills, 2002;Karl et al., 2000;Pretis et al., 2015;Reeves et al., 2007;Seidel and Lanzante, 2004;Stocker et al., 2013;Tomé and Miranda, 2004). Multivariate models of temperature and radiative forcing series provide strong evidence for a common secular trend between these variables, and help to evaluate the relative importance of its natural and anthropogenic drivers (Estrada, Perron and Martínez-López, 2013;Estrada, Perron, Gay-García and Martínez-López, 2013;Kaufmann et al., 2006;Tol and Vos, 1998). The methodological contributions of the econometrics literature to this field have been notable; e.g., Dickey and Fuller (1979), Engle and Granger (1987), Johansen (1991), Perron (1989Perron ( , 1997)), Bierens (2000), Ng and Perron (2001), Kim and Perron (2009), Perron and Yabu (2009), among many others, see Estrada and Perron (2014) for a review. Regardless of the differences in assumptions and methods (statisticalor physical), there is a general consensus about the existence of a common secular trend between temperatures and radiative forcing variables. Some of the most relevant questions about the attribution of climate change are concerned with understanding particular periods in which the rate of warming changed; e.g., rapid warming, slowdowns or pauses. An example is the apparent slowdown in warming since the 1990s, for which various methods were applied to document or reject its existence (e.g., Fyfe et al., 2016;Lewandowsky et al., 2015Lewandowsky et al., , 2016)). These issues are delicate to handle when dealing with observed temperature series, as the object of interest is the secular trend behind the observed warming. This underlying trend is affected by low-frequency natural variability (Swanson et al., 2009;Wu et al., 2011) and changes in its rate of warming are difficult to detect and attribute. New tests and approaches to investigate common features in temperature and radiative forcing can make attribution studies more relevant for climate science and policy-making by providing a better understanding of the drivers behind them. One way to tackle this problem is to devise procedures to extract the common secular trend between temperature and radiative forcing series and framing the problem in a univari-ate context where the available structural change tests can be applied. Estrada and Perron (2016) used this approach to investigate the existence and causes of the current slowdown in the warming. Their results strongly suggests a common secular trend and common breaks, largely determined by the anthropogenic radiative forcing. Their analysis is based on the results of co-trending and principal component analyses to separate the common long-term trend imparted by radiative forcing from the natural variability component in global temperature series. As they discussed, filtering the effects of physical modes of natural variability from temperature series is necessary to obtain a proper assessment of the features and drivers of the warming trend. This problem has seldom been addressed within the time-series based attribution literature (e.g., Estrada, Perron and Martínez-López, 2013;Estrada and Perron, 2016) and it constitutes a relevant research topic that requires the development of new procedures and techniques. What transpires from this research is that temperatures and radiative forcing are most likely characterized by a linear trend with two changes in the rate of growth. The first occurs in the early 60s and indicates a very large increase in the rate of growth of both temperatures and radiative forcing. This was termed as the "onset of sustained global warming". The second is related to the more recent so-called "hiatus" period, which suggests that temperatures and total radiative forcing have increased less rapidly since the mid-90s compared to the larger rate of increase from 1960 to 1990. The approach we take consists in designing statistical tests in a multivariate setting with joint-segmented trends to investigate the existence of common breaks and the "hiatus". There are two issues that until now remain unresolved, which this approach can address. The first is whether the breaks in the slope of the trend functions of temperatures and radiative forcing are common. This is important because common breaks coupled with the basic science of climate change would strongly suggest a causal effect from anthropogenic factors to temperatures. As is well known, a common linear trend can occur when the series are spuriously correlated. A common break can eliminate such concerns of spurious correlation and foster the claim that the relationship from radiative forcing to temperatures is causal. Hence, our paper first aims at developing a test for common breaks across a set of series modeled as joint segmented trends with correlated noise. The theoretical framework follows Qu and Perron (2007) and Oka and Perron (2016), building on Bai andPerron (1998, 2003). Their framework, however, precludes joint segmented trends since the regressors are functions of the break dates, which makes the problem very different. Perron and Zhu (2005) show that the limit distribution of the estimate of the break date in a single time series with a joint segmented trend follows a normal limit distribution. We build on that result to show that our common break test follows a chi-square distribution. Theoretically, our results are more general and cover a wide range of cases to test hypotheses on the break dates in a multivariate system, within or across equations. Although our model with a single equation is similar to that in Perron an Zhu (2005), our work has a different focus, namely testing, while Perron and Zhu (2005) deals with the distribution theory of the estimate. We are not aware of any test in the literature that can handle a set of general linear restrictions on break dates both across and within equations.
Furthermore, our test for an additional break is more closely related to Perron and Yabu (2009), and the studies cited therein. The attempt to increase the power of structural break tests using cross equation correlation is also new in the case of joint-segmented trends. Our results show that, once we filter temperatures for the effect of the Atlantic Mutidecadal Oscillation (AMO) and the North Atlantic Oscillation (NAO) for reasons explained in the text, the breaks in the slope of radiative forcing and temperatures are common, both for the large increase in the 60s and the "hiatus" period. The second issue relates to establishing formally via a proper testing procedure that takes into account the noise in the series, whether there was a "hiatus" period for temperatures since the mid-90s. This is important because such a test would counter the widely held view that the "hiatus" is the product of natural internal variability (Kosaka and Xie, 2013;Trenberth and Fasullo, 2013;Meehl et al., 2011;Balmaseda, Trenberth and Källén, 2013). Using standard univariate tests (e.g., Perron and Yabu, 2009), the results are mixed across various series and sometimes borderline. Our aim is to provide tests with enhanced power by casting the testing problem in a bivariate framework involving temperatures and radiative forcing. Our results indicate that indeed the "hiatus" represents a significant slowdown in the rate of increase in temperatures, especially when considering global or southern hemisphere series, for which our test points to a rejection of the null of no change for all data sources. We consider a multivariate system with n equations where the dependent variables are modeled as joint-segmented trends with multiple changes in the slope. The errors are allowed to be serially correlated and correlated across equations. Of interest is testing for general linear restrictions on the break dates including testing for common breaks across equations. The test used is a (quasi-) likelihood ratio test assuming serially uncorrelated errors, in which case it has a pivotal chi-square distribution. However, it is non-pivotal in the general case of interest. Accordingly, we also consider a corrected Wald test with a pivotal limit distribution, which can be constructed using break dates estimated one equation at a time, labelled OLS-Wald, or estimated via the complete system, labelled GLS-Wald. The latter can offer more efficient estimates when there is correlation in the errors across equations. It is, however, computationally demanding as it requires least squares operations of order O(T m ) where T is the sample size and m is the total number of breaks. The OLS-Wald test requires only least-squares operations of order O(T m 1 + ... + T mn ) where m i (i = 1, ..., n) is the number of breaks in the i th equation, making it much easier to compute in a multivariate system. Our simulation study shows that the two Wald tests have finite sample sizes close to the nominal size when the extent of the serial correlation in the noise is small. They, however, suffer from potentially severe liberal size distortions for moderate to strong serial correlation. Hence, for all three tests (since the LR is non-pivotal), we suggest a bootstrap procedure to obtain the relevant critical values. These bootstrap tests are shown to have exact sizes close to the nominal 5% level in all cases considered. We also propose a test for the presence of an additional break in some series. For instance, in our applications, the null and alternative hypotheses specify a break in the slope of the radiative forcing series (for which a rejection is easy using a single equation test), while under the null hypothesis there is no break in the temperature series but there is one under the alternative. Hence, the aim is to see whether a break in temperatures is present using a bivariate system. The theoretical framework is more general and allows multiple breaks in a general multivariate system. The test considered is a quasi-likelihood ratio test, whose limit distribution is non-pivotal and depends in a complex way on a number of nuisance parameters. Hence, we resort to a bootstrap procedure to obtain the relevant critical values. The paper is structured as follows. Section 2 presents the statistical model and the assumptions imposed. Section 3 considers testing for linear restrictions on the break dates including as a special case testing for common breaks across equations. The various tests and their limit distributions are discussed. Section 4 presents the test for the presence of an additional break in some series within a multivariate system. Section 5 presents simulation results about the size and power of the tests and introduces the bootstrap versions recommended. Section 6 presents the applications. The data are discussed in Section 6.1. Section 6.2 presents the results for the common breaks tests, while Section 6.3 addresses the issue of testing for the "hiatus" using a bivariate system with a temperature and a radiative forcing series. Brief conclusions are offered in Section 7. All proofs are in an appendix. We consider the following model with multiple breaks in a system of n variables. Each variable consists of a linear trend with multiple changes in slope with the trend function is joined at each break date. More specifically, there are m i breaks in the slope of the i th variable, and the n-variate system is:

for t = 1, ..., T and i = 1, ..., n, where b t (k

is the indicator function of the event A and k 0 ij is the j th break date for the change (with magnitudes δ ij ) in the trend of the i th variable. The vector of break dates for the i th variable is k 0 i = (k 0 i1 , . . . , k 0 im i ) ′ and k 0 is the vector of break dates for the entire system, k 0 = (k 0′ 1 , . . . , k 0′ n ) ′ .
We let m = n i=1 m i denote the total number of breaks in the system. Hence, k 0 is an m × 1 vector. For the i th variable, we have in matrix form

where

The entire system in matrix notation is

where

We are interested in the case where the number of breaks for each variable is known but the break dates are unknown. For a vector of generic break dates, we will use the notation without a 0 superscript, that is, k and k i . Also, we will simply write X to denote the regressor matrix corresponding to break dates specified by k, that is, X = diag(X(k 1 ), . . . , X(k n )). To motivate our quasi-likelihood ratio test, we first assume that u is multivariate normal with 0 mean and covariance Σ ⊗ I T . This assumption will be relaxed when we derive the asymptotic distribution of our test. For a generic break date vector k and the corresponding regressor matrix X, the log-likelihood function of y is then given by

Let θ = ( θ′ 1 , . . . , θ′ n ) ′ and Σ(k) be the maximum likelihood estimators for θ and Σ, which jointly

Thus, the maximized the log-likelihood function is

It is useful to work with break fractions; hence define λ ij = k ij /T , λ i = k i /T and λ = k/T , with the true break fractions λ 0 ij , λ 0 i and λ 0 equivalently defined. With these definitions, functions of λ such as l(λ) and Σ(λ) will be used interchangeably with l(k) and Σ(k). We now state the assumptions needed for the asymptotic analysis. Assumption 2 δ ij = 0 for j = 1, ..., m i and i = 1, . . . , n.

Assumption 3 Let u t = (u 1t , . . . , u nt ) ′ . Then, u t is stationary with E(u t ) = 0 and V ar(u t ) = Σ. In addition,

] t=1 e t +o p (1) and T -1/2 [T r] t=1 e t ⇒ W (r), where "⇒" denotes weak convergence under the Skorohod topology, W (r) is the n dimensional standard Wiener process and

Assumptions 1 and 2 are standard and simply state that the break dates are asymptotically distinct (i.e., each regime increases proportionally with the sample size T ) and the changes in the parameters are non-zero at the break dates. Assumption 3 states that u t is a stationary process and its partial sum follows a functional central limit theorem. 3 Testing for Linear Restrictions on the Break Dates

We now consider the following null and alternative hypotheses:

Some examples of this type of hypotheses are the following: 1) specific break fractions: R = I m and r = λ; 2) break fractions with a fixed distance: R = [0, . . . , 1, 0, . . . , -1, 0, . . . ] and r = c > 0; 3) common breaks: R = [0, . . . , 1, 0, . . . , -1, 0, . . . ] and r = 0. For the first case, historic events can be used to determine the value of λ. The second case can be used across equations as well as within an equation. The third case can be viewed as a special case of the second one and should be used across equations only. The LR test statistic is

Theorem 1 Let Assumptions 1-3 hold with Σ = Ψ. Then, as T → ∞, we have under H 0 in (5) that LR d → χ 2 q , where q = rank(R). Our model with broken segmented trends at unknown dates is a rare case in structural change problems with the estimates of the break dates having an asymptotic normal distribution; see Perron and Zhu (2005). Since the likelihood ratio test essentially involves quadratic forms of these estimates, it has the usual χ 2 distribution. When the short-run and long-run variances differ Σ = Ψ, the asymptotic distribution depends on nuisance parameters. It is difficult to directly modify the LR test. Instead, we develop a Wald test for the hypotheses (5) robust to serial correlation in the errors. To express the limiting distribution, define

′ dr and GG ij = g i (r)g j (r) ′ dr, with all integrals from 0 to 1. Denote by s ij the (i, j) th element of Σ -1 and by κ ij that of Σ -1 ΨΣ -1 . Now, define

We also define

GF and Γ GG equivalently, as well as

where

provide the limiting distribution of the estimate of the break fraction that maximizes the log-likelihood (4), which allows us to obtain the limit distribution of the Wald test. Theorem 2 Let Assumptions 1-3 hold, λ = k/T and k = arg min k l(k). Then, as T → ∞,

q , where q = rank(R) and

The total number of regressions required to obtain the above break date estimator λ is O(T m ). This can pose a considerable amount of computational burden especially if the procedure is to be bootstrapped. Hence, it is worthwhile to devise a test assuming diagonality in Σ because the break dates are then estimated separately from each variable and the total number of regressions required is only

While this modification lessens computational cost significantly, the resulting break date estimators can be less efficient due to the neglected correlation between variables in the system. Suppose the break dates are estimated separately for each variable from (2). Thus,

where

The limiting distribution of λ is obtained as a special case of Theorem 2 by letting Σ = I. The limiting covariance Ξ can be simplified somewhat. Define

where

q where q = rank(R) and Ξs → p Ξ s .
The result in (i) coincides with one result reported in Perron and Zhu (2005) when there is only one equation with one break. We consider a test for the null of m = n i=1 m i breaks in the system against the alternative of m + 1 breaks with the location of the additional break unspecified. We again consider a quasi-likelihood framework. Suppose that the model is now given by

where a i (h) := 1 i ⊗ b(h) for some h = k i1 , . . . , k im i and 1 i is an n × 1 vector with 1 in the i th element and 0 elsewhere. The null and alternative hypotheses are now:

The maximum of the likelihood function under H 0 is again given by ( 4) for a generic break date vector k. Similarly, the maximum of the likelihood function under H 1 is:

for a generic break date vector k and an additional break date h. The notation l (i) (k, h) and Σ(i) (k, h) is used to indicate the fact that an additional break is inserted in the i th equation. This additional break can only occur in one of the equations and we assume for simplicity that the relevant equation is known. This assumption is relaxed below. Just like we have defined λ = k/T and used l(k) and l(λ) interchangeably, we define ν = h/T and will use l (i) (k, h) and l (i) (λ, ν) interchangeably. The test statistic we consider is given by

where λ = arg sup λ l(λ) and

As is well known in the structural break literature, γ is not identified under the null hypothesis and ν must be restricted to ensure a non-divergent limiting distribution of the test statistic. The set

T imposes the relevant restrictions.1 When the equation with an additional break is not specified a priori, the LR statistic can be extended to be

In order to express the limiting distribution, we need to define additional terms. Let f i (r) and g i (r) be as defined before. Let b(r, ν) = (r -ν) + , F B i (ν) = f i (r)b(r, ν)dr, GB i (ν) = g i (r)b(r, ν)dr and BB(ν) = b 2 (r, ν)dr where all integrals are taken from 0 to 1. Recall that s ij is the (i, j) th element of Σ -1 and κ ij is that of Σ -1 ΨΣ -1 , and let

. Now, using these functions, we define

where D δ , Q GF , Q F F , Γ GF and Γ F F are as previously defined. Finally, let η (i) (ν) be a mean zero Gaussian process defined over the unit interval with

Also, denote the limit counterpart of the set

T by C (i) , replacing λij by λ 0 ij . The limit distribution of the tests are stated in the following theorem. Theorem 3 Suppose that Assumptions 1∼3 hold. Then, as T → ∞, we have under H 0 in ( 7)

The limiting distributions depend on various nuisance parameters. First, all of the true break fractions λ 0 ij for j = 1, . . . , m i and i = 1, . . . , n matter via terms such as

Second, both the short-run and longrun variances matter via the various Q and Γ matrices. Third, the trimming parameter ε ps in the set C (i) matters. Given the complexity of the limiting distributions, we neither seek a way to eliminate nuisance parameters nor attempt to tabulate critical values. Instead, we resort to using a bootstrap method to generate p-values. The various terms are evaluated in a manner similar to that of the tests for common breaks described in the next section. Note that if the short-run and long-run variances coincide and the true break fractions λ 0 ij are used instead of their estimates, V ar(η (i) (ν)) reduces to be ξ

In this section we provide Monte Carlo simulation results to assess the adequacy of the asymptotic distributions derived in Sections 3 and 4. All results are based on 1,000 replications. We focus on the test for common breaks as the tests for an additional change in slope are too computationally demanding so that a reasonable simulation experiment is prohibitive. The data generating process (DGP) is specified by:

The results are exactly invariant to the values of µ i and β i . For the slope change parameters δ 11 and δ 21 , the cases of 0.5, 1.0 and 1.5 are considered. The error terms are such that

and (e 1t , e 2t ) ′ = α(e 1t-1 , e 2t-1

The parameter ρ stands for the correlation across equations and α stands for the autoregressive parameter in each of the error terms. We consider (-0.5, 0, 0.5) for the values of ρ and (0, 0.3, 0.7) for the values of α. What influences the precision of the break date estimate is the ratio between the break magnitude and the long-run variance. In our simulations, the long-run variance is always set to be unity. The sample size is set to T = 100 and the true break dates are at mid-sample, k 0 11 = k 0 21 = 50. For each generated data, we test four null hypotheses: (1) R = I 2 and r = (0.5, 0.5) ′ , (2) R = I 2 and r = (0.525, 0.475) ′ , (3) R = [1, -1] and r = 0 and (4) R = [1, -1] and r = 0.05. Hence the null rejection probabilities for (1) and (3) pertain to the finite sample sizes (with a nominal size 5%) while those for (2) and (4) pertain to powers.
In testing these hypotheses, we consider the LR and Wald tests based on the SUR system (Theorems 1 and 2) as well as the Wald test based on estimating the break dates equation by equation using OLS regressions (Corollary 1). We refer to the two Wald tests as the GLS-Wald and OLS-Wald, respectively. For both versions, the long-run variance is estimated using a quadratic spectral window, where the bandwidth parameter is selected using Andrews' (1991) data dependent method with an AR(1) approximation. We also simulate the bootstrap version of the aforementioned tests. To obtain bootstrap p-values, we first estimate break dates imposing the null hypothesis and remove the estimated trend. Then, we fit a VAR(1) on the resulting residuals. Following Kilian (1998), we compute the bias-corrected estimates of the VAR coefficients and obtain the corresponding residuals of the VAR. From them, we construct a pseudo VAR process and add to it the estimated trend functions to obtain a bootstrap sample. Simulating bootstrap tests requires a lot of computing time especially for those based on the full system. We make use of the Warp-Speed method suggested by Giacomini, Politis and White (2013). Also, we use the feasible GLS estimate instead of the MLE in the construction of the LR and GLS-Wald tests. In Table 1, the results for the LR test are reported. Regarding the finite sample size obtained using the asymptotic critical values, two features emerge. First, the size is near or below the nominal level when α = 0, but it climbs up quickly as α increases. In the worst case, it can go up to 50%. This result is well expected since the asymptotic validity of the LR test holds only when α = 0. Second, the size inflation is more evident when δ i is small. In fact, when δ i is very large, the estimated break dates coincide with the true ones making the test statistic literally zero with a large probability. However, this type of conservativeness is of little concern because the power function does not appear to be decreasing as δ i gets large. On the other hand, the finite sample size obtained by the bootstrap procedure stays near the nominal 5% level across all simulation designs, offering significant improvements over the asymptotic test. The bootstrapped test gives slightly smaller power than the asymptotic test, but the difference is marginal especially in view of the large improvement in size. Table 2 reports the results for the GLS-Wald test. Overall, the GLS-Wald test using the asymptotic critical values are more liberal than the corresponding LR test, which is a well known characteristic of Wald tests. Despite being liberal, the GLS-Wald is less sensitive to the value of α, since we are using a robust covariance matrix. The GLS-Wald test also gets conservative as δ i gets large for the same reason as the LR test. The bootstrapped GLS-Wald test controls the size as well as the bootstrapped LR test. Table 3 shows the results for the OLS-Wald test. It performs similarly to the GLS-Wald test for both the asymptotic and bootstrap versions. Lastly, the power is close to one in almost all cases. The hypotheses that are rejected misspecifies the break dates by only 5% of the sample size, yet our tests are powerful enough to detect them. To sum up, all three tests, the LR, GLS-Wald, and OLS-Wald suffer from size distortions if the breaks are not large enough. However, the bootstrap procedures significantly help control the size without losing much power. We investigate the commonality of the break dates across the temperature and anthropogenic total radiative forcing series as well as test whether the recent "hiatus" is significant. We use two sets of global, northern and southern hemispheric temperature series, each of which will be denoted as G, N and S, respectively. The first set comes from the Climate Research Unit's HadCRUT4 (Morice et al., 2012) and the second from the NASA database (GISTEMP Team, 2015;Hansen et al., 2010). For global temperatures, we also use the data from Berkeley Earth (Rohde et al., 2013) and the dataset in Karl et al. (2015). We first discuss the data used, then the results for the common breaks tests and finally the results for the test on a change in slope in temperatures related to the so-called "hiatus" period. The annual temperature data used are from the HadCRUT4 (1850-2014) (http://www.metoffice. gov.uk /hadobs/hadcrut4/data/current/download.html) and the GISS-NASA (1880-2014) datasets (http://data.giss.nasa.gov/gistemp/). The Atlantic Multidecadal Oscillation (AMO) and the North Atlantic Oscillation (NAO) series (1856-2014) are from NOAA; http://www.esrl .noaa.gov/psd/data/timeseries/AMO/ and http://www.esrl.noaa.gov/psd/ gcos wgsp/ Timeseries/NAO/). As stated above, for global temperatures, we also use the data from Berkeley Earth (Rohde et al., 2013) and the dataset in Karl et al. (2015). We also use series from databases related to climate model simulations by the Goddard Institute for Space Studies (GISS-NASA). The radiative forcing data obtained from GISS-NASA (https://data.giss.nasa.gov/ modelforce/; Hansen et al., 2011) for the period 1880-2010 include the following (in W/m2):

well-mixed greenhouse gases, WMGHG, (carbon dioxide, methane, nitrous oxide and chlorofluorocarbons); ozone; stratospheric water vapor; solar irradiance; land use change; snow albedo; stratospheric aerosols; black carbon; reflective tropospheric aerosols; and the indirect effect of aerosols. The aggregated radiative forcing series are constructed as follows: WMGHG is the radiative forcing of the well-mixed greenhouse gases and has a largely anthropogenic origin; Total Radiative Forcing (TRF) is WMGHG plus the radiative forcing of ozone, stratospheric water vapor, land use change; snow albedo, black carbon, reflective tropospheric aerosols, the indirect effect of aerosols and solar irradiance. The temperature series are affected by various modes of natural variability such as the Atlantic Multidecadal Oscillation (AMO) and the North Atlantic Oscillation (NAO), which are characterized by low frequency movements (Kerr, 2000;Hurrell, 1995).
Since trends and breaks are low frequency features, it is important to purge them from the temperature series allowing more precise estimates of the break dates.2 Other high frequency fluctuations in temperature series do not affect the precision of the estimates of the break dates and the magnitudes of the changes in slope. Accordingly, we filter out the effect of these modes of variability by regressing each temperature series on these modes and a constant. Since the effect of natural variability might have occurred with a time lag, we choose an appropriate lag using the Bayesian Information Criterion (BIC); Schwarz (1978). The candidate regressors for the filtering are the current value and lags (up to order kmax -1) of AMO and NAO. We first work with G from HadCRUT4. We start with kmax = 2, so that the candidates are the current value and the first lag only. BIC chooses the current value of AMO and the first lag of NAO. Since the maximum lag allowed is selected, we increase kmax to 4. Then, BIC chooses the current value of AMO and the second lag of NAO (the first is not included given that we search over models not necessarily including all lags up to some chosen order). When applying the BIC, the number of observations used is limited by kmax (e.g., Perron and Ng, 2005). Having decided on the current value of AMO and the second lag of NAO, we apply the filtering to all available observations, not limited by kmax. We could repeat the same procedure to each series. However, it does not make much sense to have the same mode affect each temperature series with a different lag. Hence, we filter all series with the current value of AMO and the second lag of NAO. 3 The filtered temperature series are denoted as G, Ñ and S. Figure 1 presents graphs of the original and filtered series. For the radiative forcing variables, we use WMGHG and TRF. Figure 2 displays the time plots of WMGHG and TRF as well as the AMO and NAO. The data set used covers the period 1880-2014 and the filtering is done with the full sample. We use the data for 1900-2014 for the purpose of the common break date tests. 4 We use two sample periods: 1900-1992 and 1963-2014. The reason we split the sample is to reduce the number of breaks in the system, thereby lessening the computational burden. Our choice of sample periods is made with the intention to minimize any unwanted effect from a second break that might exist in any of the variables in the system. The break dates in WMGHG are estimated most precisely since it has a very small noise component (see Figure 2) and they are 1963 and 1992. First we consider bivariate systems with one temperature series and one radiative forcing series. The null hypothesis of a common break date is tested using the LR test from Theorem 1. Since this test is not asymptotically pivotal, we supply bootstrap p-values as well. The bootstrap is carried out in the same way as done in the Monte Carlo simulations reported earlier. We also report the GLS-Wald test from Theorem 2. Although the GLS-Wald test is asymptotically pivotal, we still supply the bootstrap p-values because of the size distortions observed in our Monte Carlo study. In addition to the 4 global, 2 northern hemispheric and 2 southern hemispheric temperature series, we also consider the average of the global, northern and southern hemispheric temperature series across different datasets. Hence, there are 11 temperature series and 11 filtered temperature series, each of which is paired with either WMGHG or TRF, yielding a total of 44 bivariate systems. Table 4 (a) reports the results for the sample period of 1900-1992 when the original temperature series are considered; the asymptotic p-values of the Wald and LR test are very small in most cases and the null of a common break date is rejected with only a few exceptions. However, the bootstrap p-values are well over the conventional thresholds 0.05 or 0.01 in all cases. Given the previous simulation results, we view the bootstrap p-values as more reliable than the asymptotic ones. Note that the break date for the radiative forcing variable is almost unchanged at either 1963 or 1966, regardless of the temperature series in the bivariate system. Also, the estimated common break date, i.e., the estimate imposing the null of common breaks, always coincide with the break date for the radiative forcing variable. On the other hand, the break date estimates for the temperature series are quite spread out. Hence, the break dates in the radiative forcing variables are being estimated more precisely than those in the temperature series. The large p-values produced by the bootstrapping is a consequence of the high level of noise in the temperature series relative to the magnitude of the break. Accordingly, we fail to reject the null of a common break date despite the large discrepancy in the estimated break dates. dates for the filtered temperature series are 1990-1992 with only a few exceptions. Also, the p-values are larger than the corresponding ones for the original temperature series. We repeat the analysis for the bivariate systems with the OLS-Wald test as a robustness check and report the results in Table 5. The results are very similar when the filtered series are considered. The break date estimates obtained equation by equation coincide with those obtained from the system or differ only by a few years at most. When the original series are used, the break date estimates for some series largely deviate from those computed by the system method. However, the fact that the bootstrap p-values are much larger than usual thresholds still applies. Thus, our conclusions remain the same. An advantage of the OLS-Wald test is that it can handle a larger system without extra computational burden.
The break date estimates will not change if we consider a larger system, because they are computed equation by equation. However, the long-run variance estimate and the bootstrap p-values can change as we consider a bigger system. We consider three five-variables systems. Each system has a global, northern and southern hemispheric temperature and the two radiative forcing variables. The first one uses the temperature series from HadCRUT4, the second one from NASA and the third one has average global, northern and southern hemispheric temperature series. The results are reported in Table 6. They show results broadly similar to those obtained bivariate systems. To provide a final robustness check, we consider testing jointly whether the two breaks near 1960 and for the "hiatus" are common. We use the OLS-Wald test, which estimates the two breaks jointly for each equation, and the period 1926-2014 to focus on the time span whereby anthropogenic factors are most likely to be the main driver affecting changes in temperatures. The results are presented in Table 7. The results obtained are broadly similar and lead to the same conclusions. When using the filtered series and the bootstrap p-values, the test fails to reject the null hypothesis of two common breaks. The first is located in the early 60's (the onset of sustained global warming) and the second in the early 90's (the "hiatus") consistent with methods dealing with one break at a time presented earlier. The literature on the existence and drivers of the current slowdown in the rate of warming has expanded quickly (Tollefson, 2014(Tollefson, , 2016)). In this section, we focus on testing if the "hiatus" can be explained by natural variability or if it is a feature of the underlying warming trend (Estrada, Perron and Martinez, 2013). A small part of the literature investigating the slowdown in the warming in global temperatures are based on formal structural change/change-point tests (Cahill et al., 2015;Pretis et al., 2015), while others are simply based on testing for trends during arbitrarily selected periods of time (Cahill et al., 2015;Foster and Rahmstorf, 2011;Lewandowsky et al., 2015). This type of approach, as well as most of physical-based studies, have typically failed to find convincing evidence for the existence of a change in the slope of global temperatures during the last decades. According to one of the commonly accepted hypothesis, the apparent "hiatus" is produced by the effects of low-frequency natural variability that result from coupled ocean-atmosphere processes and heat exchange between the ocean and the atmosphere. It has been proposed that effects of natural variability modes such as AMO, NAO and PDO were able to mask the warming trend since the 1990s, creating the illusion of a slowdown in the underlying warming trend (e.g., Guan et al ., 2015;Steinman et al., 2015;Li, Sun and Jin, 2013;Trenberth and Fasullo, 2013). However, Estrada and Perron (2016) argue that low-frequency oscillations do distort the underlying warming trend but cannot account for the current slowdown. Instead, the main effect of these oscillations has been to make it more difficult to detect the drop in the rate of warming, a real feature of the warming trend imparted by the slowdown in the radiative forcing from well-mixed greenhouse gases. Establishing whether the "hiatus" period is statistically significant is important because it would counter the widely held view that it is the product of natural internal variability (Kosaka and Xie, 2013;Trenberth and Fasullo, 2013;Meehl et al., 2011;Balmaseda, Trenberth and Källén, 2013). Using standard tests (e.g., Perron and Yabu, 2009), the results are mixed across various series and sometimes borderline. Our aim is to provide tests with enhanced power by casting the testing problem in a bivariate framework involving temperatures and radiative forcing. We consider bivariate systems with one temperature series and one forcing variable. We use the LR test for the presence of a break in temperature series given the presence of a structural break in radiative forcing series. We first estimate a break date in the radiative forcing series under the likelihood framework for a bivariate system of a pair of radiative forcing and temperature series. Then, given the estimated break date in the radiative forcing series, we apply the LR test for the null hypothesis of no break in temperatures versus one break. To obtain p-values, we apply the bootstrap procedure described in Section 5 with 1,000 samples. The results are presented in Table 8. First, as expected from previous results, the LR test is more likely to reject the null hypothesis of no break when using filtered temperature series. Hence, we report results only for that case. For the filtered series from 1900 to 1992, we reject the null at less than 5% significance level in all pairs of radiative forcing and temperature series. Hence, there is clear evidence of a break in temperatures that is near 1960 (varying between 1954 and 1966 depending on the series and the forcing). This concurs with the common break results obtained in the previous section. When using the sample 1963-2014 we reject the null in seven and eight cases out of eleven filtered series at the 10% significance level with well-mixed green-house gases and total radiating forcing, respectively. Since the sample size is small (T = 52), we might expect that the LR test may have little power, but our result suggests strong evidence for the presence of a break in temperature series. Note that the evidence for a break is stronger when using bivariate systems involving TRF. This is due to the fact that TRF exhibits a larger decrease in slope compared to WMGHG. The errors are also more strongly correlated. For instance, solar irradiance is a part of both TRF and temperatures and an important source of variations. When considering systems with TRF, the only pairs that do not allow a rejection at the 10% level are those associated with Northern hemisphere temperatures. For these, the p-values range from .14 to .19, while for all other pairs they are below 10%.
Given that the small sample size is small, the overall evidence strongly indicates a break in temperatures in the early 90s consistent with the presence of the much-debated "hiatus". We consider a multivariate system with n equations with the dependent variables modeled as joint-segmented trends with multiple changes in slope. The errors can be serially correlated and correlated across equations. We consider testing for general linear restrictions on the break dates, including testing for common breaks across equations. The test used is a (quasi-) likelihood ratio test assuming serially uncorrelated errors. Under the stated conditions, the LR test has a pivotal chi-square distribution, though it is non-pivotal in general cases of interest. Hence, we also consider a corrected Wald test having a pivotal limit distribution, which can be constructed using break dates estimated one equation at a time or via the complete system. The limit distribution of the test is standard chi-square. However, simulations show that the two Wald tests suffer from potentially severe liberal size distortions for moderate to strong serial correlation. Hence, for all three tests, we suggest a bootstrap procedure to obtain the relevant critical values. These bootstrap tests are shown to have correct size in all cases considered. Our empirical results show that, once we filter the temperature data for the effect of the Atlantic Mutidecadal Oscillation (AMO) and the North Atlantic Oscillation (NAO), the breaks in the slope of radiative forcing and temperatures are common, both for the large increase in the 60s and the recent "hiatus". We also consider a test for the presence of an additional break in some series. The theoretical framework is general and allows multiple breaks in a general multivariate system. The test considered is a quasi-likelihood ratio test. The limit distribution is shown to be non-pivotal and depends in a complex way on a number of nuisance parameters. Hence, we use a bootstrap procedure to obtain relevant critical values. In our applications, the null and alternative hypotheses specify a break in the slope of the radiative forcing series, while under the null hypothesis there is no break in the temperature series but there is one under the alternative hypothesis. Our results indicate that indeed the "hiatus" represents a significant slowdown in the rate of increase in temperatures, especially when considering global or southern hemisphere series, for which our test points to a rejection of the null of no change for all data sources considered. The statistical results are of independent interest and applicable more generally beyond the climate change applications considered. Proof of Theorems 1 and 2: Let λr = arg min λ s.t. Rλ=r | Σ(λ)| and λ = min λ | Σ(λ)|. Given the consistency of Σ( λr ), Σ( λ) and Σ(λ 0 ) under the null hypothesis, we can apply a Taylor series expansion to LR such that

Consider the first term in the above decomposition. We may write

where U λ is the matrix of estimated residuals from the infeasible GLS estimator, instead of the ML estimator, for λ. For notational simplicity, write Ω = Σ ⊗ I T and M(λ

. Then the test statistic admits the following decomposition:

Note that

where

where the last equality follows from Lemma 1 of Perron and Zhu (2005). In addition, the O(T ) term in the above equation is not o(T ) and (XX) > 0. For (UU), note that

These orders of magnitude show that (XX), which is strictly positive, is the dominant term in (A.2). When these terms are evaluated at k = kr , (XX) is still the dominant term and it contradicts to the fact that tr

The only way to avoid the contradiction is to have T 1/2 || kr -k 0 || = T 3/2 ( λr -λ 0 ) = O p (1). Since the rate of convergence is obtained as T 3/2 , the minimization can be carried out over the set

/2 for all i, j and some large constant M}. On Λ T , (UU) is asymptotically negligible and can be ignored. Thus,

The second term in (A.1) is basically the same as the first except that there is no restriction on the break dates, so that Q = Z. It follows that

and

Note that

where

For Z ′ z, let ǫ = Ω -1 u, which means that ǫ t = (ǫ 1t , ..., ǫ nt ) ′ = Σ -1 (u 1t , ..., u nt ) ′ and the partial sums of ǫ t obeys a functional central limit theorem from Assumption 3, i.e., T -1/2 [T r] t=1 ǫ t ⇒ B(r) = Σ -1 Ψ 1/2 W (r). Denote the i th element of B(r) by Bi (r). We can then express the limit of Z ′ z as the follows:

where the covariance of the vector of stochastic integrals in the above expression is given by 

). Proof of Corollary 1: Similarly to (A.2), we have the following decomposition

The results in Perron and Zhu (2005) directly applies and λ i is consistent at rate T 3/2 . Thus, we can again focus on T -3/2 neighborhood of λ 0 i , say Λ iT , and

Define

where B i (r) is the i th element of B(r) = Ψ 1/2 W (r) in Assumption 3. Furthermore, the weak convergence holds jointly in i under Assumption 3. Therefore,

where

Proof of Theorem 3: We use Σ(i) (k, h) and Σ(i) (λ, ν) interchangeably. As in (A.1),

where ν = arg min ν∈C (i)

T log Σ(i) ( λ, ν) . Define the projection matrix P (i) (λ, ν) as

where a i = a i ([T ν]) for simplicity and M(λ

we can write

From the proof of Theorem 1, recall that

To further analyze this expression, we write

The weak convergence result for the numerator of (A.3) follows from Assumption 3 and the continuous mapping theorem. To see this, consider the following decomposition.
For A 1 , note that Z ′ Z p → Ξ 1 from the proof of Theorem 1 and .5) uniformly in ν. For A 3 , let Bi (r) be as defined in the proof of Theorem 1. Then,

over ν ∈ (0, 1). Thus, A 3 is a Gaussian process with covariance matrix given by 

Therefore, combining (A.4), (A.5) and (A.6) yields the weak convergence result

over ν ∈ (0, 1). Lastly, the denominator in (A.3) is such that

uniformly in ν. The result in part (i) of the theorem follows from (A.7) and (A.8). The result in part (ii) also follows since the above result holds jointly in i.

A-6 Note: Each temperature series is denoted by two letters. For the first one, G, N and S stand for global, northern and southern hemispheric temperatures, respectively. For the second one, H, N, B and K stand for HadCRUT4, NASA, Berkeley and Karl. Avg stands for average. The filtered temperature series are denoted with a tilde. W and TRF stand for well-mixed green-house gases and total radiative forcing, respectively. Entries less than 0.01 are marked as "-". T-4 Note: Each temperature series is denoted by two letters. For the first one, G, N and S stand for global, northern and southern hemispheric temperatures, respectively. For the second one, H, N, B and K stand for HadCRUT4, NASA, Berkeley and Karl. Avg stands for average. The filtered temperature series are denoted with a tilde. W and TRF stand for well-mixed green-house gases and total radiative forcing, respectively. Entries less than 0.01 are marked as "-". T-5 Note: Each temperature series is denoted by two letters. For the first one, G, N and S stand for global, northern and southern hemispheric temperatures, respectively. For the second one, H, N, B and K stand for HadCRUT4, NASA, Berkeley and Karl. Avg stands for average. The filtered temperature series are denoted with a tilde. W and TRF stand for well-mixed green-house gases and total radiative forcing, respectively. Entries less than 0.01 are marked as "-". The numbers in parentheses are break dates. T-6 Note: Each temperature series is denoted by two letters. For the first one, G, N and S stand for global, northern and southern hemispheric temperatures, respectively. For the second one, H, N, B and K stand for HadCRUT4, NASA, Berkeley and Karl. Avg stands for average. The filtered temperature series are denoted with a tilde. W and TRF stand for well-mixed green-house gases and total radiative forcing, respectively. Entries less than 0.01 are marked as "-". The numbers in parentheses are break dates. hemispheric temperature from HadCRUT4 and the average of the two southern hemispheric temperature series, for which the estimate is still 1956. As a result, both the asymptotic and bootstrap p-values are much larger than those for the original temperature series.Tables 4 (b) presents the results for the second sample 1963-2014 when the original temperature series are used and the estimates of the break dates for the temperature series are again quite spread out, while those for the radiative forcing variables are consistently around 1990-1992. Just like in Table4(a), the asymptotic p-values suggest strong rejection of the null hypothesis, but the bootstrap p-values turns out much larger reflecting again the high level of correlated noise in the original series. Note that the estimates of the break

hemispheric temperature from HadCRUT4 and the average of the two southern hemispheric temperature series, for which the estimate is still 1956. As a result, both the asymptotic and bootstrap p-values are much larger than those for the original temperature series.Tables 4 (b) presents the results for the second sample 1963-2014 when the original temperature series are used and the estimates of the break dates for the temperature series are again quite spread out, while those for the radiative forcing variables are consistently around 1990-1992. Just like in Table4

hemispheric temperature from HadCRUT4 and the average of the two southern hemispheric temperature series, for which the estimate is still 1956. As a result, both the asymptotic and bootstrap p-values are much larger than those for the original temperature series.Tables 4 (b) presents the results for the second sample 1963-2014 when the original temperature series are used and the estimates of the break dates for the temperature series are again quite spread out, while those for the radiative forcing variables are consistently around 1990-1992. Just like in Table

hemispheric temperature from HadCRUT4 and the average of the two southern hemispheric temperature series, for which the estimate is still 1956. As a result, both the asymptotic and bootstrap p-values are much larger than those for the original temperature series. Instead of using λ under both H 0 and H 1 , it is possible to jointly minimize (λ, ν) under H 1 . However, this requires a more complex asymptotic analysis and we will focus on the simpler case. We also considered filtering the effect of the Pacific Decadal Oscilliation (PDO) index (https://www.ncdc.noaa.gov/teleconnections/pdo/) and the Southern Oscilliation Index (SOI) (http://www.cru.uea.ac.uk/cru/data/soi/); seeWolter and Timlin (1998) andMantua and Hare (2002). The results were qualitatively similar showing robustness. Hence, we do not report them. If we choose the filtering regressors for each series, G and SH from NASA and G from Berkeley require the fourth lag of NAO instead of the second. But the break date estimates change by two years at most. For the radiative forcing variables, the data is available only up to 2011. We use forecast values for 2012-2014 based on an autoregressive model with a broken linear trend. In the case of TRF the forecast is produced using an AR(2) and a broken linear trend with 1960 and 1991 as the break dates. For WMGHG, we use an AR(8) and and a broken linear trend with break dates at 1960 and 1994.