Global Average Temperature (GAT) observed since the late 1800s exhibits an upward trend widely interpreted as evidence of global warming (Mann et al., 1998;Gadea and Gonzalo, 2020;AR6-IPCC, 2021). The specific nature of the trend is an open question in the empirical literature with a non-trivial answer. Neither theoretical climate nor economic models with climatic variables help to identify a particular trend specification. From the climate quantitative perspective, two dominant strands of literature with opposing views debate. On one side, authors as Woodward and Gray (1995), Kaufmann et al. (2010), Dergiades et al. (2016), Bruns et al. (2020), and Phillips et al. (2020) approximate the GAT series by an integrated process of order 1, I(1), and advocate for the use of modelling approaches accounting for the presence of stochastic trends. On the opposing side, Seidel and Lanzante (2004), Gay-Garcia et al. (2009), Estrada et al. (2013), Estrada and Perron (2017), Mudelsee (2019), or McKitrick et al. (2023) defend that GAT can be best represented as trend-stationary process with a one-period structural break in the trend function. The assumption of deterministic trends is also common in the climate analysis as inferred from the statistical procedures implemented by the Intergovernmental Panel on Climate Change (IPCC) on its Sixth Assessment (AR6-IPCC, 2021) and previous reports. Determining whether GAT is an I(1) process or a stationary process around a non-linear trend is a key issue for detection, attribution, and impact studies of climate change (see McKitrick et al. 2023). If GAT is assumed to be I(1), attribution of global warming involves demonstrating that the stochastic trends in temperature are inherited from its association with the anthropogenic forcing from CO2 and other greenhouse gases. 1 However, the I(1) assumption also implies that exogenous temporary shocks like solar flares or volcanic eruptions generate long-lasting effects on temperature, which does not seem to be the case in the observed record. 2

If trends are assumed to be deterministic, detection of global warming is possible through traditional trend-tests (see Gadea and Gonzalo (2020) for a robust trend-detection analysis). But this assumption poses challenges for attribution and impact studies due to the problem of 'unbalanced' relations. Concretely, in attribution studies, the mismatch in the order of integration between GAT and anthropogenic forcing hinders the estimation of the climate sensitivity parameter (for a review, see Rohling et al. (2018)) using regression analysis. This is also the case for the impact studies that rely on dynamic growth equations linking the growth rate of per-capita output with temperatures (see Dell et al. (2012) and references therein). Beyond this literature, understanding the nature of the trends is helpful for producing more accurate long-term forecasts as stressed by Kaufmann et al. (2010). The existing evidence regarding the nature of trends in GAT relies on the analysis of aggregated series. This paper contributes to the debate by offering a 'micro-founded' explanation for the observed trends based on the study of the trend dynamics of individual units used in computing those averages. On aggregate, our findings are consistent with the hypothesis of Seidel and Lanzante (2004) and Gay-Garcia et al. (2009), suggesting that temperature averages are trend-stationary with a one-time permanent shock breaking the trend function. Unit root tests implemented on individual grids provide evidence in the same direction. Our empirical analysis points towards the standard aggregation process as being important in explaining the origin of trend-breaks in aggregated series when combined with an accelerated global warming. Depending on the strength of the break's signal, these situations may bias standard unit root tests towards the non-rejection zone. We illustrate our hypothesis through a set of Monte-Carlo simulations assuming linear and broken-trend individual processes and emulating the standard aggregation methods used to compute the GAT. 1 Pretis and Hendry (2013) and Bennedsen et al. (2023) study the statistical properties of radiative forcing from different components (including CO2) and find that the series follows an I(1) or even an I(2) process. 2 Notice that a shock with transitory effects on GAT can have permanent effects on other variables. the interior of South America and Africa, as well as the Antarctica. To align with standard aggregation method followed by CRU, in our empirical analysis we use directly the series of gridded temperature anomalies from the period 1961-1990. Two alternative methods to compute average temperatures are considered. Method A aggregates the full set of grids with non-missing data available each year. Denoting the resulting average as T A t , for each t:

where T it is the temperature in grid i at year t, I it is an indicator for T it non-missing, and

Due to the non-uniform observation process for grids, this approach uses data from more grids in recent years compared to the earlier part of the record. i.e., N t grows with t. Method A closely resembles the standard aggregation procedure adopted by CRU. Method B, on the other hand, selects the set of grids with non-missing data throughout the entire sample period, ensuring a stable number of grids on the computation of the average. Denoting the resulting average as T B t , for each t:

where S = {i : I it = 1, ∀t} and |S| is the cardinality of S. This is the approach adopted by Gadea and Gonzalo (2020) to estimate not only the mean but any distributional characteristic of interest. Separate averages are calculated for the Northern Hemisphere (NH) and Southern Hemisphere (SH), using data from 1880 to 2022.5 Global temperature is obtained as a weighted average of both hemispheres, with the weights accounting for the difference in land areas. The estimated series, presented in Figure 2, exhibit increasing trends that are indicative of global warming. However, two distinctions emerge depending on the estimation approach. First, the series in Panel (b) exhibits higher variability explained by the smaller number of series used each year to compute averages under method B.
Second, in Panel (a), the average in the SH follows a distinct trend compared to the NH and the globe, specially in the second part of the record; this is not observed in the series of Panel (b). The nature of the trends in aggregated temperature series is analyzed using time series methods. We initially implement standard Augmented Dickey Fuller (ADF) tests including an intercept and a linear trend and selecting the number of lags based on the Bayesian Information Criterion (BIC). The test-statistics reported in Table 1 indicate that the null of unit root cannot be rejected for the globe and NH averages estimated under method A. 6 In contrast, the null of unit root is rejected in all aggregated series computed following method B.7

A potential shortcoming of the ADF test is its sensitivity to the specification of the deterministic component under the alternative. After Perron (1989), it is well known that in the presence of a shift in the trend function, the first order autorregresive coefficient is highly biased towards unity and the unit root null is hardly rejected. To address this situation, we complement our analysis with the Kim and Perron (2009) (KP) test allowing for a break in the trend function at an unknown period. As observed in Table 1, for the globe and NH temperatures computed with method A the null of unit root is now rejected. The break in the trend is detected to occur around 1964. For the other series, the test also detects a break but the conclusion of no-unit roots remain the same as before.8

In summary, the aggregated evidence aligns with the hypothesis of Seidel and Lanzante (2004) and Gay-Garcia et al. (2009), who argue that temperature process is trend stationary with a shift in the trend function. Unit roots non-rejected by the standard ADF test for method A averages are the consequence of a misspecified model for the deterministic component under the alternative that ignores potential structural breaks. Even though the series estimated under Method B also present a trend-break, the signal is not strong enough to drive the ADF test to the non-rejection zone. All this evidence implies that the aggregation method may have an impact on the conclusion of the tests. In this section, we model the trend dynamics of the individual grids used to compute temperature averages and analyze the consequences of the aggregation method for the emergence of structural breaks. Individual tests strongly reject the presence of unit roots. Results reported in Table 2 indicate that the proportion of rejections in the ADF test is high across the three sample periods considered. For instance, in a panel of 354 grids continuously observed during the period 1920-2022, the null is rejected in 88.17% of the cases. This proportion increases further when implementing the KP test. 9

Assuming that individual trends are deterministic, rather than stochastic, we first model each individual grid through a linear-trend model of the form:

where e it = ρ i e it-1 + v it , |ρ i | < 1, and v it uncorrelated. For each sample period, we estimate 9 Panel unit root tests point towards the same conclusion. Equation 3 and generate density plots of β1 . Panel (a) in Figure 3 shows that the density shifts to the right in more recent samples, reflecting the well-known acceleration in the global warming process (Mann et al., 1998;AR6-IPCC, 2021;Hansen et al., 2023). Focusing solely on the estimated slopes for the sample period 1960-2022, the mean of the coefficient is higher for the set of grids that appear later in the record. Acceleration in global warming, combined with the standard aggregation method A, can trivially generate the breaks in aggregated series even when the individual units follow a linear model. Consider the following extreme case. During the initial part of the record, a set of N 1 units with an average trend-slope β1 are observed. Assume that at a certain period, T B, a different set of N 2 units start to be observed with average slope β2 , β2 > β1 . The slope of the average series is β1 before T B, and (N 1 β1 + N 2 β2 )/(N 1 + N 2 ) after. 10 Depending on the weights, the result of the ADF test can bias towards non-rejection. In the next section, we illustrate this point using simulations. Alternatively, we explore trend-breaks within individual grids in a model of the form:

where Aggregating individual grids with broken trends using either Method A or B results in averaged series with structural breaks, as those in Figure 2. However, as discussed previously, due to the small number of grids involved in the computation of averages under method B, the 10 For the average computed using Method B the slope is β1 the full period. signal of the break is not strong enough to wrongly drive the ADF test towards the conclusion of unit roots. For method A averages, the non-rejection of the ADF test can be attributed to both a higher signal, stemming from the utilization of more information each period, and the deepening of the break due to the acceleration in global warming that implies the inclusion of more grids in periods of higher trend. This hypothesis is supported by the simulations in the next section. The proposed 'micro-founded' explanations for the type of the trends in aggregate temperatures are validated heuristically using simulations. We consider two alternatives to simulate the nonmissing indicator I it :

• Alternative 1: A fixed number of units are non-missing during the whole sample. A second group of series are missing during the first T * periods, and non-missing from T * on. • Alternative 2: I it follows a Markov Switching (MS) process with two states (missing and non-missing) and transition-probability matrix P.

First, we analyze the case where individual units follow a linear-trend model.
Two groups of series are simulated:

where β2 1 > β1

1 and e j it = ρ j i e j it-1 + v j it , j = 1, 2. Observation alternatives 1 and 2 assume that stations entering later in the average computation are chosen randomly from groups 1 or 2. To account for warming and acceleration, we define observation alternatives 1* and 2* assuming that the series of group 2 (i.e. those with higher trend-slopes) appear as non-missing later in the record. The simulation parameters, including the transition probabilities, are set based on empirical evidence aiming to match the simulated and observed trends. Sample size is T = 150. The number of series are N 1 = 150 and N 2 = 850. The process is replicated R = 1, 000 times. We impose that 10% of units are observed from the initial period. Table 3 presents the proportion of non-rejections of the ADF and KP tests for each observation alternative and aggregation method. Consistent with our hypothesis, even though the series are generated without unit root, the ADF test do not reject the null hypothesis in 100% of the cases under observation alternatives 1* and 2* and aggregation method A. For aggregation Method B, the unit root is rejected in all cases. Moreover, the KP tests rejects the unit root in all cases. Notes: The table reports the proportion of times that the null of unit root is nonrejected. Tests at 5% of significance. Individual series simulated using the following calibration: β 1 0i ∼ N (-0.7234, 0.4266 2 ), β 1 1i ∼ N (0.0108, 0.0043 2 ), β 2 0i = -2.36, β 2 1i ∼ N (0.0271, 0.0145 2 ), ρ 1 i ∼ N (0.264, 0.1176 2 ), ρ 2 i ∼ N (0.155, 0.1747 2 ), and v j it ∼ N (0, 3 2 ), j = 1, 2. These values are based on the estimates of fitting a linear model to individual grids over different samples. Next, lets consider the case where individual units contain one break in the trend function. Similar to the previous case, we simulate two groups of series of the form:

T

where the variables and parameters are defined in Equation 4, e j it = ρ j i e j it-1 + v j it , j = 1, 2, and we impose γ2 2 > γ1 2 to the capture warming acceleration. As in the previous case, alternatives 1* and 2* assume that the series in group 2 start to count for the average computation later. The parameters are set based on empirical evidence, and the sample size, number of series, and replications are the same as before. Table 4 presents the proportion of non-rejections of the ADF and KP tests. Consistent with our hypothesis, when averages are computed using Method A, the unit root in ADF test is non-rejected in 100% of the cases under all observation alternatives. For Method B, even though the original series contain a break and it is inherited by the aggregated series, the signal of the break is not strong enough to drive the test towards the non-rejection zone. However, if the proportion of units that are observed from the beginning of the record is increased, this signal eventually becomes stronger and the ADF test concludes the presence of unit roots. ∼ N (-0.6061, 0.5642 2 ), α 1 1i , α 2 1i ∼ N (-0.5524, 0.6899 2 ), γ 1 1i , γ 2 1i ∼ N (0.0110, 0.0174 2 ), γ 1 2i , ∼ N (0.0182, 0.0233 2 ), γ 2 2i , ∼ N (0.0271, 0.0145 2 ), ρ j i ∼ N (0.1101, 0.0899 2 ), and v j it ∼ N (0, 3 2 ), j = 1, 2. These values are based on the estimates of fitting a model allowing for a break in the trend to individual grids. Aggregate and micro-founded evidence do not support the hypothesis of stochastic trends in temperature. In fact, our evidence suggests that temperatures averages are stationary around a non-linear trend, with the non-linearity being modelled as a one-time break in the trend function. The break is attributed to a combination of individual grid breaks and the standard aggregation method under warming acceleration. Moreover, the aggregation method is relevant to bias ADF-tests towards the non-rejection zone by amplifying the signal of the break. Our findings carry important empirical implications for studies on the detection, attribution, impact, and forecasting of global warming. Radiative forcing from anthropogenic greenhouse gases (such as CO2, methane, and nitrous oxide) is modelled as an I(1), or even I(2) process. It introduces a problem of 'unbalanced' regressions if the aim is to establish human influence on global warming using cointegration methods. A similar issue is present on impact studies relating economic growth and temperature data. A valid alternative includes exploring co-trending approaches (see Chen et al. (2022)) or adopting statistical methods robust to the types of trends. Notes:

Notes:

The dataset is accessible at the following URL: https://crudata.uea.ac.uk/cru/data/temperature/. Similar results are obtained if we use raw-stations data. Only 155 grids are observed over the full sample period, with the majority of them located in the NH. This is the main evidence forKaufmann et al. (2010) and related literature to approximate the trend process through stochastic trends. Standard ADF-test implemented on other distributional characteristics (quantiles, skewness, kurtosis, etc.) do not show evidence of unit roots neither. Alternatively, we considered to model the deterministic component using polynomial trends of order 2 and 3. In both cases, the variation of the series around the trend function rejects the unit root hypothesis.