ATMO/GEOS 595c
Patterns and Mechanisms of Decadal Climate Variability
Data Exercise 2: Testing time series periodicities
vs. a red noise null hypothesis
Question for discussion Monday, Feb
18th, as
a result of this data
analysis exercise:
- Are there decadal periodicities in the
Mann et al. (2000) regional reconstructions, which are significant with
respect to a red noise hypothesis?
Introduction
In Data Exercise 1, many of you
suggested that there was decadal variability in the lowpass filtered
dataseries illustrated above. In this exercise we'll develop a
Monte Carlo
estimate of the power spectra of these series we expect to see
resulting from modeling the data series as simple autoregressive
functions, as in Section 2 of
Pierce
(2001). construct a
simple decadal-timescale filter, and apply it to these data series (
here and
here,
extracted and de-headered from
here
and
here,
respectively). Given your results and reading of the papers
assigned for this week and last week, bring your answers to these
questions for discussion, Monday, 4 February 2008.
Feel free to employ any computational environment you prefer. I
will give some direction on how to do this in
matlab, which is currently
sitelicensed at the
University of Arizona (freeware
Octave
workalike is here). A short introduction to operation, nuts and
bolts coding, and displaying graphics in matlab from
another course I offer
is
here;
a rough set of computing principles is
here.
See also the links in this assignment to further explanation of
terminology and concepts.
A good general reference for time series analysis is Chatfield,
Analysis
of Time Series: An Introduction (6th ed., 2003).
Part
I. Construction
of an autoregressive model of order 1.
The order 1 autoregressive model is an appealing choice for a null
hypothesis in oceanography and climatology, for at least two
reasons. First, as we saw in Pierce (2001) and more
famously in Hasselman (1977) and Hasselman and Frankingnoul (1977),
there are good reasons why low frequency climate variations might
resemble the integration of white noise atmospheric forcing into a more
slowly varying ocean. Second, we can easily write down analytical
functions for the spectral dependence of an AR(1) process on frequency
(see here
for an introduction; also Wilks,
2006).
Let's load our data series. In matlab this is
% load the
Mann N. American-average warm season reconstructions
% data
source: http://www.ncdc.noaa.gov/paleo/ei/ei_data/namerwarm-recon.dat
% publication
url: http://www.ncdc.noaa.gov/paleo/ei/ei_cover.html
load
namerwarm_rec.dat
nyr=namerwarm_rec(:,1);
nt=namerwarm_rec(:,2);
% load the
Mann N. American-average warm season reconstructions
% data
source: http://www.ncdc.noaa.gov/paleo/ei/ei_data/namerwarm-recon.dat
% publication
url: http://www.ncdc.noaa.gov/paleo/ei/ei_cover.html
load
tropicc_rec.dat
tyr=tropicc_rec(:,1);
tt=tropicc_rec(:,2);
Let's focus on just the 1750-1900 time interval, as the 20th century is
dominated by a trend in both series:
% power
spectra: for pre-instrumental period
nt=nt(1:151);
%1750-1900 period
tt=tt(1:151);
%1750-1900 period
d1=(nt-mean(nt))./std(nt);
d2=(tt-mean(tt))./std(tt);
Calculate the power spectra of the two dataseries and store it away for
plotting later, e.g.
[pd1,fd1]=pmtm(d1,[],[],1);
The form of a first-order autoregressive function X is:
Xt = rXt-1 + N(0,sx)
with Xt the tth data point, r the lag-1
autocorrelation of X, and N the normal distribution of random numbers
with mean zero and standard deviation sx, of the series X.
So first we need the autocorrelation of X at lag 1: this is just the
correlation of X with a copy of itself shifted in time by one time
interval. For instance, consider the row vector Xt
with nt data points. Let d1a be the first nt-1 points of X, with
mean removed, and standardized by the variance of the first N-1 points
of X. Let d1b be constructed the same way as d1a, except now it
is the centered and standardized form of the second through nt points
of X. Then the AR(1) coefficient r for this series is
r = d1a*d1b'/(nt-2)
[with the prime denoting the vector transpose]
Find the AR(1) coefficient for both the North American warm season
regional average series, and the Tropical cold season regional average
series, e.g.
% find AR(1)
coefficients
l1=length(d1);
d1a=(d1(1:l1-1)-mean(d1(1:l1-1)))'./std(d1(1:l1-1));
d1b=(d1(2:l1)-mean(d1(2:l1)))'./std(d1(2:l1))';
ar1d1=d1a*d1b'./(l1-2);
Now we can construct AR(1) series having AR(1) coefficients and lengths
of our actual data series. Construct 100 AR(1) series
corresponding to each AR(1) coefficient. For instance, for series
d1, you might code
% generate AR(1)
series
nr=100;
rd1=randn(nr,l1);
s1=std(d1);
for i=2:l1
rd1(:,i)=ar1d1.*rd1(:,i-1)+s1.*randn(nr,1);
end
for i=1:nr % renormalize to zero mean and unit standard deviation
rd1(i,:)=(rd1(i,:)-mean(rd1(i,:)))./std(rd1(i,:));
rd2(i,:)=(rd2(i,:)-mean(rd2(i,:)))./std(rd2(i,:));
end
Last, calculate power spectra for each of the 100 AR(1) data series we
constructed, for each of the two cases, in the same manner as you did
for the actual data series.
prd1=zeros(nr,length(fd1));
for i=1:nr
[prd1(i,:),fr1]=pmtm(detrend(rd1(i,:)),[],[],1);
end
For each case, you now have 100 AR(1) power spectra, calculated at nf
frequencies. Sort these red noise power for each frequency by
power; in other words, if your spectral results for the AR(1) series
are 100 rows with nf columns of frequencies, then sort each column of
spectral estimates from highest to lowest, e.g.
prd1s=flipud(sort(prd1));
You can now estimate the 95th and 5th percent confidence levels for
spectral power at each frequency by plotting the 5th and 95th highest
spectral estimates for each frequency for which you calculated a
spectral estimate.
Part II. Interpretation.
For each data series, plot its power spectrum. Overlay on this
power spectrum the 5th and 95th percent confidence levels you computed
based on the corresponding AR(1) null hypothesis. Can you reject
the null hypothesis that there is no significant spectral power above
that of red noise?
Now go back to the top
of this page
and prepare your thoughts on the answer to our discussion question.
back to Syllabus/Schedule. Last updated 14 Feb 2008.
