Milankovic Theory andTime Series Analysis

Mudelsee M Institute of Meteorology University

of Leipzig Germany

Climate Statistical analysis

Data (sample) Climate

system (population, truth, theory)

Climate Statistical analysis

Data (sample) STATISTICS Climate

system (population, truth, theory)

Climate Statistical analysis Time series

analysis

Sample t(i), x(i), i 1, ..., n

Climate Statistical analysis Time series

analysis

Sample t(i), x(i), i 1, ...,

n UNI-VARIATE TIME SERIES

Climate Statistical analysis Time series

analysis

Sample t(i), x(i), y(i), i 1, ...,

n BI-VARIATE TIME SERIES

Climate Statistical analysis Time series

analysis

Sample t(i), x(i), y(i), i 1, ...,

n TIME SERIES DYNAMICS

Climate Statistical analysis Time series

analysis

Sample t(i), x(i), y(i), i 1, ...,

n TIME SERIES DYNAMICS t(i), x(i),

y(i), z(i),..., i 1 TIME SLICE STATICS

Climate Statistical analysis Time series

analysis

Sample t(i), x(i), y(i), i 1, ...,

n CLIMATE TIME SERIES

Climate Statistical analysis Time series

analysis

Sample t(i), x(i), y(i), i 1, ...,

n CLIMATE TIME SERIES o uneven time

spacing

UNEVEN TIME SPACING

Mudelsee M (in prep.) Statistical Analysis of

Climate Time Series A Bootstrap Approach. Kluwer.

ICE CORE (Vostok dD)

TREE RINGS (atmospheric ?14C)

STALAGMITE (Qunf Cave d18O)

UNEVEN TIME SPACING

Mudelsee M (in prep.) Statistical Analysis of

Climate Time Series A Bootstrap Approach. Kluwer.

Climate Statistical analysis Time series

analysis

Sample t(i), x(i), y(i), i 1, ...,

n CLIMATE TIME SERIES o uneven time

spacing o persistence

PERSISTENCE

Mudelsee M (in prep.) Statistical Analysis of

Climate Time Series A Bootstrap Approach. Kluwer.

ICE CORE (Vostok dD)

TREE RINGS (atmospheric ?14C)

STALAGMITE (Qunf cave d18O)

PERSISTENCE

Mudelsee M (in prep.) Statistical Analysis of

Climate Time Series A Bootstrap Approach. Kluwer.

Climate Statistical analysis Time series

analysis

Sample t(i), x(i), y(i), i 1, ...,

n CLIMATE TIME SERIES o uneven time

spacing o persistence

Milankovic theory

Theory Orbital variations influence Earths

climate.

Milankovic theory

Data Climate time series Theory Orbital

variations influence Earths climate.

Milankovic theory

Data Climate time series TIME SERIES

ANALYSIS TEST Theory Orbital variations

influence Earths climate.

Milankovic theory andtime series analysis

Part 1 Spectral analysis Part 2 Milankovic

paleoclimate back to the Pliocene

Acknowledgements

Berger A, Berger WH, Grootes P, Haug G, Mangini

A, Raymo ME, Sarnthein M, Schulz M, Stattegger K,

Tetzlaff G, Tong H, Yao Q, Wunsch C

Alert!

Mudelsee-bias

Part 1 Spectral analysis

Sample t(i), x(i), y(i), i 1, ...,

n Simplification uni-variate, only

x(i), equidistance, t(i) i

Part 1 Spectral analysis

Sample x(t) Time series

Part 1 Spectral analysis

Sample x(t) Time series Population X(t)

Part 1 Spectral analysis

Sample x(t) Time series Population X(t) Pr

ocess

Part 1 Spectral analysis Process level

TIME DOMAIN

X(t)

Part 1 Spectral analysis Process level

TIME DOMAIN

X(t)

FOURIER TRANSFORMATION FREQUENCY DOMAIN

Part 1 Spectral analysis Process level

TIME DOMAIN

X(t) T GT(f) (2p)1/2?T XT(t) e2pift

dt, XT X(t), T t T, 0, otherwise.

FOURIER TRANSFORMATION FREQUENCY DOMAIN

Part 1 Spectral analysis Process level

h(f) limT?8 E GT(f)2 / (2T)

NON-NORMALIZED POWER SPECTRAL DENSITY

FUNCTION, SPECTRUM

Part 1 Spectral analysis Process level

h(f) limT?8 E GT(f)2 / (2T)

NON-NORMALIZED POWER SPECTRAL DENSITY

FUNCTION, SPECTRUM ENERGY (VARIATION) AT

SOME FREQUENCY

Part 1 Spectral analysis Process level

Discrete spectrum Harmonic process Astronomy

Part 1 Spectral analysis Process level

Discrete spectrum Harmonic process Astronomy

Continuous spectrum Random process Climatic noise

Part 1 Spectral analysis Process level

Discrete spectrum Harmonic process Astronomy

Continuous spectrum Random process Climatic noise

Mixed spectrum Typical climatic

Part 1 Spectral analysis

The task of spectral analysis is to estimate the

spectrum. There exist many estimation

techniques.

Part 1 Spectral analysis Harmonic regression

HARMONIC PROCESS

X(t) Sk Ak cos(2pfk t) Bk sin(2pfk t)

e(t)

Part 1 Spectral analysis Harmonic regression

HARMONIC PROCESS

X(t) Sk Ak cos(2pfk t) Bk sin(2pfk t)

e(t) If frequencies fk known a

priori Minimize Q Si x(i) Sk Ak

cos(2pfk t) Bk sin(2pfk t)2 to obtain Ak

and Bk.

Part 1 Spectral analysis Harmonic regression

HARMONIC PROCESS

X(t) Sk Ak cos(2pfk t) Bk sin(2pfk t)

e(t) If frequencies fk known a

priori Minimize Q Si x(i) Sk Ak

cos(2pfk t) Bk sin(2pfk t)2 to obtain Ak

and Bk.

LEAST SQUARES

Part 1 Spectral analysis Periodogram

If frequencies fk not known a priori Take

least-squares solutions Ak and Bk, fk 0,

1/n, 2/n, ..., 1/2, to calculate P(fk) (Ak)2

(Bk)2.

Part 1 Spectral analysis Periodogram

If frequencies fk not known a priori Take

least-squares solutions Ak and Bk, fk 0,

1/n, 2/n, ..., 1/2, to calculate P(fk) (Ak)2

(Bk)2.

PERIODOGRAM

Part 1 Spectral analysis Periodogram

If frequencies fk not known a priori Take

least-squares solutions Ak and Bk, fk 0,

1/n, 2/n, ..., 1/2, to calculate P(fk) (Ak)2

(Bk)2. Where fk true f P(fk) has a peak.

PERIODOGRAM

Part 1 Spectral analysis Periodogram

Part 1 Spectral analysis Periodogram

Original paper Schuster A (1898) On the

investigation of hidden periodicities with

application to a supposed 26 day period

of meteorological phenomena. Terrestrial

Magnetism 31341.

Part 1 Spectral analysis Periodogram

Hypothesis test (significance of periodogram

peaks) Fisher RA (1929) Tests of significance

in harmonic analysis. Proceedings of the Royal

Society of London, Series A, 1255459.

Part 1 Spectral analysis Periodogram

A wonderful textbook Priestley MB

(1981) Spectral Analysis and Time

Series. Academic Press, London, 890 pp.

Part 1 Spectral analysis Periodogram

Major problem with the periodogram as spectrum

estimate Relative error of P(fk) 200 for

fk 0, 1/2, 100 otherwise.

Part 1 Spectral analysis Periodogram

Part 1 Spectral analysis Periodogram

More lives have been lost looking at the raw

periodogram than by any other action involving

time series! Tukey JW (1980) Can we predict

where time series should go next? In

Directions in time series analysis (eds

Brillinger DR, Tiao GC). Institute of

Mathematical Statistics, Hayward, CA, 131.

Part 1 Spectral analysis Smoothing

(No Transcript)

WELCH OVERLAPPED SEGMENT AVERAGING (WOSA)

1st Segment

2nd Segment

3rd Segment

WELCH OVERLAPPED SEGMENT AVERAGING (WOSA)

1st Segment

2nd Segment

3rd Segment

ERROR REDUCTION lt v3

Part 1 Spectral analysis Smoothing

Tapering Weight time series Spectral

leakage reduced (Hanning, Parzen, triangula

r windows, etc.)

Part 1 Spectral analysis Smoothing problem

Several segments averaged Spectrum estimate more

accurate -) Fewer (n lt n) data per

segment Lower frequency resolution -(

Part 1 Spectral analysis Smoothing problem

Part 1 Spectral analysis Smoothing problem

Subjective judgement is unavoidable. Play with

parameters and be honest.

Part 1 Spectral analysis 100-kyr problem

?t 1 fk 0, 1/n, 2/n, ... ?t d fk 0,

1/(nd), 2/(n d), ... ?f (nd)1

Part 1 Spectral analysis 100-kyr problem

?t 1 fk 0, 1/n, 2/n, ... ?t d fk 0,

1/(nd), 2/(n d), ... ?f (nd)1 BW gt

(nd)1 SMOOTHING

Part 1 Spectral analysis 100-kyr problem

nd 650 kyr ?f (650 kyr)1

Part 1 Spectral analysis 100-kyr problem

nd 650 kyr ?f (650 kyr)1 (100 kyr)1

?f (118 kyr)1 to (87 kyr)1

Part 1 Spectral analysis 100-kyr problem

nd 650 kyr ?f (650 kyr)1 (100 kyr)1

?f (118 kyr)1 to (87 kyr)1

BW wider SMOOTHING

Part 1 Spectral analysis 100-kyr problem

The 100-kyr cycle existed not long enough to

allow a precise enough frequency estimation.

Part 1 Spectral analysis BlackmanTukey

h Fourier transform of ACV

Part 1 Spectral analysis BlackmanTukey

E X(t) X(t lag) h Fourier

transform of ACV

Part 1 Spectral analysis BlackmanTukey

PROCESS LEVEL E X(t) X(t lag) h

Fourier transform of ACV

Part 1 Spectral analysis BlackmanTukey

PROCESS LEVEL E X(t) X(t lag) h

Fourier transform of ACV SAMPLE S x(t)

x(t lag) / n h Fourier transform of ACV

Part 1 Spectral analysis BlackmanTukey

Fast Fourier Transform Cooley JW, Tukey JW

(1965) An algorithm for the machine

calculation of complex Fourier

series. Mathematics of Computation 19297301.

Part 1 Spectral analysis BlackmanTukey

Some paleoclimate papers Hays JD, Imbrie J,

Shackleton NJ (1976) Variations in the Earth's

orbit Pacemaker of the ice ages. Science

19411211132. Imbrie J Hays JD, Martinson DG,

McIntyre A, Mix AC, Morley JJ, Pisias NG, Prell

WL, Shackleton NJ (1984) The orbital theory of

Pleistocene climate Support from a revised

chronology of the marine d18O record. In

Milankovitch and Climate (eds Berger A, Imbrie

J, Hays J, Kukla G, Saltzman B), Reidel,

Dordrecht, 269305.

Part 1 Spectral analysis BlackmanTukey

Ruddiman WF, Raymo M, McIntyre A (1986)

Matuyama 41,000-year cycles North Atlantic

Ocean and northern hemisphere ice sheets. Earth

and Planetary Science Letters 80117129. Tiedema

nn R, Sarnthein M, Shackleton NJ (1994)

Astronomic timescale for the Pliocene Atlantic

d18O and dust flux records of Ocean Drilling

Program Site 659. Paleoceanography 9619638.

Part 1 Spectral analysis Multitaper Method

(MTM)

Spectral estimation with optimal

tapering Thomson DJ (1982) Spectrum estimation

and harmonic analysis. Proceedings of the IEEE

7010551096. MINIMAL DEPENDENCE AMONG AVERAGED

INDIVIDUAL SPECTRA MINIMAL ESTIMATION ERROR

Part 1 Spectral analysis Multitaper Method

(MTM)

Part 1 Spectral analysis Multitaper Method

(MTM)

BETTER DIRECTLY VIA ASTRONOMY EQS.

Part 1 Spectral analysis Multitaper Method

(MTM)

Some paleoclimate papers Park J, Herbert TD

(1987) Hunting for paleoclimatic periodicities in

a geologic time series with an uncertain time

scale. Journal of Geophysical Research

921402714040. Thomson DJ (1990)

Quadratic-inverse spectrum estimates

Applications to palaeoclimatology. Philosophical

Transactions of the Royal Society of London,

Series A 332539597. Berger A, Melice JL,

Hinnov L (1991) A strategy for frequency spectra

of Quaternary climate records. Climate Dynamics

5227240.

Part 1 Spectral analysis Further points

Uneven time spacing

Part 1 Spectral analysis Further points

Uneven time spacing Use X(t) Sk Ak cos(2pfk

t) Bk sin(2pfk t) e(t) Lomb NR (1976)

Least-squares frequency analysis of

unequally spaced data. Astrophysics and Space

Science 39447462. Scargle JD (1982) Studies in

astronomical time series analysis.

II. Statistical aspects of spectral analysis of

unevenly spaced data. The Astrophysical Journal

263835853.

HARMONIC PROCESS

Part 1 Spectral analysis Further points

PERSISTENCE

Red noise

Part 1 Spectral analysis Further points

PERSISTENCE

Red noise AR1 process for uneven

spacing Robinson PM (1977) Estimation of a time

series model from unequally spaced data.

Stochastic Processes and their Applications

6924.

Part 1 Spectral analysis Further points

Aliasing

Part 1 Spectral analysis Further points

Aliasing Safeguards o uneven spacing

(Priestley 1981) o for marine records

bioturbation Pestiaux P, Berger A (1984) In

Milankovitch and Climate, 493510.

Part 1 Spectral analysis Further points

Running window Fourier Transform

Priestley MB (1996) Wavelets and time-dependent

spectral analysis. Journal of Time Series

Analysis 1785103.

Part 1 Spectral analysis Further points

Detrending

Part 1 Spectral analysis Further points

Errors in t(i) tuned dating, absolute

dating, stratigraphy. Errors in

x(i) measurement error, proxy

error, interpolation error.

Part 1 Spectral analysis Further points

Bi-variate spectral analysis For example x

marine d18O y insolation

Part 1 Spectral analysis Further points

Higher-order spectra (bi-spectra, ...)

Part 1 Spectral analysis Further points

Etc., etc.

Part 2 Milankovic paleoclimate

Part 2 Milankovic paleoclimate

Part 2 Milankovic paleoclimate

Northern Hemisphere Glaciation NHG

Part 2 Milankovic paleoclimate

Mid-Pleistocene Transition

Northern Hemisphere Glaciation NHG

Part 2 Milankovic paleoclimate Climate

transitions, trend

Part 2 Milankovic paleoclimate Climate

transitions, trend

x1, t lt t1, Xfit(t) x2, t gt

t2, x1 (t-t1) (x2-x1)/(t2-t1), t1 t

t2.

Part 2 Milankovic paleoclimate Climate

transitions, trend

x1, t lt t1, Xfit(t) x2, t gt

t2, x1 (t-t1) (x2-x1)/(t2-t1), t1 t

t2.

LEAST SQUARES ESTIMATION

Part 2 Milankovic paleoclimate Mid-Pleistocen

e Transition

Part 2 Milankovic paleoclimate Mid-Pleistocen

e Transition

100 kyr cycle

Part 2 Milankovic paleoclimate Mid-Pleistocen

e Transition

size of Barents/ Kara Sea ice sheets

DSDP 552 DSDP 607 ODP 659 ODP 677 ODP 806

Mudelsee M, Schulz M (1997) Earth and Planetary

Science Letters 151117123.

Part 2 Milankovic paleoclimate NHG

Database 24 Myr, 45 marine d18O records, 4

temperature records

benthic planktonic

Mudelsee M, Raymo ME (submitted)

NHGResults

High-resolution records

Mudelsee M, Raymo ME (submitted)

NHGResults

High-resolution records

Mudelsee M, Raymo ME (submitted)

Part 2 Milankovic paleoclimate NHG

NHG was a slow global climate change (from 3.6

to 2.4 Myr). NHG ice volume signal 0.4 .

Milankovic theory andtime series analysis

Conclusions

(1) Spectral analysis estimates

the spectrum. (2) Trend estimation is

also important (climate transitions).

G O O D I E S

Climate transitions error bars

Time series, size n

t(i), x(i) i 1,, n

t(i)

Simulated time series, x(i) ramp noise

t(i), x(i)

Simulated ramp parameters

t1, x1, t2, x2

Repeat 400 times

t1, x1, t2, x2

Ramp estimation

Take standard deviation of simulated ramp

parameters

Bootstrap errors

STD, Persistence

Noise estimation

NHG amplitudes temperature

cooling (C) in 3,606-2,384 kyr 0.12

0.47 0.62 0.29 1.0 0.5 -0.85 0.17

Mudelsee M, Raymo ME (submitted)

NHG amplitudes ice volume

Temperature calibration Dd18OT/DT -0.234

0.003 /C (Chen 1994 own error

determination) Salinity calibration Dd18OS

/DT 0.05 /C (Whitman and Berger

1992) DSDP 572 p Dd18OT 0.03 0.12

Dd18OS -0.01 Dd18OI 0.34 0.13 DSDP

607 b Dd18OT 0.15 0.07 Dd18OS -0.03

Dd18OI 0.41 0.09 ODP 806 b Dd18OT 0.24

0.12 Dd18OS - 0.05 Dd18OI 0.25 0.13

ODP 806 p Dd18OT -0.20 0.04 Dd18OS

0.04 Dd18OI 0.43 0.06 (DSDP 1085

b cooling by 1 C Dd18OI 0.35

) Average Dd18OI 0.39 0.04