Multi-scale CO2

Five time series of estimated atmospheric CO 2 with sampling intervals ranging from 0.5 million years to the relatively high frequency of one week are analysed. The yearly series shows a clear increasing trend since the beginning of the ﬁrst Industrial Revolution around 1760. The weekly series shows a clear increasing trend and also seasonal variation. In both cases, the trend is ﬁtted by a conceptual model that consists of a baseline value with an exponential trend superimposed. For the weekly series, the seasonal variation is modelled as an exponential of a sum of sine and cosine terms. The deviations from these deterministic models are treated as detrended and deseasonalised time series. Then, three sub-categories of autoregressive integrated moving average (ARIMA) models are ﬁtted to the ﬁve time series: ARMA models which are stationary; FARIMA models which are stationary but have long memory and are fractal processes, and ARIMA models which are variations on a random walk and so non-stationary in the variance. The FARIMA and ARIMA models provide better ﬁts to the data than the corresponding ARMA models. All the ﬁtted models are close to the boundary of stability, and are consistent with claims that climate change due to an increase in atmospheric CO 2 may not quickly be reversed even if CO 2 emissions are stopped. autoregressive term is non-signiﬁcant. The third autoregressive term can be attributed to approximating derivatives with centred differences, and both models allow a reasonable physically interpretation of the process as a damped oscillator, which is modelled in terms of displacement and its ﬁrst and second derivatives, with random inputs. The interpretation is of a borderline stable second order dynamic system with long memory.

recent century, and it is potentially the most difficult to avoid [18]. Global warning not only stresses ecosystems through various mechanisms, but also has impacts other areas which can lead to positive feedback. For example, Australia's bushfires, intensified through global warming, are contributing to one of the biggest increases in the concentration of atmospheric CO 2 in the past 60 years [9]. [20] claim that the climate change that takes place due to increases in carbon dioxide concentration is largely irreversible for 1,000 years after emissions stop. Cessation of emission of carbon dioxide will decrease radiative forcing, but this decrease is largely offset by slower loss of heat to the ocean. A consequence is that global atmospheric temperatures are unlikely to drop significantly for at least 1,000 years.
We fit time series models to records of atmospheric CO 2 at very different time scales, starting with a record based on pedogenic carbonate assessments extending back 420 million years, with a sampling interval of 0.5 million years [8]. At an intermediate level we consider ice core records extending back around a million years with sampling intervals of the order of 1,000 years. These are compared with relatively high frequency data from the past 2,000 and last 50 years sampled at yearly, and weekly intervals respectively. The objective is to compare the time series models and look for unifying features, with particular emphasis on stationarity. The results are discussed in the context of Solomon et al's (2009) claim that increases in carbon dioxide are largely irreversible.

Data
The main objective of this study is to compare time series models that are fitted to estimates of atmospheric CO 2 made over quite different sampling intervals, and in the case of weekly data direct measurements of CO 2 . The time intervals are: 0.5 million years; around one thousand years, yearly; and weekly (Table 1).  Foster et al. (2017) compiled the dataset using five independent techniques applied to data from 112 published studies covering the last 420 million years. In order to ensure the highest-quality compilation, they applied a set of criteria that left 1,241 observations to be analysed. Monte Carlo resampling and a LOESS fit were used to interpolate the series to a regular 0.5 million years interval. The data is available from the supplementary information from https://www.nature.com/articles/ncomms14845, with DOI: 10.1038/ncomms14845.
The Vostok ice core dataset was obtained from the National Oceanic and Atmospheric Administration (NOAA) on 10th of March 2020. It can be accessed via https://www. ncdc.noaa.gov/paleo-search/study/17975. This original publication is "Antarctic Ice Cores Revised 800KYr CO2 Data", with data set ID noaa-icecore-17975.
The original publication of Law Dome dataset is "Law Dome Ice Core 2000-Year CO 2 , CH 4 , and N 2 O Data", which was contributed by David Etheridge (Original receipt by WDC Paleo) and last updated on July 2010. This dataset was acquired from the National Oceanic and Atmospheric Administration (NOAA) on 25th of February 2020, it can be accessed via https://www.ncdc.noaa.gov/paleo-search/study/9959, with data set ID noaa-icecore-9959. This dataset has been split into two time intervals, 1 yr AD to 1760 yr AD, and 1761 yr AD to 2004 yr AD. This is because the dramatic increase due to the first Industrial Revolution, which started in 1760.
The weekly time series was acquired from the National Oceanic and Atmospheric Administration (NOAA) on 6th of March 2020, and it can be accessed via https://www.co2. earth/co2-datasets. These data are weekly carbon dioxide observations, measured as mole fraction in dry air, from the Mauna Loa Observatory in Hawaii since May 1974.

Filtering deterministic models
There is a clear increasing trend in CO 2 , that appears to be approximately exponential, since the beginning of the first Industrial Revolution, around 1760. The trend is estimated with a conceptual model that consists of an exponential increase in CO 2 superimposed on some baseline value (a), y t = a + be ct .
The model is fitted by non-linear least squares, using the nls function in R, and the residuals are considered to be the detrended time series. There is also a marked seasonal pattern in the time series of weekly data from Mauna Loa. The seasonal pattern is modelled by modulating the trend with an exponential of a linear sum of cosine and sine functions, where C m = cos(m2πt/p), S m = sin(m2πt/p) for m = 1, . . . , M , The value of p is equal to the seasonal period, 365.25/7 in the case of weekly data, and M is selected as the smallest value such that neither of the coefficients of C M +1 nor S M +1 is statistically significant at a 0.05 level. The residuals (observations less fitted values) from the deterministic models are considered to be detrended and deseasonalised time series. There is no evidence of any deterministic trend in the time series that pre-date the first Industrial Revolution, so they are assumed to be realisations of models that are stationary in the mean.
We then calculated autocorrelations, partial autocorrelations, and spectra, and applied various unit root tests that aim to detect non-stationarity in variance before fitting ARIMA and FARIMA models. We review these methods briefly in the following sections.

Akaike Information Criterion (AIC)
The Akaike information criterion (AIC) [1] is the most commonly used criterion for model selection when maximum likelihood is used to estimate parameters. The model with a lowest AIC is considered as the best fitting model. However, AIC tends to favour models with quite large numbers of parameters, and we are looking for a succinct model to emulate the dynamics of the processes. [4] suggest that a simpler model may be preferable to a model with more parameters if the simpler model has an associated increase in AIC that does not exceed two. But, this advice can be rather restrictive if applied to long time series because a small reduction, of no practical significance, in the estimated standard deviation of the errors can be statistically significant. In general, a sensible choice of model, or comparison of models, depends on the context, and it useful to consider both AIC and the estimated standard deviation of the errors (based on the unbiased estimator of the variance of the errors).

ARIMA model
An AR(p) model is a linear regression of the current value of a variable, y t say, on p past values. This can be generalised to an ARMA(p, q) model which includes a linear combination of q past errors. An ARMA(p, q) model is strictly stationary, and hence stable, if all the roots of the polynomial lie outside the unit circle. In particular AR(1) is stationary for |φ 1 | < 1, and AR(2) is stationary for |φ 2 | < 1 and φ 2 ± φ 1 < 1. The ARIMA(p, d, q) model allows for a nonstationary series to be differenced d times before fitting an ARMA(p, q) model. So, the general definition is as follows: where B is the backward shift operator, defined by By t = y t−1 .
Typically a range of ARIMA models can provide plausible fits to data from a dynamic system. The AIC is one option for choosing between them, but it is entirely empirical and tends to favour models with a large number of parameters. The AIC may be appropriate for forecasting applications, but we prefer to focus on models with some conceptual basis for modelling CO 2 . The differencing could be used to remove a deterministic trend, but in the context of a process that is stationary in the mean a d of 1 allows for an unstable process. The crucial feature of an unstable process is that it does not return to some fixed level, a well known example being the random walk. In mathematical terms, the defining characteristic of an unstable process is that the polynomial in B acting on y t has a unit root. Unless strongly indicated by model fit, we prefer AR(1) and AR(2) models, without moving average terms, because they have a conceptual interpretation in terms of elementary linear dynamic systems. The AR(2) process is a discrete time version of the second order differential equation that describes a linear oscillator subject to a sequence of random impulses. The derivatives are approximated by differences. The canonical example of a linear oscillator is a mass suspended by a spring with some damping arrangement attached (e.g. Thomson 1993, Theory of Vibration With Applications (4e), Prentice Hall)

FARIMA model
Some time series exhibit marked correlations at high lags, and they are referred to as longmemory processes [3]. The time series is considered to have long memory if the acf decays more slowly than an exponential decay. This slower decay is usually modelled as a powerlike decay [10]. The FARIMA model allows for fractional differencing, generalising the ARIMA model's integer order of differencing to allow the d parameter to take on fractional values, −0.5 < d < 0.5 [2]. The range 0 < d < 0.5 gives long memory processes, and it lies between a stationary AR(1) model and a non-stationary random walk [3]. We are interested in this model because it is close to being a random walk yet still stays bounded, and returns to a mean value if there are no error inputs. Formally, the (1 − B) d term in the ARIMA model is expanded with the generalised binomial expansion: as for as B L , where L is some cut-off point typically around 30.

Random Walk
The ARIMA(0,1,0) model is a random walk. It generally defined as follows: where w t is discrete white noise with mean 0 and variance σ 2 w .
Using the backward shift operator B, we can rewrite the defining equation as It is formally stationary in the mean, but non-stationary in the variance and so unbounded. The current value is the best predictor of future values and the model has no tendency to return to the initial value.

Augmented Dickey-Fuller (ADF) test
The model for the Augmented Dickey-Fuller test is as follows: where △ = (1 − B)y t is the first difference operator, and in the R implementation K is defined as trunc((length(y) − 1) (1/3) ). The summation term in the equation above is to correct for the presence of serial correlation.
The null hypothesis of the Augmented Dickey-Fuller test is the presence of unit root, that is H o : β = 0 against the alternative which is H a : β < 0. Under the null hypothesis, y t itself is non-stationary. Hence the ordinary central limit theorem does not apply to the least squares estimator,β, of β, and its sampling distribution does not tend to a t-distribution as the length of the time series increases. Dickey and Fuller have tabulated the asymptotic distribution of the least squares estimator for β under the null hypothesis of it being a unit root. It turns out that we can compare a t-statistic with the values of this Dickey-Fuller distribution. If the t-statistic less than the critical value, then we reject the null hypothesis [7].
There are two functions in R that can preform Augmented Dickey-Fuller test: adf.test and ur.df. The difference between the two is that for adf.test the critical value is based on the model with drift (intercept) term, whereas ur.df based on the model without drift term. Moreover, the ur.df does allow us to select the lag according to AIC or BIC.

Phillips-Perron (PP) test
The model is as follows: The null hypothesis of this test is that α = 1, against the alternative which is α < 1. They used a corrected form of t-test in order to correct for the presence of serial correlation and heteroscedasticity in the error term [21]. Again, there are two functions in R that can preform the PP test: pp.test and ur.pp.

Kwiatkowski-Phillips-Schmidt-Shin (KPSS) tests
Let y t be the observed time series for which we wish to test stationarity. We assume that we can decompose the series into the sum of a deterministic trend, a random walk, and a stationary error [11]: Here γ t is a random walk: where the u t are i.i.d. N(0, σ 2 u ). The initial value γ 0 is treated as fixed and serves the role of an intercept. The stationarity hypothesis is simply σ 2 u = 0. Since ǫ t is assumed to be stationary, under the null hypothesis y t is trend-stationary (Tang et al., 2019). The statistic they used is both the one-sided Lagrange Multiplier (LM) statistic and the Local Best Invariant (LBI) test statistic for the hypothesis σ 2 u = 0, under the stronger assumptions that the u t are normal and that the ǫ t are i.i.d. N(0, σ 2 ǫ ). Because the parameter value specified by the null hypothesis is on the boundary of the parameter space, they were interested in a one-sided LM test rather than a two-sided test [11].
Let e t , t = 1,2,..., T , be the residuals from the regression of y on an intercept and time trend. Letσ 2 ǫ be the estimate of the error variance from this regression (the sum of squared residuals, divided by T). Define the partial sum process of the residuals: Then the LM (and LBI) statistic is and high values are taken as evidence against the null hypothesis. The function kpss.test and ur.kpss in R can be used.
In summary, the ADF and PP test a null hypothesis of a unit root, which is non-stationary in variance. Conversely, the KPSS test a null hypothesis of stationarity.

Analysis
We fit AR(p), ARIMA(p, d, q) and FARIMA(p, d, q) where both p and q are no greater than three, to all five time series. These five time series include million-year series; thousand-year series; pre Industrial Revolution yearly series; detrended post Industrial Revolution yearly series; and detrended and deseasonalised weekly series. To eliminate or reduce the trend and seasonal patterns, we fit exponential trends above a threshold to the time series after the Industrial Revolution, and a seasonal model to the weekly series, prior to model fitting. The residuals, defined as observations less fitted values, are considered to be a realisation of a process that is stationary in the mean.

Million-year time scales
The time series [8] is based on pedogenic carbonate assessments extending back 420 million years at 0. At the 0.5 million year scale there are two peak values above 2000 ppm, around 400 and 200 million years ago. Current levels are around 400 ppm, but this is above the minimum value of 220 ppm that occurred around 300 million years ago. For comparison, the values around 300 and 80 million years ago are close to the value before the first Industrial Revolution (around 1760s), which is 270 ppm. The correlogram (Figure 2(a)) shows the series is highly auto-correlated, and the pacf (Figure 2(b)) shows the significant lags at 1 and 2. These suggest the series can be modelled reasonably as AR (2). But given the dominant autocorrelation at lag 1, an AR(1) model, with a coefficient close to 1, may be a reasonable approximation. Both periodogram and smoothed spectrum show that the time series is a realisation of a low frequency process. If the unit root tests are applied, neither ADF nor PP provide evidence against a unit root at the 0.05 level of significance and KPSS provides evidence against a null hypothesis of stationarity. Hence, all of those three tests consistent with a hypothesis of a unit root. Furthermore, fitting a simple AR(1) model returns a coefficient of 0.9939, which is very close to 1. However, a random walk model is somewhat unsatisfactory for a environmental process that is physically bounded, and a stationary fractional differenced model is conceptionally preferable. The results of fitting models are shown in Table 2. The residuals from AR(1) and AR(2) models still show significant autocorrelations (Figure 3(a), 3(b)). Based on Akaike information criterion, the ARIMA(3,1,2) would be chosen as best fitting model. A consequence of the differencing parameter, d=1, is that the process is unstable and so unbounded. Instability within the range of CO 2 in the time series record is quite plausible, and thresholds could be introduced to give a bounded non-linear variant to the model. The FARIMA models provide a better fit than the AR(p) models, and as they are equivalent to a high order AR model they are also stationary. Even though FARIMA(3, 0.330, 0) gives a lower AIC, as well as lower residual variance than FARIMA(2, 0.241, 0), the third autoregressive term is non-significant. The third autoregressive term can be attributed to approximating derivatives with centred differences, and both models allow a reasonable physically interpretation of the process as a damped oscillator, which is modelled in terms of displacement and its first and second derivatives, with random inputs. The interpretation is of a borderline stable second order dynamic system with long memory.   Although the difference in AICs between FARIMA(2, 0.241, 0) and ARIMA(3, 1, 2) is nearly 20, this can be contributed to the long time series, which allows identification of statistically significant terms which do not correspond to a substantial reduction in the variance of the residuals. The corresponding estimated standard deviations of the errors are 39.68 and 39.37 respectively.
The coefficients of the ARMA(2,2) are very close to the boundary of the stability region, that is |φ 2 | = 0.889, and φ 2 + φ 1 = 0.997. In contrast, if an unstable ARIMA(2,1,0) is fitted then the coefficients of the model for differences are far removed from the boundary of stability.

Thousand-year time scales
The Vostok CO 2 dataset has differing time intervals between observations ( Figure 4) and the time scale is from past to present. The usual approach to analysing time series with unequal time intervals is to interpolate to obtain equal time spacing. We considered linear interpolation and cubic spline interpolation to a time interval of 1,176 years corresponding to 372 points, which exactly equals the number of points in the original time series . The variogram describes the covariance structure over space, or time, as a function of distance h between points, and is defined as for all locations u, where Z(s) is the value of the variable of interest at location s.
Estimation of the variogram does not assume equal spacing, so variograms of the original time series, the linear interpolated series, and the spline interpolated series are compared in Figure 5.
There is reasonable agreement between the original time series and the interpolated series. The difference between the variograms for the interpolated series is not statistically significant. The linear interpolation is slightly closer to the original time series overall, although the spline interpolation is closer at the beginning and the end. As a further check of the interpolations we investigated frequency compositions. For all three time series we fitted regressions with cosine and sine terms with periods at multiples of one cycle per record length. That is regressions with 2M predictor terms of the form cos(m2πt/n) and sin (m2πt/n) for m from 1 up to M , and plotted the coefficients of determination (R-square) against M . The result is shown in Figure 6. The spline interpolation increases slowly relative to the other two series. The interpolated series have identical plots as we increase the number of harmonic terms, the cumulative periodograms, and are quite close to the regression results for the original series (not a precise periodogram because the orthogonality property is lost if the time intervals vary).
The conclusion is that both interpolation methods gave evenly spaced time series with a similar auto-correlation structure to the original time series. We chose to continue with the linear interpolation, because it is simpler and avoids a few unrealistic values at the beginning of the series where there are some relatively long time intervals between observations. The correlogram (Figure 8(a)) shows a linear decay which suggests a non-stationary time series model may be appropriate. Whilst non-stationary models might be considered physically unrealistic for the process over 420 million years, they can provide reasonable local approximation over sub-periods of, here, four hundred thousand years. The partial correlogram (Figure 8(b)) shows significant lags at lag 1 and lag 2, which suggests the time series model has two autoregressive terms. The spectrum indicates a low frequency process which is indicative of a process with a unit root or a fractionally differenced process. The results of fitting various models to the carbon dioxide time series with an interval of 1,176 years are shown in Table 3 below. The FARIMA(2, 0.064, 0) process is not statistically significantly worse than the ARIMA(3, 1, 2), and in practical terms they provide equally good fits [4]. The FARIMA model is preferred because it is simpler and it is stationary. It is also consistent with the FARIMA model fitted to the 420 million years series and a fractal process. We also investigated directionality in the thousand-year time series [13], especially whether there is a difference between the time from a threshold to a peak and the time from the peak back to the threshold. The threshold was taken as the 0.80 quantile of the marginal distribution of the time series. A peak was defined as the maximum value of the time series over the period from an upwards crossing of the threshold to the next downwards crossing. The statistic is calculated from all peaks for which there are four time steps following the up crossing and before the down crossing. The mean of the variable over the four time steps following the peak is subtracted from the mean over the four time steps before the peak. These differences are assumed to be independent and are analysed by a one-sample t-test for paired comparisons. The value of the test statistic is -2.07. So there is evidence that the rise to a peak is more rapid than the decline following it. This provides some support for Solomon et al's (2009) claim that climate change due to an increase in atmospheric CO 2 may not quickly reversed even if CO 2 emissions are ceased.

Yearly ice core data
The concentrations of CO 2 were inferred by using spline fits to the Law Dome firn [15], ice core records and the Cape Grim record (Figure 9). There is a dramatic increase from around 1800 yr AD, that can be attributed to the Industrial Revolution. We analyse the series in two parts: before the Industrial Revolution (1 yr AD -1760 yr AD) and after the Industrial Revolution (1761 yr AD -2004 yr AD). There is noticeable variation in the pre-Industrial Revolution series (Figure 10(a)), but the range for CO 2 is relatively small. In contrast, there seems to be an exponential increase for the series after 1760s (Figure 10(b)). The ADF, PP and KPSS tests in R are applied to the pre-Industrial Revolution series only, as there is a clear trend in the post-Industrial Revolution series. There is no evidence to reject the null hypothesis that there is a unit root at the 0.05 significance level, based on the ADF and PP tests. Furthermore, the KPSS test shows evidence to reject the null hypothesis of the stationarity of the series.

Pre-Industrial Revolution (IR) period
The correlogram decays very slowly down from unity, which is one of the characteristics of a random walk. The partial correlogram shows a positive partial autocorrelation at lag 1, and negative autocorrelations from lag 2 up to lag 5, and this might suggest a moving average term. The spectrum shows the property of low frequency and no dominant peak.  Table 4 shows the results for different models. There are no standard errors for the coefficients of FARIMA(1, 0.5, 0) because R is unable to compute correlation matrices for those two models. The autoregressive coefficient returned by AR(1) model suggest our series resembles a realisation of a random walk, this conclusion is consistent with the differencing operator d returned by the FARIMA(2, 0.499, 0). The ARIMA(3,1,2) again is the model that has the lowest AIC value, and although the unbounded property, a consequence of its being unstable, is not physically reasonable over a long period, it is an adequate approximation over a relatively short period. The FARIMA(2, 0.499, 0) is a succinct and more physically realistic model.
The appearance of stochastic trends, that are modelled by a random walk, may be accounted for as segments of a relatively low frequency oscillation in the underlying continuous process. The FARIMA(p,d,0) models are consistent with the underlying continuous process being fractal. Table 4: Parameter estimates of AR, ARIMA, FARIMA models for the pre-IR data.

Post-Industrial Revolution (IR) period
There is a clear exponential growth in atmospheric carbon dioxide concentration after the first Industrial Revolution, and this is fitted using the conceptual model of Equation 1. A time series plot of the residuals which are the detrended series is shown in Figure 12. Even after removing a trend, the ADF and PP unit root tests do not provide evidence to reject the null hypothesis of a unit root at the 0.05 significance level. Furthermore, the KPSS test provides evidence to reject a null hypothesis of stationarity at the 0.1 significance level.
The correlogram shows a linear decay which also suggests non-stationarity (Figure 13(a)). The partial correlogram is significant at lag 1, suggesting an AR(1) model. Based on the results of three unit root tests, we might conclude that the series is a realisation of a random walk. Also, the estimated coefficient returned by the AR(1) model is 0.9988, which is very close to 1 and consistent with the unit root tests. Both the raw periodogram and smoothed spectrum show the dominance of low frequency variation, again consistent with a random walk.    Table 5 gives the results of fitting models for the post-IR yearly time series. We observe that for this series, ARIMA (3,1,2), is again the model that has the lowest AIC. However, several other models have standard deviations of residuals that are close to 0.082. In particular, AR(2) with a standard deviation of residuals of 0.087. The estimated values of the coefficients of the AR(2) are close to the boundary of stability, and so realisations from the fitted AR(2) trend to be more oscillatory than the observed time series. The fitted FARIMA(1, 0.497, 0) model has a somewhat higher standard deviation of residuals at 0.105, but realisations from the fitted FARIMA(1, 0.497, 0) show stochastic trends and are generally qualitatively more similar to the observed series. Again, the fitted FARIMA(1, 0.497, 0) is close to the boundary of stability, and it is a plausible conceptual model that is consistent with a fractal process.

Weekly Mauna Loa Data
These data provide weekly carbon dioxide that measured as mole fraction in dry air (direct measure) since May of 1974, in Mauna Loa Observatory, Hawaii ( Figure 14). It is acceptable to use Mauna Loa as a proxy for global CO 2 levels. This is because CO 2 mixes well throughout the atmosphere. This is supported by the statistically non-significant difference in trend in Mauna Loa CO 2 levels (1.64 ppm per year) and global CO 2 levels (1.66 ppm per year) [5]. There are 20 observations that have values equal to -999.99, which indicates they are not available (NA). The reasons for them to be unavailable have not been stated. These missing values are imputed using linear interpolation, which is easily implemented as our dataset is equally spaced.
From Figure 14, we can observe the seasonal patterns, and much of this can be explained by the role of plants in the carbon cycle [6]. Plants use CO 2 as well as sunlight and water, in order to make food or other substances that they need to grow, and they release O 2 as a by-product. This process is called photosynthesis, and this can draw down the CO 2 [16]. Another process that is part of the carbon cycle is called respiration, this is the process which plants decompose and produce CO 2 .
As there is an apparent trend as well as periodicity, hence the first thing that we need to do before any analysis is to remove the trend and seasonality. This can be done by fitting the model of Section 3.1 with M = 4 (Equation (2)). There is no clear trend or seasonal pattern of the residuals after detrended and deseasonalised ( Figure 15). The correlogram (Figure 16(a)) shows highly significancy between the observations, and the partial correlogram (Figure 16(b)) shows significancy at lag 1, 2, 3 and 4. The partial correlations for lag 5, 6, and 7 are relatively small. Both the periodogram and smoothed spectrum (Figure 16(c) -(d)) show the property of low frequency.
The ADF and PP tests on the residuals show the evidence to reject the null hypothesis of the presence of unit root, at the 0.05 significance level. However, the result returned by the KPSS test is contradicted with the ADF and PP tests, it shows the evidence to reject the null hypothesis of stationarity. If an AR(1) model is fitted to the residuals, the estimated coefficient should close to 1, that is 0.821.
Finally, in Table 6 below, we looked at the weekly data. The ARIMA(2,0,1), which is equivalent to the ARMA(2,1), has the lowest AIC with standard deviation to be 0.401. Even though the difference in AICs between FARIMA(2, 0.497, 0) and ARMA(2, 1) is approximately 30, this can be contributed to the long time series, which allows identification of statistically significant terms which do not correspond to a substantial reduction in the variance of the residuals. They have standard deviation of residuals to be 0.403 and 0.401, respectively. Moreover, the FARIMA(2, 0.497, 0) is more conceptually reasonable, and it also gives a reasonable physically interpretation. Table 6: Parameter estimates of AR, ARIMA, FARIMA models for the weekly data.  The preferred models in Table 7 are not necessarily those with lowest AIC, but are those we judge to be the most suitable given a preference for relatively simple conceptual models. In particular, autoregressive models without moving average terms, AR(p), have a conceptual interpretation as a discrete approximation to a dynamic system driven by a sequence of independent shocks, are special cases of multiple regression and so relatively simple, and the pacf can provide a clear indication of a suitable order. All autoregressive models fitted to the time series were stationary and so stable, but close to the stability limit. Allowing moving average terms, ARMA(p,q), did generally improve the fit in terms of AIC with the estimates of autoregressive coefficients remaining close to the stability limit. The values of coefficients of moving average terms do not affect stability. Realisations of time series from stable AR(p), or ARMA(p,q), models that are close to the stability boundary will eventually return to some mean value if inputs are set to zero, but return to the mean will be slow. In contrast, ARIMA(p,1,q) models have a unit root and are not second order stationary, and are unstable because realisations do not return to a mean value if inputs are set to zero. This qualitative difference between ARMA(p,q) and ARIMA(p,d,q) models can be critical, for example in the context of control systems, but in the context of modelling CO 2 concentration all the models are consistent with claims that anthropogenic emissions will not easily be reversed. The lack of any underlying mean for ARIMA(p,d,q) models is conceptually appealing, but a less reasonable corollary is that they are unbounded. However, this limitation can be removed by introducing thresholds above, and below, which the process is modelled as stable. The FARIMA(p,d,0) models come between the ARMA(p,q) and ARIMA(p,1,q) models, but are still stable. The FARIMA(p,d,0) models account for a slow decay in the acf and indicate a fractal process. The notion of a fractal process fits nicely with the fact that weekly, and even yearly, sampling of an underlying continuous process with much lower natural frequency is likely to identify low frequency variation as a trend.
In Table 7, a preferred model including unit root is only shown if it has lower estimated standard deviation of errors than the model without unit root. Also, the models after IR are for errors after the conceptual trend model, with seasonal terms in the case of weekly data, is fitted. In summary, both the FARIMA and ARIMA models provide a better fit to the time series than ARMA models with the same p and q. Up until the first Industrial Revolution there is little to choose between them although the ARIMA have slightly lower estimated standard deviations of errors at two of the three scales. Formally, there is no evidence against a hypothesis of a unit root and there is evidence against a hypothesis of stationarity at the 0.1 level of significance, for all except the weekly series. These findings support a choice of ARIMA models over the corresponding ARMA models. But unit root tests do not allow for moving average terms, and ARMA models with q = 2 are close to the corresponding ARIMA models in term of AIC. However, moving average terms do not have much conceptual interpretation in this context.
The FARIMA models suggest fractal processes. The ARIMA models are not stable, which is not physically plausible at the 0.5 million years scale but may be a realistic approximation over the shorter periods. The best fitting ARMA models are close to the boundary for stability. Whichever models are chosen, the borderline stability is consistent with a general claim that a reduction in atmospheric CO 2 following a cessation of CO 2 emissions will be a protracted process. In particular, the evidence that the rise to a peak is more rapid than the decline from it (found in the thousand-year time interval series) provides some support for the claim by [20].
Following the Industrial Revolution the anthropomorphic effects appear to dominate the natural sequence of events. After fitting conceptual models for the trends, the prediction errors are well modelled as either FARIMA or ARIMA models. The former are preferred as superimposing unstable models on an exponential increase is not a convincing conceptual model. It is interesting that the FARIMA models still provide a reasonable fit as this is consistent with a fractal process at scales from half a million years down to one week.
Except for the weekly data from Mauna Loa, we have modelled time series of estimates of CO 2 rather than direct measurements. In all cases the CO 2 concentration is, nominally at least, aggregated over the sampling interval. It is possible that autocorrelations at short lags may be affected by the method of estimation. But, this cannot account for the evidence of a lack of stability in the underlying process. The borderline stability has serious consequences because a unit root corresponds to a process that does not return to some average value. If this is realistic, ceasing CO 2 emissions will not necessary lead to a return to lower levels but it will help prevent levels increasing. Empirical evidence based only on time series models might not be convincing in itself, but when it augments physical explanation (e.g. Solomon et al., 2009) the forecast should be taken seriously.
This study was funded by the Australian Research Council for a scholarship for the Master of Philosophy.

Conflicts of interest
Not applicable.
Availability of data and material The Australian Bureau of Meteorology and the National Oceanic and Atmospheric Administration provided data for analysis. Glonek. All authors provided critical feedback and helped shape the research, analysis and manuscript.