An Analysis of the Harmonic Properties of Global Ice Mass Cycles: Evidence from Deep Sea Core Sediment Records by Brian Michael Scholl

B.A. (Union College) 1996

A thesis submitted in partial satisfaction of the requirements for the degree of Master of Arts in Statistics in the GRADUATE DIVISION of the UNIVERSITY OF CALIFORNIA, BERKELEY

Committee in charge: Professor David R. Brillinger, Chair Professor David Freedman Professor Richard Muller

Spring 2003

The thesis of Brian Michael Scholl is approved:

Chair

Date

Date

Date

University of California, Berkeley

Spring 2003

An Analysis of the Harmonic Properties of Global Ice Mass Cycles: Evidence from Deep Sea Core Sediment Records

Copyright 2003 by Brian Michael Scholl

1

Abstract

An Analysis of the Harmonic Properties of Global Ice Mass Cycles: Evidence from Deep Sea Core Sediment Records by Brian Michael Scholl Master of Arts in Statistics University of California, Berkeley Professor David R. Brillinger, Chair

In this paper we examine the sample harmonic properties of time series constructed from several deep sea δ18 O sediment records under a constant sedimentation rate assumption and also the tuned records of Karner et al. (2002). Our estimates for Ocean Drilling Program Site 664 generally do not favor the Milankovitch hypothesis. The sample estimates of frequency for the northern hemisphere summer insolation series of [5], appear to correspond to the harmonics of precessional cycles. Some suggestion of obliquity is also observed. Sample harmonic properties of δ18 O time series records, however, appear to suggest a harmonic component with a period greater than 95 kyrs. The estimated coherence of Site 664 with the inclination series is greater than with summer insolation. In addition, the δ18 O data series, when band-pass filtered to preserve particular astronomical frequencies, appears to exhibit a much more stable relationship over time when the inclination component is preserved than when precession or obliquity frequencies are isolated.

2 Some support for inclination also exists among the eight other sites examined. The uncertainty in the mapping correspondence somewhat weakens these results, and all approximations are sensitive to site location and estimation method, but do not appear to support the insolation theory in any case, although we do not rule out the possibility that more complicated non-linear dynamics relating insolation to glacial accumulation may indeed be at work as discussed in [40].

Professor David R. Brillinger Thesis Committee Chair

iii

To David Brillinger my teacher, advisor, guide, To my Mom, Dad, Cathy McCarthy and the rest for their support, and to Ana, who has her own orbital properties.

iv

Contents List of Figures

v

List of Tables

vii

1 Introduction 2 Harmonic Analysis of Sediment Cores 2.1 Some Harmonic Properties of Astronomical Forcing Series 2.1.1 Inclination . . . . . . . . . . . . . . . . . . . . . 2.1.2 Insolation . . . . . . . . . . . . . . . . . . . . . . 2.2 Some Harmonic Properties of δ18 O Time Series . . . . . . 2.2.1 Untuned Site 664 Record . . . . . . . . . . . . . . 2.2.2 Tuned Site 664 Record . . . . . . . . . . . . . . . 2.2.3 Some Results for Remaining Sites . . . . . . . . . 2.2.4 Alternative Assumptions . . . . . . . . . . . . . . 2.3 Some Bivariate Relations . . . . . . . . . . . . . . . . . . 2.4 Chapter Summary . . . . . . . . . . . . . . . . . . . . . .

1

. . . . . . . . . .

6 8 8 12 15 18 23 26 31 37 40

3 Uncertainty in the Time-Depth Correspondence 3.1 Bivariate Relations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Complex Demodulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Dynamic Frequency Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

43 45 45 48

4 Concluding Remarks and New Directions

50

Bibliography

55

A

62

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

v

List of Figures 2.1

2.2 2.3

2.4 2.5 2.6 2.7

2.8

2.9 2.10

2.11

2.12 2.13

Inclination, TOP: Chiu’s R statistic with critical bound, MIDDLE: flattened periodogram with 95 percent and 99 percent confidence limits for flattened noise spectrum, BOTTOM: periodogram (dB) with bootstrapped upper confidence limit. . . . Bispectrum of inclination. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . June insolation 65’N, TOP: Chiu’s R statistic with critical bound, MIDDLE: flattened periodogram with 95 percent and 99 percent confidence limits for flattened noise spectrum, BOTTOM: periodogram (dB) with bootstrapped upper confidence limit. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Bispectrum of June insolation at 65 o N. . . . . . . . . . . . . . . . . . . . . . . . Tuned and untuned data, Site 664. . . . . . . . . . . . . . . . . . . . . . . . . . . Histogram of estimated sampling rate of the tuned δ18 O record of [30]. . . . . . . . Site 664, TOP: Chiu’s R statistic with critical bound, MIDDLE: flattened periodogram with 95 percent and 99 percent confidence limits for flattened noise spectrum, BOTTOM: periodogram (dB) with bootstrapped upper confidence limit. . . . LEFT: Site 664 (untuned), the filtered periodogram (frequency in radians), the bispectrum. RIGHT: sum-of-cosines residuals, filtered periodogram of residuals (frequency in radians), estimated signal (green) vs. time series (blue), bispectrum of residuals. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . BP filtered plots of Site 664 at various frequencies, untuned. . . . . . . . . . . . . Site 664 (tuned), TOP: Chiu’s R statistic with critical bound, MIDDLE: flattened periodogram with 95 percent and 99 percent confidence limits for flattened noise spectrum, BOTTOM: periodogram (dB) with bootstrapped upper confidence limit. Tuned record - LEFT: Site 664, the filtered periodogram (frequency in radians), the bispectrum. RIGHT: sum-of-cosines residuals, filtered periodogram of residuals (frequency in radians), estimated signal (green) vs. time series (blue), bispectrum of residuals. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . BP filtered plots of Site 664 at various frequencies, tuned. . . . . . . . . . . . . . . Residual analysis (Chiu’s R) of residuals from sum of cosines regression with bλ identified by periodogram peaks. . . . . . . . . . . . . . . . . . . . . . . . . . . .

9 11

13 14 16 17

19

20 22

24

25 27 32

vi 2.14 Results from the alternative order statistic of [47] - LEFT: Site 664 (untuned), the filtered periodogram (frequency in radians), the bispectrum. RIGHT: sum-of-cosines residuals, filtered periodogram of residuals (frequency in radians), estimated signal (green) vs. time series (blue), bispectrum of residuals. . . . . . . . . . . . . . . . . 2.15 Results from the alternative order statistic of [47] - LEFT: Site 664 (tuned), the filtered periodogram (frequency in radians), the bispectrum. RIGHT: sum-of-cosines residuals, filtered periodogram of residuals (frequency in radians), estimated signal (green) vs. time series (blue), bispectrum of residuals. . . . . . . . . . . . . . . . . 2.16 Residuals from AR fit to series . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.17 Spectrum estimates from associated AR fits to series (prewhitened). Frequency (lambda) is in cycles/kyr. TOP: untuned, BOTTOM: tuned. . . . . . . . . . . . . . 2.18 Phase and Coherence estimates for Insolation 65 o N and Site 664. Frequency is in cycles/kyr. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.19 Phase and Coherence estimates for Insolation 65 o N and Site 664. Frequency is in cycles/kyr. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 3.2 3.3

33

34 38 39 41 42

TOP ROW: coherency and phase estimates Sites 664 and 846; BOTTOM ROW: same figures Sites 664 and 982. These data are untuned and frequency is in cycles/kyr. 46 Complex demodulation (log amplitude) at λ=0.01, Site 664 . . . . . . . . . . . . . 47 Evolution of the periodogram maximizer (1/bλ). . . . . . . . . . . . . . . . . . . . 49

vii

List of Tables 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8

2.9

2.10 2.11

2.12

2.13

Spectrum Peaks (Chiu’s R Criteria): Inclination (β= 0.99 α= 0.95) . . . . . . . . . Spectrum Peaks (Chiu’s R Criteria): Insolation 65’N June (β= 0.99 α= 0.95) . . . . Spectrum Peaks (Chiu’s R Criteria): Insolation 15’S Jan (β= 0.99 α= 0.95) . . . . Spectrum Peaks (Chiu’s R Criteria): Site 664 (untuned) (β= 0.99 α= 0.95) . . . . . Results from QF technique, untuned Site 664. . . . . . . . . . . . . . . . . . . . . Spectrum Peaks (Chiu’s R Criteria): Site 664 (tuned) (β= 0.99 α= 0.95) . . . . . . Results from QF technique, Site 664 (tuned). . . . . . . . . . . . . . . . . . . . . Spectrum Peaks, and (last column) total number of peaks identified among the sites using the Benjamini-Hochberg multiple selection criteria at the given frequency. Untuned data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Spectrum Peaks, and (last column) total number of peaks identified among the sites using the Benjamini-Hochberg multiple selection criteria at the given frequency. Tuned data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . P-values for White test null hypothesis of linearity in mean . . . . . . . . . . . . . QF estimates of frequency (other parameter estimates provided in the supplemental appendix). 95% confidence interval obtained from the asymptotic distribution is provided as well as the mean and SD of the frequency estimate from 500 bootstrap simulations. K is estimated by the alternative statistic in [47]. Data is untuned. . . . QF estimates of frequency (other parameter estimates provided in the supplemental appendix). 95% confidence interval obtained from the asymptotic distribution is provided as well as the mean and SD of the frequency estimate from 500 bootstrap simulations. K is estimated by the alternative statistic in [47]. Data is tuned. . . . . AR estimate for site 664 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

A.1 Location and Approximate Depth (mcd) of Tie Points at Sites (NA indicates the depth record does not extend to the tie point). The last column is the average sedimentation rate (m/kyr) between the present and the last available tie point. This data provided by Jonathan Levine from the [30] paper. . . . . . . . . . . . . . . . . . .

8 12 14 18 21 25 26

29

30 31

36

36 37

63

viii

Acknowledgements David Brillinger, encouraged me to pursue this project at the outset even though the problems under investigation are of a very complicated nature. He always kept me motivated and prevented me from becoming discouraged. Along the way he suggested countless references, gave valuable comments, and provided guidance and wisdom beyond the scope of this paper. To David I will always be indebted. Richard Muller sparked an immediate interest in this problem during a fortuitous encounter at Berkeley’s computer store a few years ago. Rich managed to rope me into the topic with a few words in a matter of minutes. ”How on Earth can you get data for ancient climate patterns?” was the question that got it all started. A few weeks later I was desperately trying to find him on the Berkeley directory so that I could beg him for a peak. Jonathan Levine provided crucial assistance, provided much data and was always willing to talk over the course of the project. The BU Geology Department informally hosted me for a few weeks as I plundered their library and picked their brains; to Amy Dougherty, Rick Murry, and Mo Raymo I especially owe great thanks as well as Peter Huybers of MIT. David Freedman, Robert Shumway, Granville Tunnicliffe-Wilson and Bjorn Auestad provided helpful communication on their own work. Dan Karner also provided a wealth of helpful feedback on earlier drafts. And finally to my parents, family, friends and devojˇcica, for their love and support.

Brian Scholl December 3, 2002

1

Chapter 1

Introduction The widely accepted Milankovitch relationship between astronomical parameters and glacial cycles predicts that insolation (amount of solar radiation received) variations drive changes in global ice volume. Insolation has typically been associated with periodicities of approximately 41, 24, 22 and 19 thousand years (kyrs) (See [40]). Key variables used to calculate insolation are: eccentricity (the measure of orbital ellipticity), which has been associated with periods of approximately 400, 125 and 95 kyrs; obliquity (the tilt of the North Pole towards the Sun in mid-summer), which exhibits an approximately 41 kyr cycle; and precession (the lag between the Earth’s closest approach to the Sun, perihelion, and the Summer solstice), which exhibits periods of approximately 24, 22 and 19 kyrs. [40] provides more details. A cogent recent argument in support of astronomical forcing is provided by [37] based on the narrowness of observed peaks in the periodogram of climate proxy data obtained from deep sea sediment cores. However, ([39]) review shortcomings of the Milankovitch model, including the empirical absence of a 400 thousand year cycle in climate proxy records, though it dominates eccentricity variations. In turn, they relate climate variation to the inclination of the Earth’s orbit

2 through accretion of extraterrestrial particle matter. They predict a single 100 kyr cycle, which they observe in the deep-sea sediment core δ18 O data of the SPECMAP composite constructed by [28] and Site 806 of the Ocean Drilling Program (ODP) ([7]). ([37] presents similar results for additional core samples). The results are controversial and not widely accepted. There are at least three distinctive areas that may be of interest related to the statistical analysis of this topic: 1. evaluating the consistency of climate data with hypothesized astronomical forcing mechanisms given an assumed depth-time correspondence of the sediment core data record; 2. determining the sedimentation rate (i.e. the depth-time relation of the core record) using an astronomical template (this process is referred to as tuning); 3. joint determination of sedimentation rates and consistency with astronomical forcing mechanisms. The first of these items may include some examination of the harmonics of a tuned core record and an evaluation of the consistency of the data with the properties of insolation or inclination. The debate recently in the literature (e.g. [39] and [51]) has focused on re-evaluating the widely accepted astronomical forcing mechanism of insolation based on ostensibly more accurate time scale approximations. The second topic above is one of the original motivations for an interest in Milankovitch cycles during the past three decades - using the ”known” glacial forcing mechanism of Milankovitch to infer a sedimentation scale. Despite widespread criticism of tuning, early efforts were not made for the seemingly tautological purpose of detecting glacial cycles. Rather, early tuning exercises

3 were used to determine a sedimentation rate from a widely accepted climate forcing mechanism ([27], for example). The last topic has been largely neglected in the literature due to its apparent computational difficulty and with recent disagreements concerning the true astronomical forcing mechanism. However, we believe that a joint approach merits attention as the appropriate method for resolving the tuning and forcing problems: time scales are not easily approximated without some knowledge of the forcing mechanism, but the forcing mechanism is problematic due to the uncertainty of time scale approximations used to construct time series from sediment core records. Joint determination refers to the simultaneous address of the problems, rather than the sequential analysis to which the first and second problems above refer. In this paper we revisit the first problem above using δ18 O data obtained from fossilized benthic foraminifera provided to us by Jonathan Levine and Richard Muller. This data is used in [30] to construct the ”Benthic Stack” composite record and complete details, including site locations, are provided in that paper. We examine the sample harmonic properties of nine δ18 O records using untuned time scales (herein the untuned records) and the tuned records used to construct the Benthic Stack by [30]. δ18 O records have been widely employed in the paleoclimate literature (see [40] or [51]) and are generally accepted as a good proxy for historical climate features. In the present paper we will refer to δ18 O as a proxy for ice, but note that other climate features such as temperature are also somewhat reflected by the δ18 O data. Of the fourteen core samples originally considered by the authors of [30], we focused on the ten cores for which the Brunhes-Matuyama boundary (BM) was identified by those authors (please refer to Table A.1) because of a desire to include only the records containing data for the

4 entire period under consideration.1 To construct the ”untuned” records, we mapped the original, irregularly sampled, composite core records (denominated in meters composite depth (mcd)) to the time domain using the average rate of sedimentation between the first and the last tie point (the present and the BM boundary). Using a cubic spline interpolation, these irregularly sampled time series were then converted to regular time series with one observation point every kyr. In most cases, the interpolation appeared to be reasonable, but one site was dropped from the subsequent analysis due to an unreasonable fit provided by the spline interpolation function. These same nine site locations were selected from those tuned by [30] and that data will herein be referred to as the ”tuned” data. The tuned data is also interpolated to equally spaced time intervals. Our objective in this paper is to detect frequencies of interest in the regular time series constructed from the core samples and examine them in the context of hypothesized astronomical forcing series. There are limitations to the examination of δ18 O data as regular time series because of the sampling properties of the depth-core records that these are constructed from, but these problems are not explored in the current paper. The tuning and joint determination problems mentioned above are not investigated directly in the current paper, although we present a sketch of some potentially interesting future research along these lines in our concluding chapter. Our estimates for Ocean Drilling Program Site 664 generally do not favor the Milankovitch hypothesis. The sample estimates of frequency for the northern hemisphere summer insolation series of [5], appear to correspond to the harmonics of precessional cycles. Some suggestion of obliquity is also observed. Sample harmonic properties of δ18 O time series records, however, appear to suggest a harmonic component with a period greater than 95 kyrs. The estimated coherence of Site 1 Three

additional cores were later examined for eventual publication in [30], but these were not available to us at the time of our undertaking.

5 664 with the inclination series is greater than with summer insolation. In addition, the δ18 O data series, when band-pass filtered to preserve particular astronomical frequencies, appears to exhibit a much more stable relationship over time when the approximately 100 kyr harmonic is preserved than when precession or obliquity frequencies are isolated. Some support for inclination also exists among the eight other sites examined. The uncertainty in the mapping correspondence somewhat weakens these results, and all approximations are sensitive to site location and estimation method, but do not appear to support the insolation theory in any case, although we do not rule out the possibility that more complicated non-linear dynamics relating insolation to glacial accumulation may indeed be at work as discussed in [40]. This paper proceeds as follows. Chapter 2 provides an analysis of the harmonic properties of untuned and tuned δ18 O records from nine sample sites. Chapter 3 examines the uncertainty of the depth-time, or mapping, correspondence and provides some general criteria that may be considered when evaluating a mapping function as well as the related issue of the alignment of multiple time series. Chapter 4 summarizes our findings and sketches some new directions planned for future work.

6

Chapter 2

Harmonic Analysis of Sediment Cores In this chapter we provide a harmonic analysis of δ18 O records obtained from deep sea sediment cores as described in [30]. We have available δ18 O data extracted from i=1..r=9 deep sea cores at measured depths li1 , li2 , ...L, which we denote by the discrete data set yi (l j ). In making inferences relevant for glacial forcing, it is often desirable to analyze the series yi (τi (l j )), with τi (l j ) the assumed time associated with the j-th depth point of the i-th core sample for some determined depth-time, or mapping, correspondence such as that provided by assuming a constant sedimentation rate to the BrunhesMatuyama Boundary or the tuning work of [30]. As the astronomical forcing mechanisms under consideration are assumed to vary cyclically in some way determined by regular astronomical cycles, a pertinent model of the hybrid process of the i-th core record is provided by:

K

yi (τi (l j )) =

∑ ρi,k cos(λi,k τi (l j ) + φi,k ) + εi (t).

(2.1)

k=1

Henceforth, we will drop the subscript i and use t to denote time. This not only simplifies the notation, but also keeps the presentation of this chapter in the notation of the regularly sampled time-series that has been constructed for analysis by an interpolation of points on the depth

7 series. (The analysis of δ18 O records as marked point processes provides a number of interesting possibilities that we discuss briefly in our concluding chapter, but we have chosen to focus on the time series data and do not explore these possibilities in the current paper). With these notational simplifications, the model under consideration is:

K

y(t) =

∑ ρk cos(λkt + φk ) + ε(t).

(2.2)

k=1

Our prime interest in the analysis of this chapter is to provide sample estimates of the parameters K and λ1 ...λK determining the number of periodic components and the corresponding frequencies, respectively. We will also examine the agreement of our δ18 O time series with insolation and inclination forcing data. For a discussion of the sum-of-cosine modelling employed in this chapter, the reader is referred to [14], [15] and [47]. An approximation of K can be obtained by several means, but such methods rely on assumptions relating to the non-signal or noise process ε(t). One might imagine a likelihood ratio test to determine this parameter. However, even in the simplest case, [47] demonstrate the difficulty in constructing a test. Under the assumption that the non-signal process is gaussian white noise, the test statistic RT (β) developed by [23] provides a means of determining the frequencies of interest, and in turn K, from the periodogram. In [47], several information criterion statistics are developed in the spirit of the AIC statistic that has been employed to determine the order of an autoregression. These inference statistics, which evaluate the residuals from a sum of cosines regression of the form in 2.2, follow from [46] and subsequent work. [47] develop alternative statistics under the assumption that ε(t) is a process of unknown character. We will begin with the assumption that the residual process is stationary gaussian white

8

1 2 3 4 5 6

bλ 0.00125 0.00500 0.01000 0.01250 0.01625 0.02000

95% CI (0.001153, 0.001347) (0.004956, 0.005044) (0.009998, 0.010002) (0.012467, 0.012533) (0.016067, 0.016433) (0.019988, 0.020012)

1/bλ 800.0 200.0 100.0 80.0 61.5 50.0

Amp (dB) 0.8319 6.1028 19.6545 7.2217 −0.4324 0.5797

p-value 0.00 0.00 0.00 0.00 0.00 0.00

test power 0.98 0.99 1.00 0.99 0.98 0.98

Table 2.1: Spectrum Peaks (Chiu’s R Criteria): Inclination (β= 0.99 α= 0.95)

noise. A Phillips-Perron test rejects the null hypothesis of non-stationarity for all of the data series under consideration with a p-value of 0.01 or less.[43]

2.1 Some Harmonic Properties of Astronomical Forcing Series In this section we examine some basic harmonic properties of inclination and insolation, the two most prominent potential astronomical forcing mechanisms for climate.

2.1.1 Inclination Inclination has received attention of late as a possible driving force for climate cycles due to the important recalculation of [39] with respect to the invariable plane. Inclination has been associated with harmonic periods of 100 kyrs. Figure 2.3 provides plots helpful to the understanding of inclination harmonics. The results presented here are indicative of our method of investigation throughout this chapter. The panels present Chiu’s R statistic with critical region for α = .95, β =.99 [23], the flattened periodogram suggested by [15] with 95 (lower line) and 99 (upper line) percent upper confidence level for a flattened noise spectrum, and the periodogram with upper 95 percent confidence region obtained by local bootstrap as per [42]. Where relevant, the data has been tapered by a cosine bell taper prior to

9

0 −40 −100

Chiu’s R (dB)

Inclination

0.0

0.1

0.2

0.3

0.4

0.5

0.4

0.5

0.4

0.5

0

2

4

6

frequency (cycles/kyr)

0.0

0.1

0.2

0.3

−20 20 −80

dB

frequency (cycles/kyr)

0.0

0.1

0.2

0.3

frequency (cycles/kyr)

Figure 2.1: Inclination, TOP: Chiu’s R statistic with critical bound, MIDDLE: flattened periodogram with 95 percent and 99 percent confidence limits for flattened noise spectrum, BOTTOM: periodogram (dB) with bootstrapped upper confidence limit.

10 estimation (see [60], p 411). (Other tapers were considered during the analysis of individual orbital and δ18 O data, but although these may alter the periodogram estimate at a particular frequency, there was little difference in the identification of peak frequencies; that is to say the estimates of the parameter λ were not greatly affected). We define peaks in the periodogram as local maxima that exceed the critical value determined by Chiu’s test. There are indications of an approximately 100 kyr harmonic in the series. The periodogram maximizer is at approximately 100 kyr, with 95 percent confidence intervals given for the frequency estimates in Table 2.1 as per the distributional results presented in ([47], Section 3.3). For the 100 kyr cycle, the p-value associated with the null hypothesis that the peak is produced by a gaussian white process against the alternative that the frequency represents a cosinusoidal signal, according to Chiu’s [23] method, is less than 0.01, with the power of the test rounding to 1.00. Significant peaks are also detected by the periodogram corresponding to the approximate periods of 800, 200, 80, 61.5 and 50 thousand years.1 The flattened periodogram suggests an interest in the frequencies with periods of approximately 100, 50, 33, 11 kyrs, as well as some other higher frequency components. In our work, we are not interested in very high frequency components because these are unlikely to be important for the climate process, which is assumed to have a low-frequency response to forcing models. The maximum value of the periodogram corresponds to the 100 kyr cycle. Without a detailed understanding of the distribution of the noise, the frequency estimate corresponding to the periodogram maximum cannot be bettered in terms of asymptotic variance ([46]). However, the Quinn and Fernandes (QF) method treats 2.2 as a regression equation and provides an asymptotically 1 At first glace the p-values reported for p-gram peaks may seem incredulous (similar low values are reported through-

out this chapter). However, these low values are correct; they are attributable to the steep slope of Chiu’s distribution function and the fact that the amplitude of peaks relative to non-peak frequencies is large.

11

Bispectrum of Series

0.09

Frequency (cycles/kyr)

0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01 6.37e+012 0

6.37e+012 0

0.02

0.04 0.06 Frequency (cycles/kyr)

0.08

Figure 2.2: Bispectrum of inclination.

equivalent alternative estimate ([47] and [44]). The QF estimate of frequency, with K estimated by the alternative estimator in [47] that does not depend on a detailed understanding of the noise process, estimates harmonics (and approximate 95 percent confidence intervals) at bλ=.01 (0.0014, 0.0185) and .0007 (-0.0032,0.0045). These frequencies correspond to approximate 100 and 1514 kyr cycles. The general difficulty in estimating low frequencies accurately is discussed in [47].2 Finally, we note that the bispectrum of the series (Figure 2.2) is suggestive of nonlinearities. 2 Note that approximate confidence interval estimates that include zero do not imply a rejection of the frequency estimate; λ=0 is not a relevant null hypothesis for frequency estimation as discussed in [47].

12

1 2 3 4

bλ 0.02500 0.04250 0.04500 0.05250

95% CI (0.024998, 0.025002) (0.042497, 0.042503) (0.044994, 0.045006) (0.052496, 0.052504)

1/bλ 40.0 23.5 22.2 19.0

Amp (dB) 41.3794 47.7483 45.3400 43.9275

p-value 0.00 0.00 0.00 0.00

test power 1.00 1.00 1.00 1.00

Table 2.2: Spectrum Peaks (Chiu’s R Criteria): Insolation 65’N June (β= 0.99 α= 0.95)

2.1.2 Insolation The key mechanism of the Milankovitch theory is summer insolation in the northern hemisphere ([40], hence; we examine the insolation for 65 degrees North in June calculated by [5]. Eccentricity, obliquity and precession are the primary variables used to calculate insolation. Our results (Table 2.2) are suggestive of periodic components corresponding to the frequencies typically associated with precession (approximately 19, 22.2 and 23.5 kyrs) and to a slightly lesser extent, obliquity (approximately 40 kyrs) (see [40]). The QF technique estimates are bλ=.0443 (0.0432,0.0453) and 0.0423 (0.0414, 0.0431) (again, K is estimated by the alternative information criterion in [47]). The QF estimates correspond to periods of approximately 22.57 and 23.64 kyrs. Results differ somewhat for insolation at different latitudes and seasons:3 January insolation at 65 o S suggests another peak at approximately 53.3 kyrs; June insolation at 15 o N does not show an obliquity peak; and January insolation at 15 o S suggests several additional peaks including ones at approximately 400 and 100 kyrs (Table 2.3). There are some suggestions of non-linearities in the series as per the contours exhibited in the bispectrum, Figure 2.4. 3 Complete

tabulation of all results in this chapter are provided in an appendix available on request from the author.

13

−60 −120

Chiu’s R (dB)

0

Insolation 65’N June

0.0

0.1

0.2

0.3

0.4

0.5

0.4

0.5

0.4

0.5

0

2

4

6

8

frequency (cycles/kyr)

0.0

0.1

0.2

0.3

−20 −80

dB

40

frequency (cycles/kyr)

0.0

0.1

0.2

0.3

frequency (cycles/kyr)

Figure 2.3: June insolation 65’N, TOP: Chiu’s R statistic with critical bound, MIDDLE: flattened periodogram with 95 percent and 99 percent confidence limits for flattened noise spectrum, BOTTOM: periodogram (dB) with bootstrapped upper confidence limit.

14

Bispectrum of Series

0.09

Frequency (cycles/kyr)

0.08 0.07 0.06 0.05 0.04 0.03

3.27e+019 3.27e+019 3.27e+019

0.02 0.01 0

0

0.02

0.04 0.06 Frequency (cycles/kyr)

0.08

Figure 2.4: Bispectrum of June insolation at 65 o N.

1 2 3 4 5 6 7 8 9 10 11

bλ 0.00250 0.01000 0.01750 0.02000 0.02500 0.02875 0.04250 0.04500 0.05250 0.06125 0.06875

95% CI (0.002495, 0.002505) (0.009994, 0.010006) (0.017489, 0.017511) (0.019994, 0.020006) (0.024991, 0.025009) (0.028719, 0.028781) (0.042497, 0.042503) (0.044994, 0.045006) (0.052496, 0.052504) (0.061115, 0.061385) (0.068737, 0.068763)

1/bλ 400.0 100.0 57.1 50.0 40.0 34.8 23.5 22.2 19.0 16.3 14.5

Amp (dB) 19.5164 17.2861 13.7369 14.6647 12.8726 11.8990 48.0348 45.6352 44.2236 22.6269 15.6717

p-value 0.00 0.00 0.00 0.00 0.00 0.04 0.00 0.00 0.00 0.00 0.00

test power 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00

Table 2.3: Spectrum Peaks (Chiu’s R Criteria): Insolation 15’S Jan (β= 0.99 α= 0.95)

15

2.2 Some Harmonic Properties of δ18 O Time Series In this section we examine some harmonic properties of both untuned and tuned δ18 O time series records. Again, the tuned records are those used to construct the Benthic Stack in [30] and for complete details of the data and a thorough discussion of the justification of the use of δ18 O records as proxies for climate, the reader is referred to [30] and [40], respectively. For ease of discussion concerning the large amount of data considered (a total of 18 time series), we will present results in detail for Site 664 and summarize some results for additional sites. Although results do differ between individual drilling locations, no particular reason exists for focusing on Site 664. Note that frequency units are in cycles per 1000 years, all records have been de-meaned prior to analysis and, when appropriate, tapered as described in the prior section. A plot of the regular time series for Site 664 appears in Figure 2.5 (data on the right is closer to the present). The tuned series only differs slightly from the untuned series. Most noticeably, the large peaks in the data between 100 and 500 kyrs before the present are shifted back slightly. Both series exhibit occasional abrupt glacial terminations (a rapid increase in the value after a local minimum). Before continuing, the following caveat is provided. Regarding both the untuned and tuned time series records, the reader is strongly cautioned against making conclusions about results for high frequencies (higher than about 0.1 cycles/kyr). With a given time-scale assumption, the sampling properties of the original data set makes resolution of high frequency properties extremely problematic. The average sampling rate is estimated at approximately 3.4 kyrs for the tuned records of [30] and a plot of the sampling distribution is provided in Figure 2.6. A thorough treatment of this issue, however, is complex and beyond the scope of this paper.

16

−2.5

−3

−3.5

−4

−4.5

−5

−5.5 790

tuned untuned

700

600

500 400 300 200 age (kyrs before present)

Figure 2.5: Tuned and untuned data, Site 664.

100

0

17

140

Number of observations

120

100

80

60

40

20

0

0

5

10 15 20 Change in kyrs per sampling point

25

Figure 2.6: Histogram of estimated sampling rate of the tuned δ18 O record of [30].

18

1 2 3 4 5 6 7

bλ 0.00125 0.00625 0.01000 0.01875 0.02250 0.02500 0.03250

95% CI (0.001229, 0.001271) (0.006218, 0.006282) (0.009996, 0.010004) (0.018698, 0.018802) (0.022489, 0.022511) (0.024994, 0.025006) (0.032492, 0.032508)

1/bλ 800.0 160.0 100.0 53.3 44.4 40.0 30.8

Amp (dB) 2.1360 3.5385 12.8820 0.6159 3.4081 6.1102 4.0703

p-value 0.00 0.00 0.00 0.01 0.00 0.00 0.00

test power 0.99 0.99 1.00 0.98 0.99 0.99 0.99

Table 2.4: Spectrum Peaks (Chiu’s R Criteria): Site 664 (untuned) (β= 0.99 α= 0.95)

2.2.1 Untuned Site 664 Record An evaluation of Chiu’s R statistic (β = .99, α = .95) suggests several peaks for which we reject the null hypothesis of no periodic component ([23]). Figure 2.7 plots Chiu’s statistic and approximate frequencies associated with the peaks, confidence limits and other data are provided in Table 2.4. Periodic elements are suggested with periods of approximately 800, 160, 100, 53, 44, 40 and 31 thousand years. Flattened periodogram estimates suggest an interest in the frequencies corresponding to approximately 100, 89, 40, 31, 23, 18.6 17.3, 13.5, 11.7 and 10.2 thousand year cycles as well as other higher frequency harmonics. The maximum value of the periodogram estimate, 12.88 dB, corresponds to the 100 kyr cycle and a second large value of the periodogram coincides with the 40 kyr cycle. The results from the asymptotically equivalent QF estimator are b summarized in Table 2.5.4 K=11 frequencies are approximated by the order selection statistic of [46], which relies on an assumption that the residuals are white noise. An approximately 100 kyr signal is estimated along with an extremely low frequency estimate of high amplitude. Other harmonics are also detected, but for all of these estimates except the low frequency one, the associated standard errors (from [47]) are wide. The filtered periodogram of [47], which is related to the QF 4 We

results.

examined several simple modifications to the series including linearly detrending the series, with little change in

19

−40 −80

Chiu’s R (dB)

0

Site 664

0.0

0.1

0.2

0.3

0.4

0.5

0.4

0.5

0.4

0.5

0

1

2

3

4

frequency (cycles/kyr)

0.0

0.1

0.2

0.3

−20 −60

dB

frequency (cycles/kyr)

0.0

0.1

0.2

0.3

frequency (cycles/kyr)

Figure 2.7: Site 664, TOP: Chiu’s R statistic with critical bound, MIDDLE: flattened periodogram with 95 percent and 99 percent confidence limits for flattened noise spectrum, BOTTOM: periodogram (dB) with bootstrapped upper confidence limit.

20

Time Series 1

0

0

−2 790 1

kyr b.p.

0

−1

−1 790 0.1

0

κ

κ

Residuals from QF Method

2

0

1

2

3

0

−0.1

4

0

2

0

1

2

3

4

0

3.4e+011

0.01 0

0

0.005 0.01 0.015 cycles/kyr

cycles/kyr

cycles/kyr

−2 790

0

0.01 3.97e+007 0

0

0.005 0.01 0.015 cycles/kyr

Figure 2.8: LEFT: Site 664 (untuned), the filtered periodogram (frequency in radians), the bispectrum. RIGHT: sum-of-cosines residuals, filtered periodogram of residuals (frequency in radians), estimated signal (green) vs. time series (blue), bispectrum of residuals. estimate, suggests an interest in the periodic components with 392.5, 98.1, and 41.3 thousand year cycles. One might notice that the power of the filtered periodogram of the residual series (second panel on right) is substantially reduced from the original series. Also, the estimated signal seems to track the time series well - the estimated signal and series are almost indistinguishable in the third panel. Yet, the bispectrum of the residual series seems to suggest nonlinearities, although a White test of the null hypothesis of linearity in the mean of the series has a p-value of 0.3426115.

21 bλ 0.009026797 0.011512875 0.024447071 0.009982538 0.033075372 0.006444168 0.014309918 0.022276237 4.30E-07 0.026906079 0.018845661

b ρ 0.01867 0.02660 0.02005 0.03659 0.01336 0.00967 0.01382 0.01114 3.52E+04 0.00798 0.00710

CI (λ) ( -0.01258 , 0.03064 ) ( -0.00075 , 0.02378 ) ( -0.00048 , 0.04937 ) ( 0.00438 , 0.01559 ) ( -0.02544 , 0.09159 ) ( -0.30682 , 0.31971 ) ( -0.02961 , 0.05823 ) ( -0.05950 , 0.10406 ) ( 4.30E-07 , 4.30E-07 ) ( -0.11097 , 0.16478 ) ( -0.14177 , 0.17946 )

SE (BS) 0.0011 0.0005 0.0052 0.0079 0.0122 0.0114 0.0083 0.0088 0.0128 0.0191 0.0209

1/bλ 110.8 86.9 40.9 100.2 30.2 155.2 69.9 44.9 2.3E+06 37.2 53.1

Table 2.5: Results from QF technique, untuned Site 664.

Some Filtered Series Figure 2.9 plots the untuned record for Site 664 after isolating particular frequency ranges with a two-sided band-pass filter as described in [24]. The band-pass filter window widths were selected to isolate astronomical forcing frequencies related to the approximately 100 kyr cycle (associated primarily with inclination, but also with eccentricity), obliquity, and precession. The last panel is an isolation of the high frequency component obtained by preserving the frequencies with periods between 2 and 11 kyrs. Perhaps the most interesting aspect of these figures is the stability of the 100 kyr signal over the sample period as compared with the other filtered series. The variance of the obliquity signal is much larger for the most recent two-thirds of the sample. The precession component, which is the dominant mechanism determining insolation, is quite frequently changing in character over the sample. The results appear inconsistent with the insolation model because any astronomical forcing mechanism over the given sample period is expected to have a reasonably stable effect on climate.

22

100 kyr signal

Obliquity signal

0.4

0.4

0.2

0.2

0

0

−0.2

−0.2

−0.4 790 0 kiloyears before present Precession signal 0.4

−0.4 790 0 kiloyears before present noise signal 1

0.2

0.5

0

0

−0.2

−0.5

−0.4 790 0 kiloyears before present

−1 790 0 kiloyears before present

Figure 2.9: BP filtered plots of Site 664 at various frequencies, untuned.

23

2.2.2 Tuned Site 664 Record The band-pass filter tuning of [30] provides quite similar results to those of the untuned record. Turning again to Chiu’s statistic (Figure 2.10), peaks in the periodogram are detected with the same frequencies as before except that the approximately 31 kyr cycle does not appear and an approximate 36 and 30 kyr cycles are detected. Results are summarized in Table 2.6. Confidence intervals for the frequency values differ from the estimates for the untuned record because the amplitude of the peaks has changed. In Figure 2.10 the flattened periodogram suggests an interest in cycles repeating at approximately 100, 88.9, 40, 29.6, 18.6, 17.3, 14.5, and 10.3 thousand years as well as other high frequency cycles. The magnitude of the 40 kyr peak is greatest; the generally stronger amplitude of an approximate 40 kyr cycle in the tuned record is unsurprising because the tuning was designed to intensify the obliquity component in the untuned record. A single peak may be identified in the filtered periodogram with a period of approximately b=8 98.4 kyrs. Table 2.7 displays the results from the QF estimation of the frequency parameter. K frequencies are estimated. Though confidence intervals are wide for many of the frequency parameters, two high amplitude cycles are estimated with frequencies corresponding to approximately 106.4 kyrs. The high amplitude accounts for the narrow confidence intervals of these estimates. The estimated signal from the QF method seems to track the time series reasonably well, although much high frequency variation is not predicted by the model (Figure 2.11). Some suggestion of nonlinearity may be seen from the bispectrum and the peaks in the filtered periodogram of the residual series from the QF method although the p-value of a White test null hypothesis of linearity in the mean is 0.342, therefore we cannot reject the null hypothesis of linearity.

Filtered

24

−10 −40 −70

Chiu’s R (dB)

Site 664, (tuned)

0.0

0.1

0.2

0.3

0.4

0.5

0.4

0.5

0.4

0.5

0

1 2

3

4

frequency (cycles/kyr)

0.0

0.1

0.2

0.3

−20 0 −50

dB

frequency (cycles/kyr)

0.0

0.1

0.2

0.3

frequency (cycles/kyr)

Figure 2.10: Site 664 (tuned), TOP: Chiu’s R statistic with critical bound, MIDDLE: flattened periodogram with 95 percent and 99 percent confidence limits for flattened noise spectrum, BOTTOM: periodogram (dB) with bootstrapped upper confidence limit.

25

Time Series

Residuals from QF Method

2

1

0

0

0

−1

−1 790 0.1

0

κ

κ

−2 790 1

0

1

2

3

0

−0.1

4

0

0

1

2

3

4

2 0

2.92e+011 0.01 0

0

0.005

cycles/kyr

cycles/kyr

−2 790

0 2.07e+008

0.01 0

0.01 0.015 cycles/kyr

0

0.005

0.01 0.015 cycles/kyr

Figure 2.11: Tuned record - LEFT: Site 664, the filtered periodogram (frequency in radians), the bispectrum. RIGHT: sum-of-cosines residuals, filtered periodogram of residuals (frequency in radians), estimated signal (green) vs. time series (blue), bispectrum of residuals.

1 2 3 4 5 6 7 8

bλ 0.00125 0.00625 0.01000 0.01875 0.02250 0.02500 0.02750 0.03375

95% CI (0.001232, 0.001268) (0.006216, 0.006284) (0.009996, 0.010004) (0.018687, 0.018813) (0.022489, 0.022511) (0.024996, 0.025004) (0.027485, 0.027515) (0.033744, 0.033756)

1/bλ 800.0 160.0 100.0 53.3 44.4 40.0 36.4 29.6

Amp (dB) 2.7956 3.1707 11.9212 −0.2911 2.8787 7.0396 1.5878 3.6033

p-value 0.00 0.00 0.00 0.01 0.00 0.00 0.00 0.00

test power 0.99 0.99 1.00 0.98 0.99 0.99 0.99 0.99

Table 2.6: Spectrum Peaks (Chiu’s R Criteria): Site 664 (tuned) (β= 0.99 α= 0.95)

26 bλ 0.00940 0.01130 0.00940 0.02465 0.03329 0.00680 0.01446 0.02286

b ρ 41.59971604 0.320724733 41.68047996 0.178509002 0.142866718 0.093526688 0.149176778 0.110039861

CI (λ) ( 0.00939 , 0.00940 ) ( -0.00654 , 0.02915 ) ( 0.00940 , 0.00940 ) ( -0.02796 , 0.07726 ) ( -0.01997 , 0.08655 ) ( -0.01052 , 0.02411 ) ( -0.05789 , 0.08682 ) ( -0.12626 , 0.17198 )

SE (BS) 0.0008 0.0005 0.0065 0.0092 0.0113 0.0111 0.0078 0.0085

1/bλ 106.4 88.5 106.4 40.6 30.0 147.1 69.1 43.7

Table 2.7: Results from QF technique, Site 664 (tuned).

series (Figure 2.12) exhibit similar characteristics to those reported for the untuned series. The obliquity and precessional signals do not appear stable over the sample period. Some modulation in the 100 kyr signal is observed as well.

2.2.3 Some Results for Remaining Sites Tables 2.8 and 2.9, respectively, summarize the results obtained from an application of Chiu’s R statistic to the nine sample sites for the untuned and the tuned records. The final column in each table indicates the number of rejections of the null hypothesis of no cosinusoidal component according to the sequential comparison method of [4]. The procedure is designed to control the expected False Discovery Rate (FDR) and guarantees that E(FDR)≤

M m (1 − α)

≤ (1 − α), where

M is the number of hypotheses that are true (unobservable) and m is the total number of tests. In the present context m=9. α has been set at the 95 percent confidence level. In Table 2.9, one may notice that an approximately 72 kyr peak is implied for the tuned series (i.e. peaks=9=m), but that no such peak was reported above. Previously, our tables recorded only peaks - that is, local maxima - determined significant by Chiu’s criteria; in the present summary, all frequencies for which the null hypothesis of no cosinusoid is rejected are tested in the sequential comparison method. This

27

100 kyr signal

Obliquity signal

0.2

0.4

0.1

0.2

0

0

−0.1

−0.2

−0.2 790 0.4

0 kiloyears before present Precession signal

−0.4 790 0.6

0 kiloyears before present noise signal

0.4

0.2

0.2 0 0 −0.2

−0.2

−0.4 790

0 kiloyears before present

−0.4 790

0 kiloyears before present

Figure 2.12: BP filtered plots of Site 664 at various frequencies, tuned.

28 appears to be a reasonable method of inclusion both by the implications of Chiu’s test framework (which rejects the null hypothesis for periodogram ordinates in excess of the critical value) and also bearing in mind that mapping imperfections may adjust estimates of peak frequency in some sample sites (we return to this issue in the next chapter). For the untuned records, the null hypothesis of no cosinusoidal signal is rejected for all 9 sites for the approximately 114.3 thousand year harmonic. Again, this harmonic was not reported earlier because it is merely a sidelobe for Site 664; indeed only Site 502, 552 and 982 report a peak at the corresponding frequency. This is the only cycle that is detected for all sites. Eight null hypotheses are rejected for the periods corresponding to approximately 160, 89 and 73 thousand years. Somewhat interesting is the rejection of the null hypothesis at seven of the sites for 800 and 200 thousand year harmonics. Such harmonics have not previously been reported by other authors. Frequencies normally associated with precession, in the range of approximately 17-24 kyrs, are almost non-existent. Cycles consistent with obliquity (typically associated with approximately 4044 kyr cycles in the literature) are observed at many of the sites. Two sites do not reject the null hypothesis for the 100 kyr harmonic, which is possibly an indication that mapping imperfections have shifted the peak to the 114 kyr frequency for these sites. The tuned records exhibit broadly similar results. All sites reject the null hypothesis for harmonics of approximately 114, 100 and 73 kyrs. Some support remains for the 800, 200 and 160 kyr cycles. Precessional components remain largely absent. There is some detection of obliquity cycles, but again, this is to be expected given that the records were tuned to match the 41 kyr obliquity component ([30]).

29

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32

bλ 0.0013 0.0025 0.0037 0.0050 0.0063 0.0075 0.0087 0.0100 0.0113 0.0125 0.0138 0.0150 0.0163 0.0175 0.0188 0.0200 0.0213 0.0225 0.0238 0.0250 0.0263 0.0275 0.0288 0.0300 0.0313 0.0325 0.0338 0.0350 0.0375 0.0388 0.0413 0.0475

95% CI (0.0012, 0.0013) (0.0025, 0.0025) (0.0037, 0.0038) (0.0050, 0.0050) (0.0062, 0.0063) (0.0075, 0.0075) (0.0087, 0.0088) (0.0100, 0.0100) (0.0112, 0.0113) (0.0125, 0.0125) (0.0137, 0.0138) (0.0150, 0.0150) (0.0162, 0.0163) (0.0175, 0.0175) (0.0187, 0.0188) (0.0200, 0.0200) (0.0212, 0.0213) (0.0225, 0.0225) (0.0237, 0.0238) (0.0250, 0.0250) (0.0262, 0.0263) (0.0275, 0.0275) (0.0287, 0.0288) (0.0300, 0.0300) (0.0312, 0.0313) (0.0325, 0.0325) (0.0337, 0.0338) (0.0350, 0.0350) (0.0375, 0.0375) (0.0387, 0.0388) (0.0412, 0.0413) (0.0475, 0.0475)

1/bλ 800.0 400.0 266.7 200.0 160.0 133.3 114.3 100.0 88.9 80.0 72.7 66.7 61.5 57.1 53.3 50.0 47.1 44.4 42.1 40.0 38.1 36.4 34.8 33.3 32.0 30.8 29.6 28.6 26.7 25.8 24.2 21.1

BH Nulls Rejected 7 2 6 7 8 6 9 7 8 2 8 6 3 5 6 1 4 6 3 5 2 5 2 2 1 2 2 2 3 1 1 1

Table 2.8: Spectrum Peaks, and (last column) total number of peaks identified among the sites using the Benjamini-Hochberg multiple selection criteria at the given frequency. Untuned data.

30

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28

bλ 0.0013 0.0025 0.0037 0.0050 0.0063 0.0075 0.0087 0.0100 0.0113 0.0125 0.0138 0.0150 0.0163 0.0175 0.0188 0.0200 0.0213 0.0225 0.0238 0.0250 0.0263 0.0275 0.0300 0.0325 0.0338 0.0350 0.0388 0.0425

95% CI (0.0012, 0.0013) (0.0025, 0.0025) (0.0037, 0.0038) (0.0050, 0.0050) (0.0062, 0.0063) (0.0075, 0.0075) (0.0087, 0.0088) (0.0100, 0.0100) (0.0112, 0.0113) (0.0125 0.0125) (0.0137, 0.0138) (0.0150, 0.0150) (0.0162, 0.0163) (0.0175, 0.0175) (0.0187, 0.0188) (0.0200, 0.0200) (0.0212, 0.0213) (0.0225, 0.0225) (0.0237, 0.0238) (0.0250, 0.0250) (0.0262, 0.0263) (0.0275, 0.0275) (0.0300, 0.0300) (0.0325, 0.0325) (0.0337, 0.0338) (0.0350, 0.0350) (0.0387, 0.0388) (0.0425, 0.0425)

1/bλ 800.0 400.0 266.7 200.0 160.0 133.3 114.3 100.0 88.9 80.0 72.7 66.7 61.5 57.1 53.3 50.0 47.1 44.4 42.1 40.0 38.1 36.4 33.3 30.8 29.6 28.6 25.8 23.5

BH Nulls Rejected 7 1 6 7 8 6 9 9 7 4 9 3 2 5 4 4 5 7 5 5 2 6 1 4 1 1 2 1

Table 2.9: Spectrum Peaks, and (last column) total number of peaks identified among the sites using the Benjamini-Hochberg multiple selection criteria at the given frequency. Tuned data.

31

502 552 607 659 664 846 849 925 982

Untuned yt 0.8415716 0.1451898 0.8993135 0.8005026 0.3426115 0.4905442 0.9821532 0.1030922 0.4159219

∆yt 5.87E-08 1.31E-03 0.00E+00 1.54E-04 3.73E-03 5.56E-06 8.44E-03 3.56E-08 3.99E-03

Tuned yt 0.839885 0.9014542 0.8742579 0.4691265 0.3420924 0.9972838 0.7405987 0.9763824 0.2279935

∆yt 1.53E-04 4.71E-09 0.00E+00 9.00E-08 1.44E-05 7.15E-09 1.56E-03 1.78E-04 1.33E-06

Table 2.10: P-values for White test null hypothesis of linearity in mean

Nonlinearity The bispectrum of Site 664 suggested a possible nonlinear component in the series at approximately λ = .01. Table 2.10 displays p-values for the White test null hypothesis of linearity in the mean. The test does not reject the null hypothesis of linearity for any of the untuned or tuned sample sites. A test applied to the first difference in the series, however, rejects the null hypothesis at the 99 percent confidence level for all sites.

2.2.4 Alternative Assumptions An examination of the residuals obtained by a sum-of-cosines regression using the frequencies suggested by peaks in the periodogram (as previously identified) suggests that periodic components persist in the series after removal of the signal (Figure 2.13). This may merit a reconsideration of the assumption that the residual series is gaussian white noise. Figures 2.14 and 2.15 provide analytic results with K estimated by an alternative statistic that does not rely on a detailed understanding of the residual process ([47], p 80). The results are somewhat different from those

32

−40 −80

Chiu’s R (dB)

0

Residuals, Site 664

0.0

0.1

0.2

0.3

0.4

0.5

frequency (cycles/kyr)

0 −30 −60

Chiu’s R (dB)

Residuals, Site 664 (tuned)

0.0

0.1

0.2

0.3

0.4

0.5

frequency (cycles/kyr) Figure 2.13: Residual analysis (Chiu’s R) of residuals from sum of cosines regression with bλ identified by periodogram peaks.

33

2

Time Series

Residuals from QF Method 2

0

0 −2 790 0.5

0

0

−1

κ

κ

−2 790 1

0

2

0

−0.5

4

0

2

0

2

4

0

3.4e+011

0.01 0

0

0.005 0.01 0.015 cycles/kyr

cycles/kyr

cycles/kyr

−2 790

0 4.77e+010

0.01 0

0

0.005 0.01 0.015 cycles/kyr

Figure 2.14: Results from the alternative order statistic of [47] - LEFT: Site 664 (untuned), the filtered periodogram (frequency in radians), the bispectrum. RIGHT: sum-of-cosines residuals, filtered periodogram of residuals (frequency in radians), estimated signal (green) vs. time series (blue), bispectrum of residuals.

34

Time Series

Residuals from QF Method

2

2

0

0

0

−1

−2 790 0.5

0

κ

κ

−2 790 1

0

1

2

3

0

−0.5

4

0

2

0

1

2

3

4

0

2.92e+011

0.01 0

0

0.005 0.01 0.015 cycles/kyr

cycles/kyr

cycles/kry

−2 790

0 3.21e+010

0.01 0

0

0.005 0.01 0.015 cycles/kyr

Figure 2.15: Results from the alternative order statistic of [47] - LEFT: Site 664 (tuned), the filtered periodogram (frequency in radians), the bispectrum. RIGHT: sum-of-cosines residuals, filtered periodogram of residuals (frequency in radians), estimated signal (green) vs. time series (blue), bispectrum of residuals.

35 reported above. A single periodic component is detected for both the tuned and the untuned records of Site 664. For the untuned series, the QF method estimates a single harmonic corresponding to approximately 96 kyrs; bλ = 0.01031 (-0.07403, 0.09466). Results for the tuned record suggest a harmonic of approximately 96 kyrs; bλ= 0.010456817 (0.007481957,0.013431676). As before, the confidence intervals for the frequency estimates are quite wide. Estimated signals (third panel on right in Figure 2.14) provide only a partial tracking of the features of the data, with many high frequency fluctuations and many peaks in the series not captured. A peak in the bispectrum at approximately λ = 0.01 remains for both records after the estimated signal at that frequency is removed. This suggests a possible non-linear response to the 100 kyr harmonic. Results from the QF method using an alternative estimator of K that does not depend on detailed knowledge of the residual process are provided in Tables 2.11 and 2.12 along with 95 percent confidence interval from the asymptotic distribution and bootstrap statistics from a resampling b procedure as per [25]. For each site a single cosinusoid (K=1) is estimated. Generally these values are centered around the 100 kyr cycle although the corresponding frequency estimates correspond to periodicities between 88 and 158 kyrs. Bootstrapped confidence intervals tend to be more narrow than those implied by the asymptotic distribution. The bootstrapped intervals for the tuned series are generally narrower than those for the untuned series. As an alternative to the sum-of-cosines model, an autoregressive fit to the series was also attempted. Chiu’s R statistic with critical values for the residuals of the AR fit to the data are provided in Figure 2.16. We cannot reject the null hypothesis of no signal in the noise for any frequency. AR coefficients and associated standard errors are presented in Table 2.13. Prewhitened spectrum estimates (based on the autoregressive coefficients) are presented in Figure 2.17, again

36

502 552 607 659 664 846 849 925 982

bλ 0.011 0.9984 0.0102 0.0097 0.0103 0.1045 0.1131 0.0098 0.0112

1/bλ 90.72 119.16 97.68 102.67 96.76 95.64 88.37 101.37 88.56

95% CI (λ) (-.0305,0.0526) (-.0675,.0843) (-.0304,.0508) (-.0397,.0592) (-.0154,.1747) (-.0142,.0351) (-.0165,.0392) (-.0212,.0409) (-.026,.049)

BS µ(sd) .0112 (0.0014) .0083 (0.0007) .0102 (0.00005) .0097 (0.00003) .0103 (0.00006) 0.0104 (0.0001) 0.0112 (0.0003) 0.0106 (0.0038) 0.0112 (0.0024)

Table 2.11: QF estimates of frequency (other parameter estimates provided in the supplemental appendix). 95% confidence interval obtained from the asymptotic distribution is provided as well as the mean and SD of the frequency estimate from 500 bootstrap simulations. K is estimated by the alternative statistic in [47]. Data is untuned. .

502 552 607 659 664 846 849 925 982

bλ 0.00890 0.00629 0.01010 0.00980 0.01040 0.01010 0.01000 0.00980 0.00790

1/bλ 112.36 158.98 99.01 102.04 96.15 99.01 100.00 102.04 126.58

95% CI (λ) (-.0252,.0431) (-.009,.0216) (-.007,.0272) (-.0136,.0333) (-.018,.039) (-.0674,0.0877) (-.0117,.0318) (-.0075,.0271) (-.0097,.0256)

BS µ(sd) .0097 (0.0019) .0072 (0.0016) .01010 (0.00004) .0098 (0.00003) .0105 (0.0002) 0.0101 (0.00006) 0.0100 (0.00005) 0.0098 (0.00006) 0.0069 (0.0014)

Table 2.12: QF estimates of frequency (other parameter estimates provided in the supplemental appendix). 95% confidence interval obtained from the asymptotic distribution is provided as well as the mean and SD of the frequency estimate from 500 bootstrap simulations. K is estimated by the alternative statistic in [47]. Data is tuned.

37 Untuned AR Coef 2.48710064 -2.98008757 2.64650976 -2.1558607 1.76231996 -1.40376007 1.11060585 -0.83111241 0.57745748 -0.38061594 0.21362610 -0.074562270

S.E. 0.000564137 0.002456986 0.003270882 0.003331306 0.003331301 0.003331235 0.003331235 0.003331301 0.003331306 0.003270882 0.002456986 0.000564137

Tuned AR Coef 1.8618873 -1.2821619 0.5277556 -0.1355774

S.E. 0.03137174 0.06153514 0.06153514 0.03137174

Table 2.13: AR estimate for site 664

showing the contribution of low frequency components to the variance of the series. However, [47] (p 128) show that autoregressive frequency estimates are unbiased if and only if the residual variance is zero; one must be cautious of obtaining frequency estimates from such methods.

2.3 Some Bivariate Relations Thus far we have presented some results on individual forcing mechanisms and time series based on sediment core records. The preceding discussion has suggested that the sample sites may possess an approximately 100 kyr cycle, similar to the one exhibited by inclination. Insolation’s harmonics tend to resemble those of obliquity and precession, which are generally not exhibited in the sample sites. Figures 2.18 and 2.19 examine the phase and coherence relations between the sample sites, insolation at 65 o N and inclination. The series for Site 664 were lagged so that the estimated low frequency phase delay are close to zero. The estimates do not differ considerably for the untuned

38

−20 −35 −50

Chiu’s R (dB)

AR residuals untuned

0.0

0.1

0.2

0.3

0.4

0.5

0.4

0.5

frequency (cycles/kyr)

−30 −50

Chiu’s R (dB)

AR residuals tuned

0.0

0.1

0.2

0.3

frequency (cycles/kyr)

Figure 2.16: Residuals from AR fit to series

0.01

0.50

39

0.0

0.1

0.2

0.3

0.4

0.5

0.4

0.5

0.01

0.50

lambda

0.0

0.1

0.2

0.3 lambda

Figure 2.17: Spectrum estimates from associated AR fits to series (prewhitened). Frequency (lambda) is in cycles/kyr. TOP: untuned, BOTTOM: tuned.

40 and the tuned records. After lagging the series to reduce the phase angle at the low frequencies, higher frequency signals remain out of phase. Confidence intervals for the phase (provided by [32] and [8]), however, suggest that at these frequencies we can not trust the estimates of phase. The coherence estimates do not seem to suggest a strong linear relationship between insolation and Site 664. The maximum value of the coherence between the tuned record and insolation is approximately 0.40 and the upper 95 percent confidence limit is about 0.6. This value corresponds to the period of approximately 32 thousand years; a peak near this approximate frequency was observed for some sites. Inclination shows a much stronger coherency with the time series of Site 664. The maximum value of the coherence is approximately 0.75 for the untuned record, with upper confidence limit approaching 0.9. This estimate corresponds to a periodicity of approximately 750 kyrs. The coherency close to the 100 kyr cycle is slightly weaker with a value of approximately 0.58.

2.4 Chapter Summary The results of this chapter do not appear supportive of the Milankovitch hypothesis. Periodic components detected are somewhat more similar to the low-frequency periodicity of inclination. The estimated squared coherency between Site 664 and inclination is also greater than the relationship between insolation and the site. Filtered series appear to exhibit relative stability in the approximately 100 kyr cycle. These results, however, do not rule out the possibility of a more complicated non-linear response of glacial accumulation to insolation, although support for nonlinearity in Site 664 is tenuous.

41

Insol vs. Untuned

1 −3

−1

phase

2

0.8 0.4 0.0

squared coherency

3

Insol vs. Untuned

0.1

0.2

0.3

0.4

0.5

0.0

0.1

0.2

0.3

0.4

frequency

frequency

Insol vs. Tuned

Insol vs. Tuned

0.5

1 −3

−1

phase

2

0.8 0.4 0.0

squared coherency

3

0.0

0.0

0.1

0.2

0.3

frequency

0.4

0.5

0.0

0.1

0.2

0.3

0.4

0.5

frequency

Figure 2.18: Phase and Coherence estimates for Insolation 65 o N and Site 664. Frequency is in cycles/kyr.

42

Inclin vs. Untuned

1 −3

−1

phase

2

0.8 0.4 0.0

squared coherency

3

Inclin vs. Untuned

0.1

0.2

0.3

0.4

0.5

0.0

0.1

0.2

0.3

0.4

frequency

frequency

Inclin vs. Tuned

Inclin vs. Tuned

0.5

1 −3

−1

phase

2

0.8 0.4 0.0

squared coherency

3

0.0

0.0

0.1

0.2

0.3

frequency

0.4

0.5

0.0

0.1

0.2

0.3

0.4

0.5

frequency

Figure 2.19: Phase and Coherence estimates for Insolation 65 o N and Site 664. Frequency is in cycles/kyr.

43

Chapter 3

Uncertainty in the Time-Depth Correspondence In the previous chapter the sedimentation rate assumptions, and the resulting time-depth mapping, were taken as given. In the present chapter we evaluate the uncertainty of the mapping function. The intention is to not only assess the imperfections of the mapping functions employed in this paper, but also to provide general criteria for the evaluation of the performance of depth-time mappings. Some initial discussion of the potential problems of sedimentation rate assumption imperfections is warranted. The problems of ”chatter” and ”drift” are discussed in detail in [40] (ch. 5). For our purposes it sufficient to say that chatter and drift occur when the variation in the sedimentation rate occurs at a higher or lower frequency, respectively, than the underlying climate signal. Chatter tends to move periodogram peaks to higher frequencies, while drift tends to produce some broadening of peaks. In the current analysis, we do not construct our own age model, relying rather on either a constant sedimentation assumption or on the tuning of [30]. As such, we are not

44 interested in chatter and drift problems specifically, but rather combine these under the rubric of ”mapping imperfections,” which we explore below. As we are dealing with a number of series that are assumed to be signatures of the same underlying climate process, it is also of some interest to consider series alignment. Series misalignment will be an issue for a pair of series when there are mapping imperfections in one or both series. Misalignments are caused by mapping imperfections (a special case involves missing data), but the absence of misalignment does not necessarily mean that the time-depth correspondence is reasonable; that is, misalignment implies poor mapping, but poor mapping does not imply misalignment (e.g. the mapping function may be inappropriate for both series). Alignment problems themselves are not problematic for the analysis of the previous chapter because we examined the harmonic properties of the series individually and did not present results that were based on a comparison of two individual sample sites. However, for stacked (i.e. composite) records, such as SPECMAP or the Benthic Stack that was constructed from the data analyzed in the previous chapter, alignment is an issue. Asynchronous phase averaging of misaligned records can contribute to a dampening of the composite signal at a particular frequency. For example, if two series exhibit a strong 23 kyr periodic component but are out of phase, a simple averaging of the series may eliminate this harmonic from the composite series. Techniques for the alignment of multiple time series records are developed in [57] and [1]. We should note that a weak global signal may look like series misalignment along the criteria measured below.

45

3.1 Bivariate Relations The simple correlation estimates between the series (complete tables provided in the supplemental appendix) are less than 0.6 for each pair of series, with some pairs exhibiting a negative correlation. For example, Site 664 and Site 982 have an estimated correlation of -0.13. The estimated sample correlations generally improve for the tuned records, to as high as 0.71, but remain low (as low as 0.09) for many of the pairs. Figure 3.1 provides phase and coherency estimates between Site 664 and other sites. The squared coherency of the 100 kyr signal of sites 846 and 664 is approximately 0.8, but the series do not exhibit a strong coherence at other frequencies. The figure also presents the coherency between the signals at Sites 664 and 982. The coherency is much weaker for this site combination. Phase angles (right side) suggest that the signals may not be in-sync. Interestingly, coherency and phase plots for the tuned series (not shown) generally decrease the coherence and phase synchronization between the site pairs. These results suggest alignment problems and, in turn, the series may suffer from mapping imperfections. Alternatively, the results could also be an indication of a weak global signal.

3.2 Complex Demodulation Complex demodulation of the data for Site 664 raises additional questions. In Figure 3.2, the power decreases somewhat in certain ranges of the data (e.g. 400-550 kyrs before present). If the mapping function is determined correctly, an astronomically determined signal should exhibit a relatively stable relationship over time.

1 −3

−1

phase

2

0.8 0.4 0.0

squared coherency

3

46

0.0 0.1 0.2 0.3 0.4 0.5

frequency

frequency

1 −3

−1

phase

2

0.8 0.4 0.0

squared coherency

3

0.0 0.1 0.2 0.3 0.4 0.5

0.0 0.1 0.2 0.3 0.4 0.5

0.0 0.1 0.2 0.3 0.4 0.5

frequency

frequency

Figure 3.1: TOP ROW: coherency and phase estimates Sites 664 and 846; BOTTOM ROW: same figures Sites 664 and 982. These data are untuned and frequency is in cycles/kyr.

47

−7 −4 −1

untuned

0

200

400

600

800

600

800

kiloyears before present

−7 −4 −1

tuned

0

200

400

kiloyears before present

Figure 3.2: Complex demodulation (log amplitude) at λ=0.01, Site 664

48

3.3 Dynamic Frequency Analysis Our examination of the bivariate relationships among the samples in the previous sections, suggested that the mapping produced misalignment of the series. The identification of periodogram peaks adjacent to widely accepted cycles of interest (as per inclination, eccentricity, precession, obliquity and insolation frequencies) in several of the sites may be an indication of chatter and drift problems associated with the mapping function. In ([40], 157), the authors note that the posited astronomical origins of ice ages imply that the harmonic properties obtained from core data should not vary over time for the past million years when the mapping correspondence is correctly specified. Thus, if we consider subsets of the time series data, there should be little observed difference in the frequency properties of different time slices. A high resolution trace of the evolution of the dominant harmonic signal over the sample period is provided in Figure 3.3. The figure plots the period associated with the periodogram maximizer over time (a discussion of dynamic frequency analysis is provided by ([56], 5.10)). A window width was fixed and the periodogram was computed on a rolling (i.e. overlapping) window of observations. We focus on the periodogram maximizer because of the favorable properties of that estimate of frequency presented in [47]. A plot of the evolution of the periodogram maximizer suggests some imperfections in the mapping function. The periodogram maximizer does not appear to be stable over the sample period for either the untuned or the tuned series although changes are infrequent and short. These results are sensitive to the size of the data window: doubling the size of the window stabilizes the periodogram maximizer of the untuned record at the 100 kyr cycle, but the results for the tuned record remain similar to those reported in Figure 3.3.

49

105 90 75

kyrs/cycle

untuned

0

100

200

300

400

500

600

400

500

600

index

105 90 75

kyrs/cycle

tuned

0

100

200

300 index

Figure 3.3: Evolution of the periodogram maximizer (1/bλ).

50

Chapter 4

Concluding Remarks and New Directions This paper has examined the harmonic properties of nine δ18 O deep sea core records as proxies for the evolution of climate during the past 790 thousand years. We examined regularly sampled time series constructed by assuming a constant rate of sedimentation between the identified Brunhes-Matuyama geomagnetic reversal tie point and the present, and also the tuned series constructed by [30]. Our intention has been to apply simple stochastic models and widely accepted time series techniques to approximate some general frequency properties of the series and their match to hypothesized forcing data. Results do not appear to support the Milankovitch hypothesis of summer insolation forcing of climate. These results are somewhat tenuous, however, because of some imperfections in the assumed time-depth correspondence for both the tuned and the untuned series. Possible misalignments among the series are a potential indication of a poor mapping correspondence. Also, an examination of the evolution of the harmonic properties generated by rolling (overlapping) data windows pro-

51 vides some indication that the harmonic properties of the series do not exhibit the stable relationship over time that might be expected when the process is driven by astronomical forces. These results may be an indication that greater refinements in the tuning process are required before a stable relationship between astronomical parameters and climate can be determined. Yet, we must note that according to several measures, the assumed sedimentation rates appear to be broadly reasonable and that series misalignments may be of greater relevance to the construction of stacked series than to the method of investigation employed in the current paper. Therefore the general inconsistency with the Milankovitch hypothesis can not be discounted entirely, although we do not rule out the possibility that more complicated non-linear dynamics relate insolation to glacial accumulation. We reserve the final comments of this paper for a brief discussion of some anticipated future research.

Point Process Analysis In the present paper the marked point processes of δ18 0 records were converted to a regularly sampled time series for some added convenience in subsequent analysis. This transformation of the series is not essential, and indeed many interesting possibilities arise from keeping the data in the form of a point-process, with the notation as described at the beginning of Chapter 2. Pertinent techniques for frequency domain analysis of point-process data are described in [15]; a discussion of the fit of sum-of-cosines models to point-process data is provided by [14]. An interesting possibility for a parametric tuning is provided by the hybrid process notation. One could imagine tuning the time scale of the data by minimizing the difference in some function of the series, subject to time scale parameter restrictions. We might begin by specifying a discrete time-depth function such as τ(l j ) = al j + by(l j ) + ρ cos(ωl j + ψ), with l j representing the

52 cumulative (total) depth below the sea floor, and the second and third terms on the right hand side are some values that might correlate with deviations from the average sedimentation rate. The periodic component may capture cycles of biogenous productivity, while the inclusion of the sampled series may capture the relationship that is likely to exist between glacial cycles and sea floor accumulation in the traditional spirit of the tuning problem. The basic concept would be to minimize the difference in some function of the δ18 O and orbital records subject to the restrictions of a mapping correspondence such as the one above. The major difficulty implied by the above method is that it would seem any change in θs = {a, b, ψ, ρ, ω} implies a change in the data, especially in the case where one seeks to use an interpolation of the continuous process. That is, if working with a regularly sampled time series, a change in a parameter estimate implies a change in the depth-time correspondence and requires that one reconstruct the time series by interpolation along the new time scale prior to the next iteration in the estimation method. Also, tuning adjustments to the time scale make specifying a time series relationship (i.e. a function to minimize) between the orbital series O(t) and the δ18 O, y(t), difficult because if one time scale is changed, the relation between y(τ(li )) and O( t j ) may not be as relevant as the relation between y(τ(lk )) and O( t j ). These problems are especially difficult if one deals with a regular time series constructed from the depth records because much of the data used may be interpolated, not directly observed data. Consider a frequency domain specification of the minimization function, applied directly to the hybrid data process. For example, let dby = ∑Lj=1 y(τ(l j ))e−iτ(l j )λ be the DFT of the (perhaps normalized in some way) marked point process of δ18 O values, with some similar specification for the orbital series (see [15]). Then consider an adjustment to the time scale by a change in θs .

53 The adjustment changes the depth-time correspondence in a clearly specified functional way and, most importantly, does not require a recomputation of the marked point process data because the frequency domain transform is applied directly to the observed sample data, not the interpolated. Comparisons between the sample sites and orbital records take place in the frequency domain so only the DFT magnitudes change with no adjustment in the relative position of the data as mentioned above.

Threshold Modelling The analysis of Chapter 2 suggests some possibility of a non-linear component in the δ18 O. In addition, the Site 664 data plotted in Figure 2.5 exhibits periods of rapid warming at the terminus of a long period of slow cooling. In the paleoclimate literature this is referred to as a glacial termination (see [40]). These properties suggest that some form of threshold model may merit investigation. An informal threshold model was developed by [52], but was criticized in [40]. No threshold value is obvious from an examination of the series, so one might consider a model with such value endogenously determined within the stochastic model.

Joint Determination The joint determination of harmonics and sedimentation rates does not appear to have received much attention in the literature. By joint determination, we mean devising sedimentation rates without an astronomical template while making inferences regarding the consistency of the resulting series with astronomical forcing mechanisms. Instead of joint determination, most work devises a sedimentation scale and then evaluates the harmonics of the resulting series. The work of [3] is one example that makes use of the assumption that spectrum peaks in the resulting series

54 should be narrow in devising sedimentation rates and so does not confine the data to a predetermined astronomical origin such as insolation. The joint determination problem may involve some functional specifications of the interplay between astronomical forces and climate. One could also provide a parametric specification of the sedimentation process similar to that discussed above. One difficulty, again, is that any timedomain adjustment to the sedimentation function requires a reinterpolation of the entire discretevalued δ18 O series, so that a match to various astronomical function specifications is more difficult because the entire data series may completely change with any adjustment. Yet, a seemingly analogous problem exists in the econometric literature. [20] derive a multi-equation business cycle model of the macroeconomy. They use the Generalized Method of Moments to estimate model parameters with some twenty orthogonality conditions that may be derived from their model. The interesting and relevant aspect of their approach is that they have a time series of capital depreciation that exhibits similar problems to our unknown sedimentation rates. The depreciation time series is assumed to be measured correctly, on average, but quarterly deviations from the average are unobserved; depreciation is fundamental to determining the total stock of capital, which is a key variable in the overall macroeconomic model. Their model uses the implications of the orthogonality conditions to approximate parameters for depreciation, which are then used to construct a depreciation time series and subsequently reconstruct the entire capital stock series with each iteration of the estimator.

The forgoing discussion has not delved into some specific difficulties encountered in the implementation of the techniques. We are currently at work resolving some of the practical limitations to implementing these methods and look forward to their publication in future research.

55

Bibliography [1] Bjorn Auestad, Robert H. Shumway, Dag Tjostheim, and Kenneth Verosub. Nonlinear alignment of time series with applications to varve chronologies. Unpublished manuscript, University of California, Davis, courtesy Robert Shumway, 2001. [2] A. Banerjee, J. J. Dolado, J. W. Galbraith, and D. F. Hendry. Cointegration, Error Correction, and the Econometric Analysis of Non-Stationary Data. Oxford University Press, Oxford, 1993. [3] John Belcher and Granville Tunnicliffe Wilson. Time scale estimation by tracking parameter variation. Journal of Time Series Analysis, 21(3):237–248, 2000. [4] Y. Benjamini and Y. Hochberg. Controlling the false discovery rate: A practical and powerful approach to multiple testing. Journal of the Royal Statistical Society, Series B, 50:289–300, 1995. [5] A. Berger and M.F. Loutre. Insolation values for the climate of the last 10 million years. Quaternary Sciences Review, 10(4):297–317, 1991. [6] W.H. Berger, L.W. Kroenke, L.A. Mayer, and et al., editors. Proceedings of the Ocean Drilling Program, Scientific Results, Volume 130. Ocean Drilling Program, College Station, TX, 1993. [7] T. Bickert, W.H. Berger, S. Burke, H. Schmidt, and G. Wefer. Late quaternary stable isotope

56 record of benthic foraminifers: Sites 805 and 806, ontong java plateau. In Proceedings of the Ocean Drilling Program, Scientific Results, 130:411–420, 1993. [8] Peter Bloomfield. Fourier Analysis of Time Series. Wiley, New York, 2000. [9] Edward A. Boyle. Characteristics of the deep ocean carbon system during the past 150,000 years: σco2 distributions, deep water flow patterns, and abrupt climate change. Proceedings of the National Academy of Sciences of the United States of America, 94(16):8300–8307, 1997. [10] David R. Brillinger. The spectral analysis of stationary interval functions. In L.M. Le Cam and J. Neyman and E.L. Scott (eds.), Proceedings of the Sixth Berkeley Symposium on Mathematical Statistics and Probability, pages 483–513, 1972. [11] David R. Brillinger. The analysis of time series collected in an experimental design. Multivariate Analysis - III, pages 241–256, 1973. [12] David R. Brillinger. Time Series: Data Analysis and Theory. Holden-Day, San Francisco, 1981. [13] David R. Brillinger. In the finite fourier transform of a stationary process. In David R. Brillinger and P.R. Krishnaiah (eds.), Handbook of Statistics, Vol. 3, pages 21–37, 1983. [14] David R. Brillinger. Fitting cosines: Some procedures and some physical examples. In Applied Probability, Stochastic Processes, and Sampling Theory, pages 75–100, 1987. [15] David R. Brillinger. Time series, point processes, and hybrids. The Canadian Journal of Statistics, 22(2):177–206, 1994.

57 [16] David R. Brillinger. Examining an irregularly sampled time series for whiteness. Resehas IME-USP, 4(4):423–431, 2000. [17] David R. Brillinger. Point Process Analysis Lecture Notes. University of California, Berkeley, 2001. [18] David R. Brillinger and P.R. Krishnaiah, editors. Handbook of Statistics, Vol. 3. Elsevier, Amsterdam, 1983. [19] Peter J. Brockwell and Richard A. Davis. Time Series: Theory and Methods. Springer-Verlag, New York, 1991. [20] Craig Burnside and Martin Eichenbaum. Factor-hoarding and the propagation of businesscycle shocks. The American Economic Review, 86(5):1154–1174, Dec 1996. [21] L.M. Le Cam, J. Neyman, and E.L. Scott, editors. Proceedings of the Sixth Berkeley Symposium on Mathematical Statistics and Probability. University of California Press, Berkeley, CA, 1972. [22] M. Casdagli and S. Eubank, editors. Nonlinear Modeling and Forecasting. Addison-Wesley, New York, 1992. [23] Shean-Tsong Chiu. Detecting periodic components in a white gaussian time series. Journal of the Royal Statistical Society, Series B (Methodological), 51(2):249–259, 1989. [24] Lawrence Christiano and Terry Fitzgerald. The band pass filter. Unpublished manuscript, Northwestern University, 1999.

58 [25] David A. Freedman and Stephen C. Peters. Bootstrapping a regression equation: Some empirical results. Journal of the American Statistical Association, 79(385):97–106, 1984. [26] T. Hasan. Complex demodulation: Some theory and applications. In D. R. Brillinger and P.R. Krishnaiah (eds.), Handbook of Statistics, Vol. 3, pages 125–156, 1983. [27] J. Hays, J. Imbrie, and N. Shackleton. Variations in the earth’s orbit: Pacemaker of the ice ages. Science, 194:1121–1132, 1976. [28] J. Imbrie, J. Hays, D. Martinson, A. McIntyre, A. Mix, J. Morley, N. Pisias, W. Prell, and N Shackleton. The orbital theory of pleistocene climate: Support from a revised chronology of the marine δ18 o record. In A. Berger et. al. (eds.), Milankovitch and Climate, Part I, pages 269–305, 1984. [29] G.R. Johnson. Brunhes-matuyama magnetic reversal dated at 790,000 yr b.p. by marineastronomical correlations. Quaternary Research, 17:135–147, 1982. [30] Daniel Karner, Jonathan Levine, Brian Medeiros, and Richard Muller. Constructing a stacked benthic δ18 o record. Paleoceanography, 17:2002, 3. [31] Stephen M. Kay. Modern Spectral Estimation: Theory and Application. Prentice Hall, Englewood Cliffs, NJ, 1988. [32] L.H. Koopmans. The Spectral Analysis of Time Series. Academic Press, New York, 1974. [33] P.R. Krishnaiah, editor. Multivariate Analysis - III. Academic, New York, 1973. [34] T.H. Lee, H. White, and C.W.J. Granger. Testing for neglected nonlinearity in time series models. Journal of Econometrics, 56:269–290, 1993.

59 [35] Richard A. Muller. Glacial cycles and orbital inclination. Lawrence Berkeley Laboratory Report, LBL-35665, 1994. [36] Richard A. Muller and Gordon J. MacDonald. Glacial cycles and orbital inclination. Nature, 377:107–8, 1995. [37] Richard A. Muller and Gordon J. MacDonald. Glacial cycles and astronomical forcing. Science, 277:215–218, 1997. [38] Richard A. Muller and Gordon J. MacDonald. Simultaneous presence of orbital inclination and eccentricity in proxy climate records from ocean drilling program site 806. Geology, 25:3–6, 1997. [39] Richard A. Muller and Gordon J. MacDonald. Spectrum of 100-kyr glacial cycle: Orbital inclination, not eccentricity. Proceedings of the National Academy of Sciences of the United States of America, 94(16):8329–8334, 1997. [40] Richard A. Muller and Gordon J. MacDonald. Ice Ages and Astronomical Causes. Springer, New York, 2000. [41] H. Joseph Newton and Marcello Pagano. A method for determining periods in time series. Journal of the American Statistical Association, 78(381):152–157, 1983. [42] Efstathios Paparoditis and Dimitris N. Politis. The local bootstrap for periodogram statistics. Journal of Time Series Analysis, 20(2):193–222, 1999. [43] P. Perron. Trends and random walks in macroeconomic time series. Journal of Economic Dynamics and Control, 12:297–332, 1988.

60 [44] M. B. Priestley. Detection of periodicities. In T. Subba-Rao, M.B. Priestley and O. Lessi (eds.), Applications of Time Series Analysis in Astronomy and Meteorology, pages 65–88, 1997. [45] M.B. Priestley. Spectral Analysis and Time Series. Academic Press, New York, 1981. [46] Barry G. Quinn. Estimating the number of terms in a sinusoidal regression. Journal of Time Series Analysis, 10:71–75, 1989. [47] Barry G. Quinn and Edward J. Hannan. The Estimation and Tracking of Frequency. Cambridge University Press, Cambridge, U.K., 2001. [48] T. Subba Rao. Analysis of nonlinear time series (and chaos) by bispectral methods. In M. Casdagli and S. Eubank (eds.) Nonlinear Modeling and Forecasting, pages 199–226, 1992. [49] T. Subba Rao, editor. Developments in Time Series Analysis. Chapman and Hall, New York, 1993. [50] T. Subba Rao, M.B. Priestley, and O. Lessi, editors. Applications of Time Series Analysis in Astronomy and Meteorology. New York, 1997. [51] M.E. Raymo. The timing of major climate terminations. Paleoceanography, 12(4):557–585, 1997. [52] M.E. Raymo. Glacial puzzles. Science, 281:1467–1468, 1998. [53] W.F. Ruddiman, M.E. Raymo, D.G. Martinson, B.M. Clement, and J. Blackman. Pleistocene evolution: Northern hemisphere ice sheets and north atlantic ocean. Paleoceanography, 4:353–412, 1989.

61 [54] J.J. Shackleton, A.L. Berger, and W.R. Peltier. An alternative astronomical calibration of the lower pleistocene timescale based on odp site 677. Transactions of the Royal Society of Edinburgh, Earth Sciences, 81:251–261, 1990. [55] Chad Shafer. Point Process Analysis Labs. University of California, Berkeley, 2001. [56] Robert H. Shumway and David S. Stoffer. Time Series Analysis and Its Applications. Springer, New York, 2000. [57] Robert H. Shumway and Kenneth L. Verosub. State space modeling of paleoclimatic time series. Proceedings of the 5th International Meeting on Statistical Climatology, June:22–26, 1992. [58] James Stock. Estimating continuous-time processes subject to time deformation: An application to postwar u.s. gnp. Journal of the American Statistical Association, 83:77–85, 1988. [59] T. Teraesvirta, C. F. Lin, and C. W. J. Granger. Power of the neural network linearity test. Journal of Time Series Analysis, 14:209–220, 1993. [60] W.N. Venables and B.D. Ripley. Modern Applied Statistics with S-Plus. Springer, New York, 1999. [61] Timothy J. Vogelsang. Two simple procedures for testing for a unit root when there are additive outliers. Journal of Time Series Analysis, 20(2):237–252, 1999.

62

Appendix A

63

Site DSDP 502 DSDP 522 DSDP 607 ODP 659 ODP 664 ODP 677 ODP 846 ODP 849 ODP 925 ODP 982 RC13-110 RC13-229 ODP 805 ODP 806

Location 11 ◦ 30’N, 79 ◦ 23’W 56 ◦ 03’N, 23 ◦ 14’W 41 ◦ 00’N, 32 ◦ 58’W 18 ◦ 05’N, 21 ◦ 02’W 00 ◦ 06’N, 23 ◦ 14’W 01 ◦ 12’N, 83 ◦ 44’W 03 ◦ 06’S, 90 ◦ 49’W 00 ◦ 11’N, 110 ◦ 31’W 04 ◦ 12’N, 43 ◦ 29’W 57 ◦ 31’N, 15 ◦ 53’W 0 ◦ 6’N, 95 ◦ 39’W 26 ◦ S, 11 ◦ E ◦ 1 14’N, 160 ◦ 33’E 0 ◦ 19’N, 159 ◦ 22’E

Present

GT II

GT V

BM Bound

m Rate ( kyr )

0 0 0 0 0 0 0 0 0 0 0 0 0 0

2.2 2.4 5.8 4.2 4.6 6.1 5.4 4.0 5.1 2.9 3.0 3.2 1.9 2.6

7.8 6.9 17.4 13.8 15.7 17.1 16.3 13.1 16.7 9.9 8.8 9.1 6.3 8.8

16.2 14.2 31.7 24.8 28.8 30.8 28.7 22.5 28.9 20.0 NA NA NA NA

0.0205 0.0179 0.0401 0.0313 0.0364 0.0389 0.0363 0.0284 0.0365 0.0253 0.0202 0.0209 0.0145 0.0202

Table A.1: Location and Approximate Depth (mcd) of Tie Points at Sites (NA indicates the depth record does not extend to the tie point). The last column is the average sedimentation rate (m/kyr) between the present and the last available tie point. This data provided by Jonathan Levine from the [30] paper.