According to the most recent International Panel on Climate Change report, humanity is unlikely to prevent global warming by 1.5 • above pre-industrial levels. Still, aggressive curbing of greenhouse-gas emissions and carbon extraction from the atmosphere could limit its rise and even bring it back down (IPCC (2022)). But this window is rapidly closing, and, above the 1.5 • threshold, the chances of tipping points, extreme weather, and ecosystem collapse will become even more sizeable. An environmental tipping point is when small climatic changes might trigger large, abrupt, and irreversible environmental changes and lead to cascading effects. Recent IPCC assessments suggest that tipping points might arise between 1 • and 2 • warming, and likely to manifest at current emissions levels if they have not already occurred. Well-known tipping points concern the Greenland and the West Antarctic ice sheets, the Atlantic Meridional Overturning Circulation (AM OC), thawing permafrost, EN SO, and the Amazon rainforest. Recent evidence suggests that melting ice sheets is accelerating because of warming air and ocean temperatures and less snowfall. Some studies indicate that the irreversible disintegration of the Greenland ice sheet could occur at 0.8 • and 3.2 • warming (Wunderling et al. (2021)). An unstoppable ice sheet melting in Antarctica would manifest at 2 • warming (DeConto et al. ( 2021)). Ice sheets melting adds fresh water to the North Atlantic, weakening the AM OC, one of the main global ocean currents, which is already in its weakest state in 1,000 years (Caesar et al. (2021)). Its shutdown would cause significant cooling along the US east coast and Western Europe, alter rainfall and cause more drying. At the current global warming pace, a 50% weakening of AM OC is expected by 2100, and a tipping point between 3 • and 5.5 • warming. Moreover, the Arctic is warming twice as faster as the planet on average, and it has already warmed 2 • , causing permafrost thawing, which releases CO2 and methane into the atmosphere. Available estimates point to 1400 billion tons of carbon frozen in the Arctic's permafrost, twice as much carbon already in the atmosphere, and a 2 • warming could even cause the thawing of 40% of the world's permafrost. The El Niño-Southern Oscillation or EN SO cycle is an oscillating warming and cooling pattern affecting rainfall intensity and temperatures in tropical regions. It can strongly influence weather in many parts of the globe. El Niño and La Niña are the warm and cool phases of the EN SO cycle, respectively. Oceans warming can trigger a tipping point in the EN SO cycle, increasing its variability and intensity and shifting its teleconnection eastward (Cai et al. (2021)). Extreme rainfalls and droughts will no longer occur in tropical regions but throughout the earth due to the destabilization of these natural oscillations. The Amazon rainforest has already lost about 17% of its tree cover. At the current rate of deforestation, the loss could reach 27% by 2030. Lovejoy and Nobre (2018) estimate the dieback of the Amazon Forest at 20%-25%; beyond this deforestation threshold, the rainforest would transform into a savannah, potentially releasing up to 90 gigatons of CO2. Some climate models already indicate that the Amazon will be a net generator of C02 by 2035, setting the dieback threshold at 3 • warming. Further uncertainty on the compound effect of the above phenomena arises from their potential interaction, allowing tipping points to occur even below 2 • warming. Overall, greenhouse gases generated by human activity over the last two centuries have driven the global trend temperature up. This temperature warming has widely impacted the natural environment and has raised the risk of irreversible changes of state with catastrophic consequences (see also Schellnhuber (2008)) and Solomon et al. (2009)). In this paper, we link the concept of an environmental tipping point to the statistical property of time irreversibility. A stationary process {Y t } T t=1 is said to be time-reversible if its statistical properties are independent of the direction of time. In other words, the vectors (Y 1 , Y 2 , . . . , Y T ) have identical joint distributions as (Y -T , Y -(T -1) , . . . , Y -1 ) for every integer T . Hence, a time-reversible process (TR) exhibits a temporal symmetry in its probabilistic structure. In the alternative circumstance, we have time irreversibility when the stochastic process behaves differently according to the direction of time considered. TR has been under investigation in various fields over the years, for instance, in the different branches of physics, where researchers have been investigating whether time has some preferred direction in explaining physical phenomena (see Wald (1980), Levesque and Verlet (1993), Holster (2003)). This univocity along the time direction appears to be a tipping point property, as once a tipping point is reached, the system undergoes an irreversible state change. This paper aims to investigate whether TR has the potential to offer insight into the process of climate change and its implications for the natural environment. Studying TR in the context of climate change is motivated by the possibility of answering the following questions: are there divergences between the forward-time and backward-time joint probability distributions for the process of climate change and global warming? Are these processes symmetric over time? Is this property similarly present in natural oscillations that temperature warming might have permanently impacted, inducing changes in their frequencies and intensity of occurrence? Irreversibility in this context might carry insights into the event of state changes. This paper then introduces new strategies to detect whether a stochastic process is time-reversible. There are already several tests for TR in the econometric literature. See, for instance, Ramsey and Rothman (1996), Hinich and Rothman (1998), Chen et al. (2000), Belaire-Franch andContreras (2003), andProietti (2020). The shortcoming of many of these approaches is that they usually impose strong restrictions on the model or are not trivial to apply.
Our new strategies are grounded on mixed causal and noncausal models (see Gourieroux and Jasiak (2016)). Unlike causal models, which only consider the relationship between present and lagged values, mixed casual and noncausal models also compute the relationship between present and future values. This framework leads to nonlinear conditional expectations (e.g., Gourieroux and Jasiak (2022)). The connection between these models and TR gives rise to our testing strategies. Furthermore, similarly to Proietti (2020), we can test for TR on non-stationary time series using a novel approach. We extract the trend component using the Hodrick-Prescott (HP) filter imparted in a time-reversible closed-form solution. Then, the cyclical component, which records the process's oscillations around its trend, is responsible for the potential time-irreversibility feature of the stochastic process. The rest of the paper is as follows. Section 2 summarizes the properties of time-reversible processes and reviews the existing methods to detect TR. Section 3 introduces our new TR strategies. Namely, we show how our new approaches exploit the properties of mixed causal and noncausal models. We then evaluate their performance through Monte Carlo experiments. Section 4 extends our framework to non-stationary time series, and Section 5 presents the empirical assessment of some relevant climate variables. Finally, Section 6 concludes. (1975) shows that if a Gaussian error term characterizes an ARMA model, then the process is time-reversible. Indeed, Gaussian processes are entirely defined by their second-order moments, which have the property of being time symmetrical. Hallin et al. (1988) consider two-sided linear models of the form:

where the stationary condition

t=1 is time-reversible, then either t is a Gaussian white noise, or there exists a k and s ∈ {0, 1} such that θ 2k+j = (-1) s θ 2k-j . However, t has to be a sequence of i.i.d. zero-mean random variables with finite moments of all orders. It is an unrealistic assumption for non-Gaussian processes and many time series. Breidt and Davis (1992) extend Weiss's results to non-Gaussian processes assuming milder conditions than Hallin et al. (1988). They take the following ARM A(p, q) process into account:

where L indicates the backshift operator, φ(z) has r roots outside and s roots inside the unit circle (r + s = p), and t has a finite variance. For simplicity, we set the polynomial θ(L) = 1, such that (2) can be rewritten as:

where φ + (L) has r roots outside the unit circle while φ -(L) has s roots inside. It is well known that (3) has a unique stationary solution given by a two-sided moving average representation, as expressed in (1). Breidt and Davis (1992) claim that if φ(z) and φ(z -1 ) have different roots, then Y t is reversible if and only if the error term is Gaussian. In the other case, that is when the two polynomials φ(z) and φ(z -1 ) have the same roots, (1) (or equivalently (3)) is time-reversible regardless of the distribution of t . Indeed, if p > 0 and φ(z) and φ(z -1 ) have the same roots, 1/φ(z) has the Laurent expansion

with θ -p/2-j = θ -p/2+j , for j = 0, 1, . . . (see Breidt and Davis (1992)). This implies that the result of Hallin et al. (1988) is a consequence of the conclusion that the two polynomials φ(z) and φ(z -1 ) have the same roots. Moreover, unlike Hallin et al. (1988), Breidt and Davis (1992), only assume that the error term must have finite variance. Ramsey and Rothman (1996) define the stationary stochastic process {Y t } T t=1 is time-reversible only if:

for all i, j, k ∈ N + . This is a sufficient condition for TR, but not a necessary one since it only considers a proper subset of moments from the joint distributions of {Y t }. Since it is impractical to show that (5) holds for any i, j, and k, they adopt a restricted definition of TR by imposing i + k ≤ m and k ≤ K. In particular, they restrict m = 3 so that the symmetric-bicovariance function is given by:

for all integer values of k. Ramsey and Rothman (1996) claim that i + j = 3 is sufficient to provide a valid indication of time irreversibility. Ramsey and Rothman (1996) also introduced a new procedure to test TR that became a standard approach to investigating business cycle properties such as asymmetry. It amounts to a TR test statistic distributed as a standard normal distribution:

with:

and:

for various integer values of k. Under the null hypothesis, we have a timereversible process. The pre-requisite of the test is that the data must possess finite first sixth moment. If the distribution lacks this property, the test size can be seriously distorted (see Belaire-Franch and Contreras (2003)). Chen et al. (2000) propose a new class of TR tests, which, unlike Ramsey and Rothman (1996), does not require any moment restrictions. This class of tests relies on the fact that if {Y t } T t=1 is a time-reversible process, then for every

The drawback of this approach is that it allows for testing the symmetry of X t,k for each value of k, but not jointly for a collection of k values, which would require a portmanteau test.2 Moreover, its implementation is not trivial. Similar reasoning is followed by Proietti (2020) since also his test is based on the idea that X t,k has to be symmetric for every k > 0. However, Proietti (2020) uses a weaker definition of TR as {Y t } T t=1 can also be non-stationary. This Section introduces new strategies to assess TR in stationary stochastic processes, exploiting the properties of mixed causal and noncausal models. Breidt et al. (1991) introduce mixed causal and noncausal models as expressed in equation (3). They define the polynomial φ -(z) as noncausal and the polynomial φ + (z) as causal.
A required condition for identifying the causal from the noncausal component is the non-Gaussianity of the innovation term. Lanne and Saikkonen (2011), rewriting the noncausal polynomial in (3) as a lead polynomial, start with a mixed causal and noncausal model expressed as:

where L -1 produces lead such that L -1 Y t = Y t+1 . A mixed causal and noncausal model represented in this way is denoted as MAR(r,s), where ϕ(L -1 ) is the noncausal polynomial of order s and φ(L) is the causal polynomial of order r. Exactly as representation (3), r+s = p is true even in this case. Purely causal and purely noncausal models are obtained setting respectively ϕ(L -1 ) = 1 and φ(L) = 1 (see Gouriéroux et al. (2013), Hencic and Gouriéroux (2015), Hecq et al. (2016), Fries and Zakoian (2019), Hecq and Voisin (2021), Giancaterini and Hecq (2022), and Fries ( 2021)). In ( 8), both causal and noncausal polynomials have their roots outside the unit circle, such that:

The tests for TR that we propose have the common feature of extending the results obtained by Breidt and Davis (1992) to the MAR(r,s) representation (8). This is possible if and only if Condition 3.1 is true. Condition 3.1 A stochastic process that can be expressed as a MAR model is time-reversible if and only if φ(z)ϕ(z -1 ) have the same roots as φ(z -1 )ϕ(z). Namely, when: r = s and φ i = ϕ i , f or i = 1, . . . , s.

This implies that MARs are time-reversible if and only if the causal polynomial has the same order and the same coefficients as the noncausal polynomial and vice versa. Remember that it is impossible to identify a MAR model under the Gaussianity of the innovation term. Hence, in that case, we have a time-reversible process (see Weiss (1975)). The first strategy aims to evaluate whether a stochastic process meets Condition 3.1. In particular, it uses a procedure similar to the one used to identify MAR models (see Lanne and Saikkonen (2011) and Hecq et al. (2016)). The procedure is as follows:

1. We estimate a conventional autoregressive process (also called pseudocausal model) by OLS, and the lag order p is selected using information criteria (for instance, AIC or BIC). 2. We test the normality in the residuals of the AR(p). If the null hypothesis of Gaussianity is not rejected, we cannot identify a MAR(r,s) model, and for the reasons above, we have a time-reversible process. Moreover, if the null hypothesis of normality is rejected and the estimated p is an odd number, the condition r = s can never be satisfied. According to Condition 3.1, this result would allow us to identify our process as timeirreversible. However, the selection of p might not be univocal and depend on the information criterion employed. As such, to have more robust results before proceeding to the next step, we increase p by one unit so that r = s is still possible. In the alternative case that p is an even number, we directly proceed to the next step. 3. We select a model among all MAR(r,s) specifications with r + s = p if p is an even number; otherwise r + s = p + 1. This step is performed using a maximum likelihood approach (see Giancaterini and Hecq (2022) and references therein). In the selection procedure, we also include the model given by the restricted likelihood that imposes commonalities in causal and noncausal parameters (the model with the same restrictions as in Condition 3.1). Note that when we compute the information criteria of the model with restricted likelihood, instead of estimating p parameters (or p + 1 if p is an odd number), we estimate p/2 of them (or (p + 1)/2), implying a smaller penalty term. Finally, we choose the model with the smallest information criteria. Consider a short example to illustrate how the strategy works. We suppose that we estimate a conventional AR model by OLS, and we reject the Gaussian hypothesis of the residuals, for instance, using the Jarque-Bera test. Furthermore, we assume we select the number of lags p equal to 2. To analyze whether our process is time-reversible, we then compute the log-likelihoods and then the information criteria of the following four models: MAR(2,0), MAR(0,2), MAR(1,1) as well as the MAR(1,1) with the restriction φ = ϕ. If the model with the smallest information criteria is the one with the restriction, we have a time-reversible process. We have a time-irreversible process in the alternative case where another model is selected. This approach allows knowing with a limited number of steps whether the process is time-reversible. Its shortcoming is that information criteria are very sensitive to the sample size, and model selection might not be robust to sample update or trimming. Moreover, model selection can depend on the information criterion employed, i.e., AIC rather than BIC, HQ, or others. Finally, even for the same information criterion, the value used for model selection can only slightly differ from values shown by either lower or higher-order alternative models. The second strategy we introduce is more robust concerning the sample and slight differences in the value of information criteria when models are compared. However, more steps are required to identify the TR of the process than for the previous approach. It requires the following steps: steps 1 and 2 are identical to what we described in 3.1; 3. We select a model among all MAR(r,s) specifications with r + s = p if p is an even number (otherwise r + s = p + 1). Then, we choose the one with the largest likelihood (since we are considering models with the same number of parameters). 4. If the selected model is the one with r = s (in our previous example, it was the MAR(1,1)), we compute a likelihood ratio test, taking into account the same restrictions as in Condition 3.1. If we do not reject the test's null hypothesis, we have TR.
On the other hand, if we reject the null hypothesis, we identify the process under investigation as time-irreversible. We now analyze the performance of these two strategies using Monte Carlo experiments. We take into account data-generating processes (dgp) defined by an error term with a skewed Student's-t distribution, generated by joining two scaled halves of the Student's-t distribution (see Fernández and Steel (1998)):

where I( ) and I(-) stand for the indicator function:

for the density function of a symmetric Student's-t, and γ ∈ R + . In case γ = 1, we have f ( ) = g( ), hence (10) is a symmetric Student'st with ν degrees of freedom. The assumption that the error term follows a Student's-t is not a particularly strong hypothesis. It is a distribution that offers a good summary of the features of other (non-Gaussian) fat-tailed and symmetric distributions. Furthermore, our Monte Carlo experiments consider N = 1000 replications, four different sample sizes, T = (100, 200, 500, 1000), and the following combinations of causal and noncausal coefficients:

• MAR(1,1) :φ 0 = 0.8, ϕ 0 = 0.8; time-reversible process;

• MAR(1,1) :φ 0 = 0.8, ϕ 0 = 0.5; time-irreversible process;

• MAR(1,1) :φ 0 = 0.8, ϕ 0 = 0.1; time-irreversible process;

• MAR(1,0) :φ 0 = 0.8; time-irreversible process. In our Monte Carlo study, we also include results obtained by Ramsey and Rotham's test, setting k = 2. Table 1 shows the frequencies with which the two new strategies and the test proposed by Ramsey and Rotham detect the processes as time-irreversible when γ = 1, p is known, and r and s are unknown. In particular, columns Strategy 1 and Strategy 2 indicate the percentage of times the stochastic processes are identified as time-irreversible when the strategies from Sections 3.1 and 3.2 are implemented. The last column, RR (1996), indicates how often we reject the null hypothesis of TR when the methodology proposed by Ramsey and Rothman (1996) is used. The Bayesian Information Criteria (BIC) is used in Strategy 1. The results exhibit that Strategy 1 detects TR with greater precision, but is "undersized" for large T . This is because the penalty terms can differ from a tiny number in a large sample. On the other hand, Strategy 2 looks consistent and performs better when the processes under investigation are time-irreversible (frequencies are not size-adjusted, though, which makes the results of Strategies 1 and 2 not easy to compare). Finally, the test proposed by Ramsey and Rothman (1996) clearly shows size distortion problems. This is not an unexpected result since, as previously stated, the test can show a seriously distorted size if the distribution lacks a finite sixth moment. The Student's-t distribution has a finite sixth moment for ν > 6. As a consequence, the power of the test also performs poorly for RR (1996). Our simulation studies also consider cases where the error term is characterized by γ = 1. In these scenarios, we simulate a process with a skewed error term and proceed as if γ = 1: Strategies 1 and 2 are followed assuming a symmetric Student's-t distributed error term. The results obtained under these new circumstances are similar to those in Table 1. This suggests that the test size and power are not sensitive to the eventual asymmetry of the error term. The outcomes are available upon request. Table 2 shows different results when p is assumed unknown. In this case, before implementing our strategies, we estimate a pseudo-causal model in each replica of our simulation study to capture the dynamics p. Since there is more uncertainty under these new conditions, the results are less precise with small sample sizes (T = (100, 200)). However, the table displays that the outcomes align with Table 1 for large values of T . The percentages displayed in the column RR(1996) of Table 2   Finally, to analyze the result sensitivity to the persistence level, we imple-ment new Monte Carlo experiments considering as dgp a MAR(1,1) with the following combinations of causal and noncausal coefficients:

• MAR(1,1) :φ 0 = 0.95, ϕ 0 = 0.95; time-reversible process;

• MAR(1,1) :φ 0 = 0.95, ϕ 0 = 0.5; time-irreversible process;

• MAR(1,1) :φ 0 = 0.95, ϕ 0 = 0.1; time-irreversible process;

• MAR(1,0) :φ 0 = 0.95; time-irreversible process. Even in this case, the outcomes are similar to those displayed in Tables 1 and2 and available upon request. The goal of this section is to detect TR in non-stationary processes {Y t } T t=1 that can be expressed as:

where f Y is a generic trend function, and cc Y is a stationary process that captures the cyclical fluctuations of Y around f Y . We show that whenever the trend component is computed using the HP filter, then f Y can be expressed as a time-reversible process. As a consequence, the potential time irreversibility of process Y would be captured by its cyclical component cc Y . In other words, whenever f Y is estimated by using the HP filter, model ( 11) is time-irreversible (or reversible) if and only if its cyclical component is irreversible (or reversible). The HP filter estimates the trend component through the following minimization problem (see Hecq and Voisin (2021)):

According to De Jong and Sakarya (2016), the optimization problem (11) has the following closed-form solution:

for t = 3, . . . , T -2. The λ parameter penalizes the filtered trend's variability; therefore, the higher its value, the smoother the trend component:

with either i = 2 (see Backus and Kehoe (1992)) or i = 4 (Ravn and Uhlig (2002)). It can be shown that ( 12) can be rewritten as:

where ψ 1 (λ) = 4λ λ+1 , and ψ 2 (λ) = -λ. For instance, for annual data, we can adopt λ = 6.25, implying

The results underline that the filter of the trend component is given by a timereversible MAR(2,2) polynomial minus a constant value.
Since the latter does not affect the symmetry over time of our process, and our goal is to investigate the time reversibility of f Y , we do not consider the constant term in our investigation. As a consequence, we can approximate f Y as follows:

Using the Laurent expansion as in (4), we have:

where because of the identity of the lead and lag polynomials, δ is symmetric over time. Hence, even if f Y is a non-stationary process, we can apply a weaker definition of TR and define it as time-reversible. This result implies that the potential time irreversibility (or reversibility) lies with the cyclical component of Y . To illustrate how our new strategies perform under the new conditions, we implement new Monte Carlo experiments where f Y is a random walk with drift:

with η ∼ N (0, 1), and cc Y as a MAR(1,1):

Finally, the process {Y } T t=1 is obtained by the sum of the two processes, that is:

Alternatively, we could have considered the following process as dgp:

However, the reason not to consider such a process is that, as expressed in (20), the resulting dgp implies that the first difference process (∆Y t ) is a MAR(1,1). This is not a realistic assumption because MARs are typically used to capture explosive bubbles, and the first difference operation eliminates most locally explosive behaviors (see Hecq and Voisin (2021)). In each replica of our Monte Carlo experiment, we simulate the non-stationary process {Y t } T t=1 , remove the trend component using the HP filter, and then apply our strategies on cc Y . The coefficients used for the cyclical component cc Y are the same as in the previous section. Table 3   The results are similar to those displayed in Tables 1 and2, with the difference that the power of the strategies is less accurate under these new conditions, especially when the sample size considered is small (T = (100, 200)). In our empirical investigation, we analyze annual data for the global land and ocean temperature anomaly (GLO), the global land temperature anomaly (GL), the global ocean temperature anomaly (GO), solar activity (SA), emissions of greenhouses gas (GHG), emissions of nitrous oxide (N 2O). When available, we also use monthly data to control for potential small sample distortions in our statistics, as revealed by the simulation results. In particular, we consider the following monthly series: the Southern Oscillation Index (SOI), the North Atlantic Oscillation Index (N AO), the Pacific Decadal Oscillation Index (P DO), the global mean sea level (GM SL), the Northern Hemisphere sea ice area (N H), the Southern Hemisphere sea ice area (SH), the global component of climate at a glance (GCAG), and, finally, the global surface temperature change (GIST EM P ). 3  GLO, GL, GO, GCAG, and GIST EM P measure global warming. They provide the difference between the current temperature from a standard benchmark value. Positive anomalies show that the observed temperature is warmer than the benchmark value, and negative temperatures show that the observed temperature is colder than the benchmark value. In particular, GCAG provides global-scale temperature information using data from NOAA's Merged Land Ocean Global Surface Temperature Analysis (NOAAGlobalTemp), which uses comprehensive data collections of increased global coverage over land (Global Historical Climatology Network-Monthly) and ocean (Extended Reconstructed Sea Surface Temperature) surfaces. SOI is one of the most important atmospheric indices for determining the strength of El Niño and La Niña events and their possible effects on weather conditions in the tropics and various other geographical areas. El Niño events are characterized by sustained warming of the central and eastern tropical Pacific, whereas La Niña events show sustained cooling of the same areas. These changes in the Pacific Ocean and its overlying atmosphere occur in a cycle known as the El Niñ o-Southern Oscillation (EN SO). High values of SOI indicate La Niña events, whereas negative values indicate El Niño events. The N AO determines the westerly winds' speed and direction across the North Atlantic and the winter sea surface temperature. When the N AO index is far above average, there is a greater likelihood that seasonal temperatures in northern Europe, northern Asia, and South-East North America will be warmer than usual. In contrast, seasonal temperatures in North Africa, North-East Canada, and southern Greenland will be cooler than usual. The opposite is true when N AO is far below average. P DO is a climatic cycle that describes anomalies in sea surface temperature in the Northeast Pacific Ocean. The P DO has the power to influence weather patterns all over North America. Finally, GM SL, N H, and SH are climate indicators providing information on how much of the ice land is melting, and their connection with global warming is straightforward. Figure 1 presents the data. GLO, GL, GO, SA, GHG, and N 2O range from 1881 to 2014, SOI and N AO from January 1951 to December 2021, P DO from January 1854 to December 2021, GCAG and GIST EM P from January 1880 to December 2016, GM SL from January 1880 to December 2015, and, finally, N H and SH from January 1979 to 2021. 3 GLO, GL, and GO are obtained from https://www.ncdc.noaa.gov/cag/global/time-series. SOI and N AO are obtained from https://www.cpc.ncep.noaa.gov/data/indices/soi and https://www.cpc.ncep.noaa.gov/products/precip/CWlink/pna/norm.nao.monthly.b5001.current.ascii.table, respectively. P DO is obtained from https://www.ncdc.noaa.gov/teleconnections/pdo/. For GHG and SA, the source is Hansen et al. (2017). For N 2O, we use the historical reconstruction computed in Meinshausen et al. (2017) (data available at https://www.climatecollege.unimelb.edu.au/cmip). N H and SH are obtained from https://psl.noaa.gov/data/timeseries/monthly, GCAG and GIST EM P from https://datahub.io/core/global-temp and, finally, GM SL from https://datahub.io/core/sea-level-rise. As shown in Figure 1, time series (a)-(i) are characterized by a positive trend. Hence, according to the strategy introduced in Section 4, their potential time-reversibility (or irreversibility) lies with their cyclical component. For this reason, we can remove their trend and extract their cyclical fluctuations using the HP filter.
Figure 2    The goal is to investigate the time-reversibility of the variables displayed in Figure 2 and the latter five in Figure 1. We estimate autoregressive models (see Sections 3.1 and 3.2) for each time series. We use the BIC to identify the number of lags (p). Next, we test the normality of the residuals of the nine AR(p) models. Since for cc GLO , cc GL , cc GO , cc SA , and SH we do not reject the null hypothesis of normality (significance level 0.05) of the Shapiro-Wilk test (p-values equal to 0.83, 0.59, 0,24, 0.08, and 0.45, respectively) and the Jarque-Bera test (p-values equal to 0.64, 0.25, 0.15, 0.13, and 0.35 respectively), we identify them as time-reversible processes. On the other hand, in cc GHG , cc N 2O , cc GCAG , cc GIST EM P , cc GM SL , SOI, N AO, P DO, and N H, we reject the null hypothesis of Gaussianity of both the Shapiro-Wilk test (p-values are close to zero for cc GHG , cc N 2O , cc GCAG , cc GIST EM P , cc GM SL , P DO, N H and 0.0362, 0.0020 for SOI and N AO, respectively) and the Jarque-Bera test (p-value equal to 0.0307 for SOI and close to zero for all the other variables) at a significance level of 0.05. We can then fit MAR models to our data, identifying cc GHG , cc GCAG ,and cc GIST EM P as MAR(2,0), cc N 2O as MAR(4,0), cc GM SL as MAR(6,2), P DO as MAR(0,4), SOI as MAR(2,2), N AO as MAR(1,1), and N H as MAR(12,2) (Table 4). Since the condition r = s is not met, GHG, N 2O, GCAG, GIST EM P , GM SL, and P DO are timeirreversible. However, TR is still a possible outcome for SOI and N AO; hence, we implement the next steps of Strategies 1 and 2. Since the information criteria of the restricted MAR(2,2) (BIC=2705.795) is larger than the one provided by the MAR(2,2) with no restrictions (BIC=2659.974), Strategy 1 identifies SOI as time-irreversible. The same result follows from Strategy 2: the null hypothesis of TR is rejected since the estimated likelihood ratio test statistic equals 52.57. Contrastingly, N AO is identified as time-reversible from both strategies: the information criteria of the restricted MAR(1,1) is lower (BIC=2450.267) than the one provided by the unrestricted MAR(1,1) (BIC=2456.317), and the estimated likelihood ratio test statistic is equal to 0.6985. Even if this last time series rejects the null hypothesis of Gaussianity, it is very close to the Gaussian case since the estimated degrees of freedom (ν) equals 96.2. However, identifying N AO as non-Gaussian does not affect our conclusions as both strategies identify it as time-reversible. In summary, our findings identify GLO, GL, GO, SA, N AO, and SH as time-reversible and GHG, N 2O, GCAG, GIST EM P , GM SL, SOI, P DO, and N H as time-irreversible. The time irreversibility of cc GHG and cc N 2O is a noticeable property of variables that account for the warming trend in global temperatures (IPCC (2014), Morana and Sbrana (2019)). We expect time irreversibility also to be present in other variables affected by greenhouse gas emissions, as, statistically, a linear combination of time-irreversible and timereversible variables is also time-irreversible. This result can explain why GCAG, GIST EM P , GM SL, SOI, P DO, and N H are time-irreversible. In particular, these results underline how global warming might have exerted feedback effects on natural oscillations, temperatures, and the environment in general. Among others, Morana and Sbrana (2019) show that GHG emissions are the key determinant of the warming trend in global temperatures. The irreversibility of GCAG and GIST EM P further corroborates these findings. Yet the evidence is inconclusive as the cyclical components of GO, GL, and GLO are time-reversible. However, whether this latter result might be an artefact due to their shorter sample is plausible. Morana and Sbrana (2019) also document Atlantic hurricanes' increasing natural disaster risk and destabilizing impact on the EN SO cycle. Indeed, oceans warming can trigger a tipping point in the EN SO cycle, increasing its variability and intensity and shifting its teleconnection eastward (Cai et al. (2021); see also Cai et al. (2014), andCai et al. (2015)). Global warming can also profoundly affect P DO, shortening its lifespan and suppressing its amplitude (Li et al. (2020)). Moreover, the melting of land ice and warming ocean waters cause rising sea levels affecting coastal shorelines. High-tide flooding is increasing in magnitude and frequency: minor floods occur multiple times per year; major floods might occur even yearly. Even if GHG emissions stopped, the sea level would continue to rise. Finally, the Arctic is warming twice as faster as the planet on average, and finding irreversibility in N H but not in SH, is interesting in this respect. Despite not being conclusive, the results might indicate that some irreversible environmental changes are ongoing. This paper links the concept of an environmental tipping point to the statistical concept of time irreversibility. A tipping point signals an environmental change that is large, abrupt, and irreversible and generates cascading effects. A tipping point is a point of no return, which we associate with a temporal asymmetry in a phenomenon's probabilistic structure, whereby it behaves differently according to the direction of time considered. This univocity along the time direction signals that the system has undergone an irreversible change. Well-known tipping points concern the Greenland and the West Antarctic ice sheets, the Atlantic Meridional Overturning Circulation (AM OC), thawing permafrost, EN SO, and the Amazon rainforest. Recent IP CC assessments suggest that tipping points might occur even between 1 • C and 2 • C warming relative to pre-industrial temperature averages. Therefore, they are likely to arise at current emissions levels if they have not already occurred.
We then introduce two new strategies, grounded on mixed causal and noncausal models, to detect whether a stochastic process is time-reversible (TR). Unlike existing approaches, our methods do not impose strong restrictions on the model and are straightforward to implement. Moreover, similarly to Proietti (2020), they can also be applied to non-stationary processes and, therefore, useful to assess some key variables, such as temperature anomalies and GHG emissions, which appear to exhibit this property. Our simulation studies show that the strategies perform accurately and have a solid ability to detect TR. In the empirical analysis, we have considered fourteen climate time series, i.e., annual and monthly global temperature anomalies (GLO, GL, GO; GCAG, GIST EM P ), solar activity (SA), natural oscillations (N AO, SOI, P DO), the global mean sea level (GM SL), the Northern (N H) and Southern (SH) Hemisphere sea ice areas, global sea levels, greenhouse gas emissions (GHG, N 2O). We detect time irreversibility in GHG and N 2O emissions, SOI and P DO, GM LS, N H, and the monthly temperature anomaly series. Yet not in the annual temperature series, SH, and N AO (and SA). The time irreversibility of GHG emissions is a noticeable property of variables that are well-known causes of global warming. It may then explain the time irreversibility of GM SL, N H, global temperature, and some natural oscillation indices such as P DO and SOI, and therefore signal that some potentially irreversible environmental changes are ongoing. This evidence might not be apparent from annual temperature data due to the relatively smaller sample size available for annual than monthly data. Recent studies suggest global temperature has already warmed by 1.3 • C and could cross the 1.5 • C threshold within a decade. While not conclusive, our findings urge the implementation of correction policies to avoid the worst consequences of climate change and not miss the opportunity window, which might still be available, despite closing quickly. with which time irreversibility is detected when the error term has a symmetric Student's-t distribution (γ=1) and ν 0 = 3. Finally, r and s are assumed as unknown and p as known. Y : M AR(1, 1); φ 0 = 0.8, ϕ 0 = 0.1; ν 0 = 3, γ = 1 Y : M AR(1, 0); φ 0 = 0.8; ν 0 = 3, γ = 1

with which time irreversibility is detected on non-stationary time series; r and s are assumed as unknown and p as known. Data are considered to have quarterly frequency (λ = 1600). Chen et al. (2000) state that to jointly test X t,k for a collection of k values, a portmanteau test is required.