Digital Signal Processing analysis of global temperature data time series suggests global cooling ahead

This DSP engineer is often tasked with extracting spurious signals from noisy data. He submits this interesting result of applying these techniques to the HadCRUT temperature anomaly data. Digital Signal Processing analysis suggests cooling ahead in the immediate future with no significant probability of a positive anomaly exceeding .5°C between 2023 and 2113. See figures 13 and 14. Code and data is made available for replication. – Anthony

Guest essay by Jeffery S. Patterson, DSP Design Architect, Agilent Technologies

Harmonic Decomposition of the Modern Temperature Anomaly Record

Abstract: The observed temperature anomaly since 1900 can be well modeled with a simple harmonic decomposition of the temperature record based on a fundamental period of 170.7 years. The goodness-of-fit of the resulting model significantly exceeds the expected fit to a stochastic AR sequence matching the general characteristic of the modern temperature record.

Data

I’ve used the monthly Hadcrut3 temperature anomaly data available from http://woodfortrees.org/data/hadcrut3vgl/every as plotted in Figure 1.

Figure 1 – Hadcrut3 Temperature Record 1850-Present

To remove seasonal variations while avoiding spectral smearing and aliasing effects, the data was box-car averaged over a 12-month period and decimated by 12 to obtain the average annual temperature plotted in Figure 2.

Figure 2 – Monthly data decimated to yearly average

A Power Spectral Density (PSD) plot of the decimated data reveals harmonically related spectral peaks.

Figure 3 – PSD of annual temperature anomaly in dB

To eliminate the possibility that these are FFT (Fast Fourier Transform) artifacts while avoiding the spectral leakage associated with data windowing, we use a technique is called record periodization. The data is regressed about a line connecting the record endpoints, dropping the last point in the resulting residual. This process eliminates the endpoint discontinuity while preserving the position of the spectral peaks (although it does extenuate the amplitudes at higher frequencies and modifies the phase of the spectral components). The PSD of the residual is plotted in Figure 4.

Figure 4 – PSD of the periodized record

Since the spectral peaking is still present we conclude these are not record-length artifacts. The peaks are harmonically related, with odd harmonics dominating until the eighth. Since spectral resolution increases with frequency, we use the eighth harmonic of the periodized PSD to estimate the fundamental. The following Mathematica (Mma) code finds the 5^th peak (8^th harmonic) and estimates the fundamental.

wpkY1=Abs[ArgMax[{psdY,w>.25},w]]/8

0.036811

The units are radian frequency across the Nyquist band, mapped to ±p (the plots are zoomed to 0 < w < 1 to show the area of interest). To convert to years, invert wpkY1 and multiply by 2p, which yields a fundamental period of 170.7 years.

From inspection of the PSD we form the harmonic model (note all of the radian frequencies are harmonically related to the fundamental):

(*Define the 5^th order harmonic model used in curve fit*)

model=AY1*Sin[wpkY1 t+phiY1]+AY2*Sin[2*wpkY1* t+phiY2]+AY3*

Sin[3*wpkY1* t+phiY3]+AY4*Sin[4*wpkY1* t+phiY4]+AY5*

Sin[5*wpkY1* t+phiY5]];

vars= {AY1,phiY1,AY2,phiY2,AY3,phiY3,AY4,phiY4,AY5,phiY5 }

and fit the model to the original (unperiodized) data to find the unknown amplitudes, AYx, and phases, phiYx.

fitParms1=FindFit[yearly,model,vars,t]

fit1=Table[model/.fitParms1,{t,0,112}];

residualY1= yearly- fit1;{AY1→-0.328464,phiY1→1.44861,AY2→-

0.194251,phiY2→3.03246,AY3→0.132514,phiY3→2.26587,AY4→0.0624932,

phiY4→-3.42662,AY5→-0.0116186,phiY5→-1.36245,AY8→0.0563983,phiY8→

1.97142,wpkY1→0.036811}

The fit is shown in Figure 5 and the residual error in Figure 6.

Figure 5 – Harmonic model fit to annual data

Figure 6 – Residual Error Figure 7 – PSD of the residual error

The residual is nearly white, as evidenced by Figure 7, justifying use of the Hodric-Prescott filter on the decimated data. This filter is designed to separate cyclical, non-stationary components from data. Figure 8 shows an excellent fit with a smoothing factor of 15.

Figure 8 – Model vs. HP Filtered data (smoothing factor=3)

Stochastic Analysis

The objection that this is simple curve fitting can be rightly raised. After all, harmonic decomposition is a highly constrained form of Fourier analysis, which is itself a curve fitting exercise that yields the harmonic coefficients (where the fundamental is the sample rate) which recreate the sequence exactly in the sample domain. That does not mean however, that any periodicity found by Fourier analysis (or by implication, harmonic decomposition) are not present in the record. Nor, as will be shown below, is it true that harmonic decomposition on an arbitrary sequence would be expected to yield the goodness-of-fit achieved here.

The 113 sample record examined above is not long enough to attribute statistical significance to the fundamental 170.7 year period, although others have found significance in the 57-year (here 56.9 year) third harmonic. We can however, estimate the probability that the results are a statistical fluke.

To do so, we use the data record to estimate an AR process.

procY=ARProcess[{a1,a2,a3,a4,a5},v];

procParamsY = FindProcessParameters[yearlyTD["States"],procY]

estProcY= procY /. procParamsY

WeakStationarity[estProcY]

{a1→0.713,a2→0.0647,a3→0.0629,a4→0.181,a5→0.0845,v→0.0124391}

As can be seen in Figure 9 below, the process estimate yields a reasonable match to observed power spectral density and covariance function.

Figure 9 – PSD of estimated AR process (red) vs. data Figure 9b – Correlation function (model in blue)

Figure 10 – 500 trial spaghetti plot Figure 10b – Three paths chosen at random

As shown in 10b, the AR process produces sequences which in general character match the temperature record. Next we perform a fifth-order harmonic decomposition on all 500 paths, taking the variance of the residual as a goodness-of-fit metric. Of the 500 trials, harmonic decomposition failed to converge 74 times, meaning that no periodicity could be found which reduced the variance of the residual (this alone disproves the hypothesis that any arbitrary AR sequences can be decomposed). To these failed trials we assigned the variance of the original sequence. The scattergram of results are plotted in Figure 11 along with a dashed line representing the variance of the model residual found above.

Figure 11 – Variance of residual; fifth order HC (Harmonic Coefficients), residual 5HC on climate record shown in red

We see that the fifth-order fit to the actual climate record produces an unusually good result. Of the 500 trials, 99.4% resulted in residual variance exceeding that achieved on the actual temperature data. Only 1.8% of the trials came within 10% and 5.2% within 20%. We can estimate the probability of achieving this result by chance by examining the cumulative distribution of the results plotted in Figure 12.

Figure 12 – CDF (Cumulative Distribution Function) of trial variances

The CDF estimates the probability of achieving these results by chance at ~8.1%.

Forecast

Even if we accept the premise of statistical significance, without knowledge of the underlying mechanism producing the periodicity, forecasting becomes a suspect endeavor. If for example, the harmonics are being generated by a stable non-linear climatic response to some celestial cycle, we would expect the model to have skill in forecasting future climate trends. On the other hand, if the periodicities are internally generated by the climate itself (e.g. feedback involving transport delays), we would expect both the fundamental frequency and importantly, the phase of the harmonics to evolve with time making accurate forecasts impossible.

Nevertheless, having come thus far, who could resist a peek into the future?

We assume the periodicity is externally forced and the climate response remains constant. We are interested in modeling the remaining variance so we fit a stochastic model to the residual. Empirically, we found that again, a 5^th order AR (autoregressive) process matches the residual well.

tDataY=TemporalData[residualY1-Mean[residualY1],Automatic];

yearTD=TemporalData[residualY1,{ DateRange[{1900},{2012},"Year"]}]

procY=ARProcess[{a1,a2,a3,a4,a5},v];

procParamsY = FindProcessParameters[yearTD["States"],procY]

estProcY= procY /. procParamsY

WeakStationarity[estProcY]

A 100-path, 100-year run combining the paths of the AR model with the harmonic model derived above is shown in Figure 13.

Figure 13 – Projected global mean temperature anomaly (centered 1950-1965 mean)

Figure 14 – Survivability at 10 (Purple), 25 (Orange), 50 (Red), 75 (Blue) and 100 (Green) years

The survivability plots predict no significant probability of a positive anomaly exceeding .5°C between 2023 and 2113.

Discussion

With a roughly one-in-twelve chance that the model obtained above is the manifestation of a statistical fluke, these results are not definitive. They do however show that a reasonable hypothesis for the observed record can be established independent of any significant contribution from greenhouse gases or other anthropogenic effects.