Autocorrelation and trends

by Richard Saumarez

A recent debate at Climate Etc. has involved autocorrelation and trends.   Essentially, the argument is that the temperature signal is perturbed by a random signal, and the influence of this random fluctuation persists for an unknown length of time and so the temperature at any particular instant is influenced by its past history.

This is a presentation of some of the issues that were raised by the papers recently presented by Ludecke, Link and Ewert [see here, here, here and here].

Why does persistence produce trends?

It may not be intuitively obvious why trends should occur in persistence models.  Figure 1 shows a resting point (stage 0) at zero.  A random, normally distributed value is added to it and it moves to point 1.

Figure 1 An Illustration of how persistence leads to trends.  Upper trace shows how persistence stops alters the probability of the temperature returning to baseline, creating a trend.   The bottom trace shows how a lack of persistence does not bias the probability of the next step.

The value at point decays slowly, it has persistence, so that after a time interval it reaches point 1a.  A second random value is then added to it.  Notice that probability of it returning to zero or below is the area of the distribution curve that is below zero, i.e.: very small.  This illustrates why trends occur with persistence.

As the value is increased, note that the rate of decay increases – compare the decay rate between points1 and 1a with that between 5 and 5a. Therefore the observed value will remain above zero until there has been a run of negative disturbances that will return it.

The lower trace shows a non-persistence model.  In this case, the value has returned to zero before the next disturbance and the observed values will have a mean of zero.

This is the basis of the argument that “persistence” will cause trends in temperature.  If an impulse of radiation is stored as heat, it carries over so that at the next impulse, the temperature will be slightly raised and become raised still further, resulting in a temperature trend.  The controversy is how significant is this effect and can observed temperature trends be explained on the basis of a simple heat reservoir and a random change in radiative fluxes?

What models of Heat storage or Persistence can we use?

Any simple model begs the question of what mechanism makes temperature persistence conform to a particular mathematical relationship.

There are two basic models, one is a power law model in which the temperature decays by a power of time and is often called the Hurst law, who defined in terms of Hydrology.

T = (Ht)^g

(JC formatting note:  the g superscript should be gamma).

(t is time, H is a scaling parameter and g, the exponent, describes how the process varies with time – g is negative).

Using this idea, Figure 2 shows a model, consisting of a leaky reservoir, which is surrounded by a dam that is more porous at the top than at the bottom.  If it is filled suddenly with water at a high temperature, water will escape rapidly from the top and as the level subsides, there will be a slower leak through the bottom levels.  Thus the level will decay at a slower rate as the reservoir empties, not only because the pressure head is smaller, but also because the rate of transfer is reduced by passing through a smaller area of less porous material.

Figure 2 A model of a power law system with the response shown (black). The responses for g=0.6-0.8 are shown in blue and the responses to the reservoir being filled to different levels are shown in red.  

I have solved this model to give the decay curve and corresponds to a gamma of 0.7.  It is important that the response to an impulse of water entering the reservoir will decay differently depending on how full the reservoir is (red curves) and therefore this system is non-linear. When computing a response, the system can be treated as quasi-linear by assuming that the actual level of the water does not vary much from its mean level.  This allows linear methods of processing to be applied.  If however, there is a wide range of levels, because the range of inputs is large compared to volume of the system, the behaviour of the system must be solved explicitly through its governing differential equations.  Whether one can really treat the temperature system as quasi-linear, or not, depends on measuring the inputs (if we knew what they were) and the responses of a series of temperature records and determining the coherence between them.

The alternative model is a multi-compartment model that has exponentially decaying compartments as shown in figure 3.  It is well known that temperature persistence does not conform to a simple, single exponential decay, but a system with a fast  component and a large volume slow component can produce much longer levels of persistence.  As shown in figure 3, this model produces a decay that is very similar to a power law decay and, if communication between the compartments can occur, the decay curve becomes almost indistinguishable from a power law model.

Figure 3 The response of a multi-compartment model.  The power law responses for g=0.7 and 0.8 are shown in blue. Note that a 3 compartment model has a very similar response to a power law model.

There is a very important assumption in the use of persistence models to calculate global temperature trends.

The Earth is assumed to be driven by a random change in flux and a component of this affects all parts of the Earth equally.  In other words, some of the “input” at each measuring point is correlated.  If they were independent, global trends would not occur.

This assumption allows one to assume an “average model” for persistence and to drive it in order to produce a simulated temperature record. 

Determination of the model and its relationship to trends

This was done in a number of steps:

1)    Temperature records were extracted from the HADCRUT III, UTAH and GHCNM databases and all those at least 100 years long with less than 5% of the data missing were selected.  The missing data was interpolated linearly across the gaps.

2)    The yearly (“seasonal”)  fluctuations were removed.

3)    This data was used to generate an autocorrelation and signal power model.

4)    A system response was created using either a bi-exponential or a Hurst model.

5)    Gaussian noise was used to create a sequence 1364 years long that have the same variance as the residual signals.  A one hundred year record was randomly selected from each sequence, in order to account for measured temperature trends to start at a random time with respect to any trend in temperature. The amplitudes and duration of the trends were extracted from 100000 sequences.

Figure 4

 Autocorrelation functions

The ACF is a measure of how a sample in the record is correlated to its neighbouring samples and gives a measure of the persistence as it demonstrates the effects of past samples in determining the current sample.  In principle, a copy of the signal is made and is shifted in time by one sample. The correlation coefficient between the shifted copy and the original version determines the relationship between the signal and that at one sample apart.  This is repeated for all possible shifts of the copy of signal and the correlation coefficient is plotted against the value of the shift (Figure 5).  In practice this can be done efficiently through the discrete Fourier Transform of the signal.

The long-term behaviour of a system can only be estimated from about 1/5th of the length of the record

There are a number of problems in calculating the ACF.  The first is that as the signal copy is shifted and each sample multiplied with the sample at the same position in the original signal, the samples at the end of the shifted signal “drop off” the end and so do not contribute to the correlation at that shift, making the estimate less reliable.  Therefore only shifts of about 1/5^th of the total signal length should be used, in our case about 20 years.

The method used by L and LLE is slightly different from the autocorrelation function, in that they calculate a cumulative sum of the temperature excursion and manipulate this sum as a function of time by fitting a second order polynomial as a smoothing and de-trending step.  They then derive a persistence factor a =1-g/2.  Hence an alpha, which they quote, of, say, 0.6 corresponds to a g of 0.8.

Extracting the ACF from temperature records is quite difficult and may lead to artefacts.

Any signal processing on data may alter the underlying statistics of the data.  It is important to use methods that do not distort the statistic in question or, at least, use a method that distorts it in a predictable manner.

The underlying temperature signal is dominated by the yearly variation as can be seen from the amplitude spectrum of the signal shown in figure 6, and this has to be removed to estimate the long-term persistence. The method by which this done can lead to quite different results for the ACF.  One approach might be to simply low-pass filter the signal with a cut-off point well above 1 year^-1 as shown by the purple filter response in figure 6.  This has a huge effect on the ACF as shown by the purple broken line in figure 5.  A selective band-pass filter can be used but this appears to produce “ringing” in the short term ACF and to alter its persistence in the tail.  The method, used here, which causes least disturbance to the ACF is to interpolate across the spectral bands in question and the effect of this filter on random noise is shown as the green line in figure 5.

Most signal processing methods contain compromises and the method has to carefully tailored to question that is being asked about the data.

Figure 5 The Mean ACF of the temperature records shown as blue squares with a 2 compartment model (black) and a power law model (red) fitted to the data.  The dashed purple line shows the effect of using a low pass filter to eliminate seasonal trends (i.e.: disastrous!) and the green line shows the effect of the seasonal de-trending method having very little effect on the ACF

Once the data has the yearly variation removed, it is normally distributed.

The Yearly ACF values are plotted as blue squares in figure 5 and the best-fit power law  (red g=0.875) and bi-exponential (black) models are shown.  In this data set, he bi-exponential model appears to fit slightly better, however, when the log ACF is plotted against the log delay, it is clear that there are considerable errors involved and it is not possible to say how closely the data conforms to a power law model.

Figure 6 Mean amplitude spectrum of the temperature records with peaks at 1,2,3 year^-1.  The red line shows a selective band pass filter and the purple line a low pass filter designed to remove seasonal trends.  Interpolation across these peaks gives least artefact in the ACF.

Figure 7 Mean ACF plotted on a log-log plot.  The vertical blue lines show the 95% limits of the distribution and the red dashed lines show a bi-compartment model.

It is not possible with this data set to be certain of the model to describe persistence, but one can use a power law or multi-compartment model with different assumed values  and determine their abilities to produce trends.

Figure 8 Simulated record (blue) with low pass filtered version (black) and trends superimposed in red.

Figure 8 shows the residual temperature (light blue) using a power-law model with g=0.875.  [This has been calculated using the ODE, not through the assumption of linearity and convolution]. The magnitude of disturbance has been drawn from a distribution so that the power in the model signal is drawn from the same distribution as in the observed seasonally detrended signals.  This is repeated 100000 times to give a distribution of the expected residual signals.   This has been repeated for gamma of 0.7.

The global temperature signal is considered to be the sum of a number of temperature signals drawn from stations that are perturbed with a highly correlated disturbance.  Since we have only one long-term temperature record, we cannot say anything about the statistics of this signal.  But, we can determine the ACF of the simulated signals and these show considerable variability as shown in figure 9.  The mean ACF and its 95% limits are shown for a power law model of g=0.7 (red) and for a 2 compartment model which approximates g=0.9.  Irrespective of how persistence is extracted from these signals, there is a considerable range of g fitted to the signal and the power law and 2 compartment model are indistinguishable given a single temperature record. The 95% limits for a true g =0.7 is 0.58-0.87.

Therefore, the model suggests that if we determine the persistence model from one set of global temperature records, there is likely to be considerable errors in the estimate.

Figure 9 Mean ACFs and their 95% error limits for a power law model (Red-g=0.7) and a 2 compartment model (Black)with an approximate g of 0.9.

Estimation of trends in the simulated signal.

Estimation of trends in a noisy signal is a difficult process and can be done in a number of ways. The results depend to some extent on how one defines a trend.

Figure 10

The black line in Figure 8 shows a low-pass filtered version of the signal with trend lines plotted in red.  There are clearly shorter-term trends within the signal and, if one decomposes a signal into very short-term trends, it is not possible to distinguish a long-term trend.  The approach used here, is to low pass filter the signal heavily, so as to eliminate the rapid fluctuations and use the maxima and minima of this signal to determine the ends of trends.  The trend line is then fitted to the data between these intervals.

This is shown more clearly in figure 10 where the low pass filtered signal is shown, and the trends extracted, using a filter with cut off points at 1/30 years (solid red line) and at /100 years (dotted red line).  In general, the magnitude of the trend is less as longer trends are identified and there will be more variability around the trend.

Figure 11 Histograms of trend lengths and excursions.  Power law processes with , A,  g=0.9  (upper) and, B, g=0.7 (lower). The trends > 100 years are shaded red and the region between 60 and 100 years is also shaded

Using this approach, the distribution of trends for power law gammas of 0.9 and 0.7 (i.e.: a=0.55 and 0.65) can be calculated and are shown in figure 11.  According to these calculations the probability of a 100 year trend in either direction with g=0.9 is 0.02 with a mean value of 0.27^o C.  When g =0.7, the probability of a 100 year trend is 0.09 with a mean value of 0.4^o , while the probability of a >100 year trend of >0.7^oC is  0.02.

Conclusions

The purpose of the analysis presented by Ludecke, Link and Ewert is to characterise the Earth’s temperature as a system, in particular the lag in temperature with energy changes.

Any system that involves lag can produce trends.  Whether the lag is imposed via a Hurst process, a multi-compartment system or anything else is a secondary consideration.  However, the problem is to characterise the process sufficiently accurately to make predictions about the magnitudes of trends is another matter.  The data available is strictly limited and one cannot easily distinguish between different physical models that may underlie the trends.  Satellite temperature and radiation flux data will eventually lead to much better temperature series but this data would be better incorporated into physical models that incorporate the mechanisms of energy storage and would render the approach taken here unnecessary.

One of the major problems in this analysis is that temperature signals at different points on the Earth are not independent.  If heat is stored in the ocean and transferred by currents to another point, or by some other mechanism, the temperature rise at this point is not an “external trend” as stated, which has to be eliminated, but is an important part of the mechanism of the creation of regional trends within the system as a whole.

Figure 12 A model of energy transfer between different regions with peristence

The addition of energy into a region, as shown in figure 12, will give a false value for the persistence.  However as system similar to that in figure 12 is analysable, in principle, provided the system is linear[1]. If one is using a power law representation of the persistence, the problem is far more difficult.  Essentially, one writes the energy flux equation for each component in the system as a set of coupled differential equations, or in the complex frequency domain, and solves then to obtain the parameters of the persistence in each element.  This can be done via the ACF of the temperature at each point.  In practice, this is practically impossible with the data available and the solution is highly unstable and very sensitive to small errors in data, even for a restricted subset of the problem.  However, LL&E have attempted to address this problem by de-trending the data, which is an arbitrary approach and smoothes the data, which may not necessarily give a better value for the calculated persistence.

A major problem is the quality of the data on which a model is constructed.  As I have pointed out, using correlation methods to determine persistence, long-term records are required and the number of “clean” long term records, especially in the tropics, are strictly limited and, to mind, global conclusions based on this data should treated sceptically.

However, applying a pure persistence model, it is conceivable that the increase in temperature of 0.7oC over 100 years is due to persistence but unlikely (<0.01).  On the other hand, trends of 40 years of 0-3 –0.4oC overlaid on a steadily increasing temperature are likely and may contribute to some of the relatively short-term fluctuations seen in the temperature record.