Milankovic Theory and Time Series Analysis

1
Milankovic Theory andTime Series Analysis
Mudelsee M Institute of Meteorology University
of Leipzig Germany
2
Climate Statistical analysis
Data (sample) Climate
system (population, truth, theory)
3
Climate Statistical analysis
Data (sample) STATISTICS Climate
system (population, truth, theory)
4
Climate Statistical analysis Time series
analysis
Sample t(i), x(i), i 1, ..., n
5
Climate Statistical analysis Time series
analysis
Sample t(i), x(i), i 1, ...,
n UNI-VARIATE TIME SERIES
6
Climate Statistical analysis Time series
analysis
Sample t(i), x(i), y(i), i 1, ...,
n BI-VARIATE TIME SERIES
7
Climate Statistical analysis Time series
analysis
Sample t(i), x(i), y(i), i 1, ...,
n TIME SERIES DYNAMICS
8
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
9
Climate Statistical analysis Time series
analysis
Sample t(i), x(i), y(i), i 1, ...,
n CLIMATE TIME SERIES
10
Climate Statistical analysis Time series
analysis
Sample t(i), x(i), y(i), i 1, ...,
n CLIMATE TIME SERIES o uneven time
spacing

11
UNEVEN TIME SPACING
Mudelsee M (in prep.) Statistical Analysis of
Climate Time Series A Bootstrap Approach. Kluwer.
12
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.
13
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

14
PERSISTENCE
Mudelsee M (in prep.) Statistical Analysis of
Climate Time Series A Bootstrap Approach. Kluwer.
15
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.
16
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
17
Milankovic theory
Theory Orbital variations influence Earths
climate.
18
Milankovic theory
Data Climate time series Theory Orbital
variations influence Earths climate.
19
Milankovic theory
Data Climate time series TIME SERIES
ANALYSIS TEST Theory Orbital variations
influence Earths climate.
20
Milankovic theory andtime series analysis
Part 1 Spectral analysis Part 2 Milankovic
paleoclimate back to the Pliocene
21
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
22
Alert!
Mudelsee-bias
23
Part 1 Spectral analysis
Sample t(i), x(i), y(i), i 1, ...,
n Simplification uni-variate, only
x(i), equidistance, t(i) i
24
Part 1 Spectral analysis
Sample x(t) Time series
25
Part 1 Spectral analysis
Sample x(t) Time series Population X(t)
26
Part 1 Spectral analysis
Sample x(t) Time series Population X(t) Pr
ocess
27
Part 1 Spectral analysis Process level
TIME DOMAIN
X(t)
28
Part 1 Spectral analysis Process level
TIME DOMAIN
X(t)
FOURIER TRANSFORMATION FREQUENCY DOMAIN
29
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
30
Part 1 Spectral analysis Process level
h(f) limT?8 E GT(f)2 / (2T)
NON-NORMALIZED POWER SPECTRAL DENSITY
FUNCTION, SPECTRUM
31
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
32
Part 1 Spectral analysis Process level
Discrete spectrum Harmonic process Astronomy
33
Part 1 Spectral analysis Process level
Discrete spectrum Harmonic process Astronomy
Continuous spectrum Random process Climatic noise
34
Part 1 Spectral analysis Process level
Discrete spectrum Harmonic process Astronomy
Continuous spectrum Random process Climatic noise
Mixed spectrum Typical climatic
35
Part 1 Spectral analysis
The task of spectral analysis is to estimate the
spectrum. There exist many estimation
techniques.
36
Part 1 Spectral analysis Harmonic regression
HARMONIC PROCESS
X(t) Sk Ak cos(2pfk t) Bk sin(2pfk t)
e(t)
37
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.
38
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
39
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.
40
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
41
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
42
Part 1 Spectral analysis Periodogram
43
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.
44
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.
45
Part 1 Spectral analysis Periodogram
A wonderful textbook Priestley MB
(1981) Spectral Analysis and Time
Series. Academic Press, London, 890 pp.
46
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.
47
Part 1 Spectral analysis Periodogram
48
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.
49
Part 1 Spectral analysis Smoothing
50
(No Transcript)
51
WELCH OVERLAPPED SEGMENT AVERAGING (WOSA)
1st Segment
2nd Segment
3rd Segment
52
WELCH OVERLAPPED SEGMENT AVERAGING (WOSA)
1st Segment
2nd Segment
3rd Segment
ERROR REDUCTION lt v3
53
Part 1 Spectral analysis Smoothing

Tapering Weight time series Spectral
leakage reduced (Hanning, Parzen, triangula
r windows, etc.)
54
Part 1 Spectral analysis Smoothing problem
Several segments averaged Spectrum estimate more
accurate -) Fewer (n lt n) data per
segment Lower frequency resolution -(
55
Part 1 Spectral analysis Smoothing problem
56
Part 1 Spectral analysis Smoothing problem
Subjective judgement is unavoidable. Play with
parameters and be honest.
57
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
58
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
59
Part 1 Spectral analysis 100-kyr problem

nd 650 kyr ?f (650 kyr)1
60
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
61
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
62
Part 1 Spectral analysis 100-kyr problem
The 100-kyr cycle existed not long enough to
allow a precise enough frequency estimation.
63
Part 1 Spectral analysis BlackmanTukey
h Fourier transform of ACV
64
Part 1 Spectral analysis BlackmanTukey
E X(t) X(t lag) h Fourier
transform of ACV
65
Part 1 Spectral analysis BlackmanTukey
PROCESS LEVEL E X(t) X(t lag) h
Fourier transform of ACV
66
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
67
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.
68
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.
69
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.
70
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
71
Part 1 Spectral analysis Multitaper Method
(MTM)
72
Part 1 Spectral analysis Multitaper Method
(MTM)
BETTER DIRECTLY VIA ASTRONOMY EQS.
73
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.
74
Part 1 Spectral analysis Further points
Uneven time spacing
75
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
76
Part 1 Spectral analysis Further points
PERSISTENCE
Red noise
77
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.
78
Part 1 Spectral analysis Further points
Aliasing
79
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.
80
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.
81
Part 1 Spectral analysis Further points

Detrending
82
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.
83
Part 1 Spectral analysis Further points
Bi-variate spectral analysis For example x
marine d18O y insolation
84
Part 1 Spectral analysis Further points
Higher-order spectra (bi-spectra, ...)
85
Part 1 Spectral analysis Further points
Etc., etc.
86
Part 2 Milankovic paleoclimate
87
Part 2 Milankovic paleoclimate
88
Part 2 Milankovic paleoclimate
Northern Hemisphere Glaciation NHG
89
Part 2 Milankovic paleoclimate
Mid-Pleistocene Transition
Northern Hemisphere Glaciation NHG
90
Part 2 Milankovic paleoclimate Climate
transitions, trend
91
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.
92
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
93
Part 2 Milankovic paleoclimate Mid-Pleistocen
e Transition
94
Part 2 Milankovic paleoclimate Mid-Pleistocen
e Transition
100 kyr cycle
95
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.
96
Part 2 Milankovic paleoclimate NHG
Database 24 Myr, 45 marine d18O records, 4
temperature records
benthic planktonic
Mudelsee M, Raymo ME (submitted)
97
NHGResults
High-resolution records
Mudelsee M, Raymo ME (submitted)
98
NHGResults
High-resolution records
Mudelsee M, Raymo ME (submitted)
99
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 .
100
Milankovic theory andtime series analysis
Conclusions
(1) Spectral analysis estimates
the spectrum. (2) Trend estimation is
also important (climate transitions).
101
G O O D I E S
102
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
103
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)
104
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