Time-series analysis of the Pleistocene deep-sea paleoclimatic record

Marine Geology, 36 (1980) 1--22 © Elsevier Scientific Publishing Company, Amsterdam -- Printed in The Netherlands

1

TIME-SERIES ANALYSIS OF THE PLEISTOCENE DEEP-SEA PALEOCLIMATIC RECORD

MADELEINE BRISKIN and JAMES HARRELL

Department of Geology, University of Cincinnati, Cincinnati, Ohio 45221 (U.S.A.) (Received August 30, 1978; revised and accepted July 12, 1979)

ABSTRACT Briskin, M. and Harrell, J., 1980. Time-series analysis of the Pleistocene deep-sea paleoclimatic record. Mar. Geol., 36: 1--22. Pleistocene deep-sea cores V16-205, V28-239 and composite core RC11-120/E49-18 were selected for a time-series analysis of climatic periodicities which correspond to secular changes in the earth's orbital parameters. These data were analyzed using the Periodic Regression with Cyclic Descent -- a technique which offers three major advantages over more conventional approaches because (1) the observation interval may be unequally spaced with respect to time, (2) resolution of all the discrete periodic components present in the data, including those with low frequencies, is made possible, and (3) the associated Monte Carlo experiments permit the evaluation of the importance of the observed periodic components vis-a-vis "red" climatic noise. Analytical results indicate periodicities which correspond to secular changes in the earth's orbital parameters. In particular, the results strongly suggest a quasi-periodicity of 413,000 years recorded in deep-sea globigerina oozes. Clearly this periodicity is reminiscent of a beat in the eccentricity of the earth's orbit, first detected in tropical North Atlantic core V16-205 by Briskin and Berggren (1975). INTRODUCTION In r e c e n t years m a n y geologists and climatologists have t u r n e d their attent i o n t o w a r d the o c e a n in an a t t e m p t to unravel the p a t t e r n o f Pleistocene and H o l o c e n e global climatic variations. Analysis o f the climatic o v e r p r i n t i n g in deep-sea s e d i m e n t o l o g i c r e c o r d s m a y provide us with the m u c h s o u g h t - a f t e r c o m p r e h e n s i v e secular climatic m o d e l . Briskin and Berggren ( 1 9 7 5 , fig.10) n o t e d a s t r o n g q u a s i - p e r i o d i c i t y o f 4 3 0 , 0 0 0 to 4 5 0 , 0 0 0 years in the estimated w i n t e r t e m p e r a t u r e m i n i m a o f deep-sea core V 1 6 - 2 0 5 . F u r t h e r m o r e , Briskin ( u n p u b l i s h e d d a t a ) observed a s t r o n g c o r r e s p o n d e n c e b e t w e e n these m i n i m a and the t h e o r e t i c a l e c c e n t r i c i t y m a x i m a o f the e a r t h ' s orbit. The observed p e r i o d i c i t y was clearly r e m i n i s c e n t o f the d o m i n a n t 4 1 3 , 0 0 0 y e a r beat in the orbital e c c e n t r i c i t y . Classical F o u r i e r analysis p r o v e d i n a d e q u a t e at the time and a r e s o l u t i o n o f these lowf r e q u e n c y events had to await t h e a p p l i c a t i o n o f a d i f f e r e n t a p p r o a c h to time-series analysis - - Periodic Regression with Cyclic Descent.

The main purpose of this study was to analyze several deep-sea cores for climatological periodicities and demonstrate their possible correlation to the secular changes in the earth's orbital parameters as recalculated by Andr~ Berger {1976}. This was accomplished by the application of the periodic regression with cyclic descent technique to three deep-sea cores: composite core R C l l - 1 2 0 / E 4 9 - 1 8 (Hays et al., 1976), V28-239 (Shackleton and Opdyke, 1976) and V16-205 (Briskin and Berggren, 1975; Van Donk, 1976). Because the periodic regression and cyclic descent approach used in this study is a different application of time-series analysis, it is necessary to describe the technique in detail in order that it may be clearly available to other workers interested in its application. Periodic Regression -- Background

It is often of interest to determine the periodic or cyclic repetitions present in a "serial" data set. Such a data set is referred to as a time-series, signal or record. It is composed of observations made on a dependent variable (e.g., sea surface temperature, oxygen isotope ratio, etc.) at specific points along a time line, where time constitutes the independent variable. Thus for a core, we have a suite of samples and for each sample we have paired observations on the dependent variable and down-core depth (and hence, time). One approach to time-series problems is by way of regression analysis. In this sense, the analysis involves fitting, in a "least squares" manner, sinusoidal curves to the data. The wavelengths, amplitudes and phases (see Fig.l) of those "best fit" sinusoidal waveforms represent, respectively, the frequencies, the relative contributions to the total variance, and the locations of the maxima and minima of the "periodic c o m p o n e n t s " in a time-series. Of course, it is important to note that the significance of least-squares fits to time-series data rests not only on significant variance reduction, b u t also on the absence of significant autocorrelations in the residuals formed after subtraction of the best-fit sinusoidal waveforms from the d a t a . It is often desired (1) to analyze serial data having observations which are unequally spaced with respect to time, and (2) to be able to resolve all periodic components actually present in the data, especially those with the longer wavelengths. The standard approaches to time-series analysis -- classical Fourier analysis and the method of the "truncated autocovariance function" (Blackman and Tukey, 1958) -- are incapable of achieving these goals. These methods were designed to treat series of equally spaced observations and, without additional manipulations ( cf., Walker, 1971), will resolve only f~equency intervals or bands. These bands are delimited b y frequencies which are harmonically related to the length of the time-series in the case of the first approach and to the maximum lag of the autocovariance function in the case of the second. Furthermore, these approaches lack resolution in the low to intermediate frequency portion of the power spectrum and require seemingly endless "fiddling" with the data and results -- filtering, window carpentry, window closing, etc.

,

!

g Y~

L

oi t~ -

T I M E

2~r T -

-

Fig.1. A t t r i b u t e s o f a s i n u s o i d a l curve. See t h e d i s c u s s i o n o f eq.1 f o r d e f i n i t i o n s o f t h e symbols.

In this paper we discuss an alternative approach to time-series analysis -Periodic Regression. Periodic regression allows us to achieve the two goals enumerated above -- obtaining estimates of any discrete waveforms present in a time-series where the observation interval may be nonuniform. The variances associated with each of these waveforms will sum to the original variance of the time-series; and the waveforms, when combined, will be practically indistinguishable from the Fourier series representation o f a continuous record. Thus any time-series may be decomposed into a set of sinusoidal waveforms which (1) have arbitrary integer frequencies which are harmonics of the length of the record (classical Fourier analysis) or (2) have noninteger frequencies and which correspond to the actual periodic components which might be present in the data (periodic regression). The results of periodic regression are depicted in a line spectrum reminiscent of the Schuster "periodogram". Other approaches to time-series analysis which are similar to periodic regression have been described by Schickendanz and Bowen (1977) and James (1966) PERIODIC REGRESSION

Development

of the approach

Consider the sinusoidal curve of Fig.1. The equation for this curve is: (f2~ t i ) Yi = Y - - A S + A c ° s ( f Oi - P ) = ~ - A S + A c o s \ ~ --P

where 0 i = the ith value of time, in radians (i = 1, 2 . . . . t i = the ith value of time, in units of time; 0 i a n d t i are related by the identity: 0i

ti -

2~

T

27rti , or0i

=

-

T

, N);

(1)

where T = the total length of the abscissa or record in units of time and N is the number of paired observations (samples) Yi = the value of y, the dependent variable, corresponding to the ith value of 0 or t; = the mean of the y's; f = frequency, i.e., the number of times a waveform repeats itself in the interval T (i.e., frequency = the number of wavelengths in T, f = T/W; where W is the wavelength or period of the sinusoid). Frequency is a dimensionless quantity which may assume noninteger values. A = the amplitude of the sinusoid, in units of y; = the average value of the cosine term; P = the phase (i.e., lateral offset) of the sinusoid; on the left side of eq.1 phase is expressed in radians and on the right side in units of time. We develop the derivation further by considering the following trigonometric identity:

k cos ( u - v )

=k

cos(v) cos(u) + k sin(v) sin(u)

Thus eq.1 may be rewritten as:

Yi = ~ - - A ~ + A c o s ( P )

cos\

~[2"ti i T :

ff2"ti i

+ A sin (P) sin \ ~ !

Letting ~ = A cos(P),/3 = A sin(P) and I = :~ -- A ~, eq. 1 then becomes: Yi

=

I + c~ cos

+ ~ sin

(9.)

In time-series analysis, the general strategy is to evaluate a large number of sinusoids of prespecified frequency. It is hoped that among this number will be those frequencies actually present in the data. These will be readily identifiable from their relatively large amplitudes. We want to obtain a " b e s t " fit to the data for each of these sinusoids by selecting the optimal values for the amplitude A and phase P. The fitting and selection process can be achieved through multiple regression analysis. From eq.2 we can now define the required regression model: Yi = ~ + o~] cos

+ ~j. sin

+ e

(3)

where c~j and ~j. are the regression coefficients, and/] the intercept term, corresponding to the ]th frequency, ~.; Yi is the y value for the ith time, ti; and e is the residual (random error). In applying eq.3 it is desirable that the time-series be " s t a t i o n a r y " , that is, its statistical properties are independent of time. In order to effect approximate stationarity in the mean, it is usually considered sufficient to remove any linear trend present in the data. The det~ending is aceomplished by fitting

a linear regression line to the original data and using the residuals (Yi -- Yi, where y~ is the predicted or estimated y value) from this line in all subsequent calculations. This is the practice followed in this study. In the discussion which follows, the Yi'S will correspond to the residuals from the linear trend rather than, as previously, the original values of the dependent variable. For the residuals, the mean (Yr) is approximately zero, and the variance (s~) is equal to that of the original observations. In general, :~r will not be equal to zero because of the nonuniform observation interval and the, perhaps inappropriate, assumption of a linear trend. Formulas for the unknown regression coefficients a and t3, and intercept I are derived through application of the "least squares" criteria. Accordingly, we seek to minimize the sum of the squared residuals: N ~, ( Y i - - Y i ) ~ i=1

Letting: cij = cos ( ~ )

and sij = sin ( ~ )

Then defining the function F to be: N N F = Z (yi - yi) = Z [yi - (Ij + o~j cij + [3j sij)] 2 i=1 i=l

we require that: aF

aF -

aF -

-

o

Taking the partial derivatives we obtain the following system of simultaneous linear equations (the subscripts have been dropped to simplify the presentation, and the summations are for i = 1 , . . . , N): aF aI - --F'(Y) + N I + a Y~(c) + [3 Y,(s) = 0 aF as

--F,(yc) + I Y_.(c) + a Y~(c 2 ) + [3 ~ ( s c ) = 0

aF aft

--Y~(ys) + I F,(s) + ~ F,(sc) + [3 Y~(s 2 ) = 0

Rearranging and expressing the equations in matrix notation, we arrive at the solution for the unknowns, --

(4)

LX(s) ~(sc) ~(s ~)j

LX(ys)J

The portion of the total variance, s~, of Yi accounted for b y the sinusoidal regression curve with frequency f1 and corresponding period or wavelength T/fi is equal to the Coefficient of Determination, R ] , the square of the multiple Pearson p r o d u c t - - m o m e n t correlation coefficient: N

i=1 R] - N

~_ 2

(yi __ 4) 5

Sy

2

Sy

i=1 where all symbols are as previously defined. The quantity, R ] , is known in the terminology of spectral analysis as the " p o w e r at frequency/)". It may be thought of simply as the ratio of the explained variation to the total variation. If/i happens to be an integer and the observation interval is uniform, R] is identical to the power of t h e / j t h harmonic in classical Fourier analysis. The "phase" is the distance, in time, from the beginning of the timeseries to the nearest positive time crest (i.e., point of maximum amplitude) of the sinusoidal regression curve (see Fig.l). We determine the phase by scanning the regression eurve (eq.3) until we find the t i corresponding to the first maxima of y~. The periodic regression approach to time-series analysis may be outlined as follows: (1) The original time-series is first detrended, in order to achieve stationarity, via linear regression. In the case of a non-linear trend, detrending with higher order polynomials may be necessary. The residuals from the regression line will constitute the "residual time-series" which will be used in all subsequent calculations. The mean of the residual time-series is approximately zero. (2) The residual time-series is decomposed into a set of nonharmonie waveforms which correspond to the periodic components actually present in the data. The decomposition method employed is that of "Cyclic Descent" (Bloomfield, 1976, pp.20--25). Cyclic descent is a strategy rather than a specific algorithm. The general idea is to resolve each waveform independent of all others present in the time-series. One does this by approximately removing all periodic components from the time-series other than the one being resolved. Waveforms which have similar frequencies tend to be mutually interfering with the consequence that the period, power and phase of each is inaccurately determined (Bloomfield, 1976, p.24). For example, if these three parameters were determined for a range of wavelengths, each wavelength differing from those adjacent by a small and constant amount, serious resolution problems are to be expected. This was the procedure followed by Schickendanz and Bowen (1977). It is then this interference and, hence, error that we wish to minimize, and cyclic descent is the best means of doing this. There are many different ways in which the cyclic descent strategy may be implemented. Bloomfield (1976, p.23) presents an algorithm and F O R T R A N program for one version of cyclic descent. Our algorithm is different. The procedural steps are as follows:

(a) A range of wavelengths is scanned at some small, constant increment (i.e., calculations are made for each wavelength using eqs.4 and 5 and the procedure for phase determination). The shortest wavelength should correspond to the Nyquist frequency and be taken, for unequally spaced data, at twice the average observation interval. The longest wavelength is arbitrary and limited only b y practical considerations. The single wavelength having the greatest power is identified from the results of the scan. The corresponding waveform is subtracted from the residual time-series. Using this "modified" residual time-series, the same range of wavelengths is scanned, and the single most powerful waveform is again determined. This is subtracted from the time-series. Up to this point then, we have subtracted from the original residual time-series the two most powerful waveforms. We continue the process of resolving waveforms, one at a time, and subtracting them successively from whatever remains of the time-series. This process is stopped upon reaching some prespecified limit. This first step of the cyclic descent provides us with approximate values for the wavelengths, amplitudes, and phases of all periodic components in the data. (b) Next, the initial estimates of the q waveforms resolved in the previous step are further refined. To do this, each waveform is considered individually in the following manner. If we are dealing with the jth waveform ~] with wavelength Wj, all q -- 1 other waveforms are subtracted from the original residual time-series. A range of wavelengths, (Wj + ~ Wj) to (Wi -- ¼~), is scanned at a small, constant increment. The single most powerful waveform resolved is our revised estimate of hi. The same procedure is followed for all q -- 1 other waveforms. (c) The procedure of step (b) could be iteratively repeated until some prespecified convergence criterion on the parameters is satisfied. We chose instead to stop after repeating step (b) a second time, accepting these estimates for wavelengths, powers and phases of all periodic components in the data. (3) The results of periodic regression are graphically displayed in a line (power) spectrum where the normalized power R :n (i.e., R:nj = (R] /Rmax) 2 • 100 and j = 1 , . . . , q) is plotted against the correspondin~ frequency (expressed as the number of cycles per 10,000 years; f = 10,000/W) or wavelength. The graphics associated with this spectrum will also depict, in some cases, an evaluation of the possible significance of the observed powers (e.g., see Fig. 4). A listing of the F O R T R A N IV computer program MADI is available from the authors upon request. The program performs all the calculations for periodic regression as well as those for the Monte Carlo experiments to be described later. By appropriate selection of the starting and terminating scan wavelengths (step 2a above), periodic regression may be employed in a manner analogous to "band-pass filtering" (Craddock, 1957; Brier, 1961). The purpose of such filtering is to eliminate the contribution to the original time-series of those waveforms which are outside the wavelength range of interest. Such edited

time-series, when graphically displayed, illustrate the amplitude and phase characteristics of a particular frequency band. The process is the same as "tuning-in" to y o u r favorite radio station. In an ideal band-pass filter we would require that the "frequency response" be unity within the wavelength range of interest and zero outside this range. Unlike conventional digital filters, periodic regression provides us with a nearly ideal frequency response without any loss of data points at the ends of the time-series.

Evaluating the significance of frequencies obtained from periodic regression The results of periodic regression will indicate the existence of definite periodicities in a time-series. It has next to be decided which of these are the product of deterministic cyclic climatic mechanisms and which of these are the result of "noise" (i.e., those which stem from stochastic processes internal to the climatic system coupled with analytical errors). One way of making this decision is to compare the power spectra for independent parameters (e.g., oxygen isotope ratio and coarse fraction) in the same deep-sea core, or better, the same parameter in two or more widely separated cores. If a particular periodic c o m p o n e n t is consistently present and powerful, then it is not unreasonable to assume that it probably reflects the influence of a deterministic, as opposed to stochastic, process. However, the deterministic origin of periodic components can only be asserted with high confidence when their frequency and phase exhibit an acceptable correspondence with those predicted by a particular deterministic model, or the exact nature of the stochastic noise is known. In the absence of such models deterministic events may go unrecognized. Most students of paleoclimate think that the stochastic noise is characteristically "red", that is, in its power spectrum the power increases in some fashion with decreasing frequency {Gilman et al., 1963; Mitchell, 1965, 1976; Hasselman, 1976). This noise is envisaged as the product of numerous independent stochastic processes which are all internal to the climatic system (i.e., the atmosphere, cryosphere and oceans). Each is "characterized by its own unique time constant and each {adds) to the variability of climate on all time scales longer than that time constant (Mitchell, 1976, p.492)." These disparate stochastic processes arise from the differences in the rate and intensity of response, and m e m o r y capacity of the various components of the climatic system. The red stochastic noise may be said to exhibit the statistical property of "persistence". Such a p h e n o m e n o n arises from the Markovian dependency (i.e., lack of mutual independence} between the time-adjacent observations in a paleoclimatic record. According to Mitchell {1963, p.14), " . . . when persistence is present in the time-series, the higher frequencies of variation become damped because successive values of the series are restrained by their dependence on adjacent values from varying as widely as they otherwise could. At the same time the lower frequencies of variation are comparatively unaffected by persistence, which can impose only slight constraints on the relation

between values that are greatly separated in time. The effect of persistence is therefore to selectively filter out the higher frequencies, i.e., to 'redden' the spectrum." Superimposed on the stochastic background noise is narrow-band climatic variability (i.e., the power is confined to certain narrow frequency intervals) contributed by deterministic forcing processes external to the climatic system. The Milankovitch Theory of cyclic climate change is a well-known example of such a process. It must be kept in mind, however, that for a particular frequency band the portion of the total variance attributable to a deterministic process may be small compared to that contributed by the noise, where here we include in the noise, in the case of deep-sea core records, errors stemming from sampling and analytical procedures, and sample chronology. Thus, the signals generated by the deterministic processes may be deeply imbedded in the noise, and as such, are irresolvable. It should be clear from the preceding discussion that the contribution to a time-series from deterministic and stochastic processes consists in both cases of periodic components. Thus, in most circumstances, the two sources of climatic variability will be indistinguishable. The distinction can only be made in those cases where (1) the exact nature of the noise is known; (2) the frequencies are of substantial power and are consistently present in the timeseries of independent parameters in the same core, or the same parameter in two or more widely separated cores; or (3) the specific frequencies (with substantial power) and phase relationships predicted by the deterministic model are found to exist. The first approach is exceedingly difficult, if not impossible, in practice; and only the second and especially the third, offer us a reasonable solution to the problem. We must, however, in the latter two approaches, contend with the question of what constitutes a "substantial" power. Unless the frequencies of interest are characterized by "large" peaks in the power spectrum, we should be very cautious in accepting them as deterministic because noise also exists at these particular frequencies. What is needed then is some means of judging the substantiality of the powers of the observed periodic components: do they rise significantly above the general level of the surrounding noise? There is no formal statistical test of significance that we can invoke because such tests require that the e x a c t nature of the noise be known, or at least, assumed. We can, however, roughly estimate the persistence in our time-series and thus model the background level of the stochastic red noise. This estimate will aid us in passing judgment on the substantiality of powers. COMPENSATING F O R THE E F F E C T S OF PERSISTENCE

Our m e t h o d of estimating the background level of the stochastic red noise involves Monte Carlo experimentation. The idea here is to create a suitably large number, m, of simulated time-series each modeled after the original residual time-series (i.e., they have the same mean, variance and level of persistence), but otherwise consisting of random normal deviates (i.e., white

10 noise). Periodic regression is then applied to each simulated time-series for those q frequencies earlier found to be present in the original residual timeseries. A comparison for a given frequency, of the m " r e d " powers with the actual power for the paleoclimatic data should provide us with a means of judging the substantiality of the latter. This is n o t a statistical test of significance. The purpose of this procedure is to permit the investigator to evaluate, in a semi-quantitative manner, the power of the observed periodic components vis-a-vis an estimate of the " r e d " climatic noise. Consider the case where a given periodic component is found to be more powerful in the geologic timeseries than in a high percentage of the simulated time-series. Such a result strongly suggests rather than demonstrates (in a probablistic sense) that the periodic component observed in the geologic time-series is the product of a deterministic and not a stochastic process. The mathematical development is as follows: We assume the persistence is approximately a first order Markov process: X i = rl X i _ 1 + e i

(6)

where X i is the value of the dependent variable at time t i in the simulated time-series (i = 1, 2 , . . . , N); X i _ 1 is the value of the dependent variable at time t(i_l) ; r~ is the level of persistence and is a constant which is equal to the population value of the lag 1 serial autocorrelation coefficient of the Yi residual time-series (0 ~< r~ ~< 1); and e i is the error term and corresponds to an independently distributed (i.e., random) variable with mean Yr (approximately equal to zero) and variance (1 --r~) sy2 (Mitchell, 1966); where Yr and s~ are, respectively, the mean and variance of the residual time-series. The observation interval [i -- (i -- 1)] is assumed to be constant. We use as our estimate of rl the lag 1 autocorrelation for a Fourier series representation of the original residual time-series. The autocorrelation can be computed only from equally spaced observations and thus could not be calculated directly from our paleoclimatic records. It was necessary then to apply some form of interpolation to the original data in order to obtain the equally spaced observations. Accordingly, we sampled, at a uniform interval, the Fourier series representation of the residuals which was based on N / 2 [if N is even, otherwise ( N - - 1)/2] harmonic frequencies (i.e., in eq. 3 the frequency f1 assumed only integer values: 1, 2 , . . . , N / 2 ) . The ei's were computed in the following manner: e i = z i ~ / ( 1 --r, ) s~ + Yr 2

(7)

where :Yr and Sy are, respectively, the mean and variance of the original residual time-series; and z i is a random normal deviate produced by the application of the "central limit t h e o r e m " to the values obtained from the IBM 370 congruential pseudorandom number generator RANDU. According to the central limit theorem, the means 0i of large samples taken from a non-normally distributed population will tend toward a normal distribution. RANDU is used to provide sample values v1 from a uniformly distributed population where all values fall within the range 0 to 1. Thus,

11

zi

=

_

_

0 u

=

(8)

Ou

where uj is a single value generated by RANDU (k of these are used in the computation of each z i ) ; # u is the mean of the uj's (for our uniform distribution,~ u = l&); and o u is the standard deviation of the uj's (ou = ~ / x / - k ) . Eq.8 is the standardization formula and thus the N z i ' s are sequentially random but overall are normally distributed with a mean of zero and variance of one. We use k equal to 24. If rl equals zero, then X i = e i and the time-series consisting of the N X i values will constitute "white noise". In the power spectrum of white noise, the power at all frequencies will tend to be equal. As rl approaches 1 (i.e., the persistence or first order Markovian dependence increases), the power spectrum of X i progressively reddens (Mitchell, 1963, fig.l). It should be kept in mind that although in the present model of persistence X i is explicitly dependent upon X i _ l , it is also indirectly dependent on X i _ 2, X i _ a, etc. This is true because X i _ 1 is directly dependent on X i _ 2 which in turn is directly dependent upon Z i _ 3 and so on. From the above procedures, we generate m different simulated time-series, each consisting of N paired values of X i a n d t i. The values of the dependent variable X i are obtained from eq.6, and the corresponding ages t i are separated by a constant time increment. This increment is equal to the average spacing between the observations in the original residual time-series. The X i ' s are random normal deviates with a built-in persistence equal to rl. Further, the mean and variance of the X i ' s are equal respectively to ~ (approximately zero) and sy. 2 In our Monte Carlo experiments we use m equal to 99. The starting value X~ for the simulated time-series was set equal to e~ in eq.7. A simulated time-series with the same unequal spacing between observations as in the residual time-series could have been obtained by representing the former by a Fourier series which would then be sampled at a nonuniform interval. We think this extra step with the large amounts of additional computing time that it requires is not necessary in our case because there should be no significant differences in the power spectra of two simulated time-series which are identical in all respects except for a moderate difference in the spacing of their observations. In this study, the paleoclimatic records analyzed all have highly uniform observation intervals. Using the periodic regression eqs.4 and 5 we compute, for each of the m simulated time-series, powers for the same q frequencies observed in the original residual time-series. If the m red powers, for a given frequency, are in general small compared to the corresponding observed power in the paleoclimatic data, we consider the latter substantial and, hence, more likely the result of deterministic rather than stochastic processes. These comparisons of powers are graphically depicted in the lower portions (the PS values) of Figs.2, 4, 5, 6 and 7. PS is the percent of the m red powers which are smaller than the corresponding single observed power from the geologic time-series.

12 R ~ Rf 14

IIOC ,oc

1 9~

12

8~

I0

7( 6d

8

-

R~

Rf 100

,4

,[OO9o,z,4 t

!

~018

f°t

8O

C.dovis

70 I0

7O

F°° t

6O

:

5O

6"

4O 3O

42O

3 4

2O

I(

I0

,

i 0

, I 2

f 3~ 2q

5

6

.7 .8 #

0

r;

2 .3 ,4 .5 .6 .7 8 f

20 40

20 40

60 70

60 70

I 80 90

i

l

,

i

i

i

J

I 0

0 .I 2 3 4 5 6 7 .8 f

80 90

95

95

ps

PS

I00

o1

I00 ps

Fig. 2. Time-series analysis of ~ O t', T s and percent C. d a v i s i a n a (the PATCH ELBOW data set) from composite core RC11-120/E49-18. The power spectra published in Fig.5 of Hays et al. (1976) are contrasted with the results of periodic regression. The former are depicted as continuous curves and were computed using the method of the "truncated autocovariance function" (Blackman and Tukey, 1958). The data records have an average observation interval of 3200 years and span the period 0--487,000 years B.P. The wavelength interval scanned for periodic regression was 10,000--450,000 years. The frequency (f) corresponds to the number of cycles per 10,000 years. In the figures, the PS values are plotted using a logarithmic scale ( 1 0 Ps/100) in order to e m p h a s i z e the differences a m o n g the high PS values. Fig.8 illustrates s o m e o f the characteristics o f the results o b t a i n e d from the M o n t e Carlo experiments. S u p e r i m p o s e d on the periodic regression line spectrum for the o x y g e n i s o t o p e data o f core V 2 8 - 2 3 9 ( S h a c k l e t o n and O p d y k e , 1 9 7 6 ) are t w o sets o f curves. The solid and dashed curves o f each set correspond, respectively, to the mean p o w e r and the 8 0 t h percentile p o w e r c o m p u t e d from the m simulated time-series. For a given f r e q u e n c y , the mean p o w e r is the arithmetic average o f the m red p o w e r s and the 8 0 t h percentile p o w e r is that p o w e r b e l o w w h i c h o c c u r 8 0 percent o f the m red powers. The 8 0 t h percentile p o w e r is analogous to the one-sided 8 0 percent c o n f i d e n c e limit in classical statistics. The light-lined curves are based on m = 99. These are the same simulated time-series used to obtain t h e PS values in Fig.4. The heavy-lined curves are based on m = 5 0 0 . When m is very large, the mean and 8 0 t h percentile curves should closely approach the theoretical " p o p u l a t i o n " curves. Three i m p o r t a n t points can be m a d e from Fig.8. First,

13

RZn

IOOj, 40 /

ECCENTRICITY -

50 20 I0 0 0

iI i i I i i I .2 .3 .4 5 .6 7 .8

f=)w XI04

R2n I00~,

OBLIQUITY

20 / 15

I0 5 0

i

0 RZn I00

1

I

l

r

i

i

.2 .3 .4 .5 ,6 .7 .8 f=

~wX 104

PRECESSION

80 60 40 20 0

r i i r i 0 .I .2 .3 .4 .5 ,6 .7 .8 f= ~.~X tO4

Fig. 3. Periodic regression analysis of the astronomical data (precession, obliquity and eccentricity) of Berger (1977a) for the period 0--487,000 years B.P. The observation interval is uniform at 1000 years. The wavelength interval scanned for periodic regression was 10,000--450,000 years. The frequency (f) corresponds to the number of cycles per 1 0 , 0 0 0 years.

the two sets of curves are essentially coincident. Thus m = 99 is a sufficiently large number of simulated time-series for the Monte Carlo experiments. Second, the shape of the curves conforms nicely to the idea of a "red" noise power spectrum. Third, for both sets, the 80th percentile power curve does not parallel the mean power curve, the difference between the two becoming progressively inflated towards the low frequency end of the power spectrum. The ratio of the former to the latter, however, is nearly constant at 1.6 at all frequencies. Thus judgments on the deterministic or stochastic identity of the low frequency periodic components will be conservative. In periodic regression there is no a priori knowledge of h o w many periodic components will be resolved. The resolution process is terminated when some prespecified

14

R2n 100 90

8 0 le V 2 8 - 2 5 9

80 70 6O 50

L

40 30 2O I0 0

]

(X 103 YRS)

50 0 2O 40 6O 70 80 90 95 I00 PS

I00

150

200

250

300

350

400

450

500

550

600

I

11

Fig.4. Periodic regression analysis of 6018 from core V28-239 (Shackleton and Opdyke, 1976). The data record has an average observation interval of 5000 years and spans the period 5000--2,042,000 years B.P. The wavelengths (W) are resolved to the nearest 1000 years. The wavelength interval scanned for periodic regression was 10,000--500,000 years.

R2n tOO

-

/ 50-

S O 18 V I 6 - 2 0 5

40-

30-

20-

I0-

Oi 0 0 20 4,0 6O 7O 80 90

(X 103 YRS)

PS

i

50

I00

150 I

i

t

t

i

200

250

300

350

i

400

t

i

i

450

500

550

600

I

II'Jl

Fig.5. Periodic regression analysis of 8 0 TM from core V16-205 (Van Donk, 1976). The data record has an average observation interval of 16,300 years and spans the period 0--2,208,000 years B.P. The wavelengths (W) are resolved to the nearest 1000 years. The wavelength interval scanned for periodic regression was 25,000--600,000 years.

15 RZn

I00

% COARSE FRACTION V I 6 - 2 0 5 60

50

40

30

20

IO

0 (X 103 YRS)

i

0

i

5O I

0 20 4.0 6O 70 8O 9O 95 ~00 PS

I00

i

i

I

I

I

I

I

150

200

250

300

350

400

450

500

I

I

I

I

I

l

I

I

550

600

I

Fig.6. Periodic regression analysis o f "coarse fraction" from core V16-205 (Briskin and Berggren, 1975). The data record has an average observation interval o f 17,000 years and spans the period 0--2,244,000 years B.P. The wavelengths are resolved to the nearest 1000 years. The wavelength interval scanned for periodic regression was 25,000--600,000 years. R2n I00 413,000 YRS 80 -

9 6 , 0 0 0 YRS i 2 6 , 0 0 0 YRS

60 -

40

-

:~0

-

O~ 0

, ,fal,

------

ECCENTRICITY OBLIQUITY PRECESSION

i

50

I00

150

200 250 (X I03 YRS)

300

i

i

350

400

i

450

500

Fig.7. Periodic regression analysis of the astronomical data (precession, obliquity and eccentricity) of Berger (1977a) for the period 0--2,100,000 years B.P. The observation interval is uniform at 1000 years. The wavelengths (W) are resolved to the nearest 1000 years. The wavelength interval scanned for periodic regression was 10,000--600,000 years.

16 R2 I0 ,5018 V 2 8 - 2 3 9

9 8 7 6 5 4 3 2 i 0

~--, I0

20

50

,--, 40

50

60

1 70

--m=500 --m=99 i

i 80

i

i,t,fTl,,,i 90 I00

f Fig.8. Periodic regression line spectrum for 601~ from core V28-239 (Shackleton and Opdyke, 1976). The superimposed curves correspond to the mean power (solid) and 80th percentile power (dashed) for m = 99 (light line) and m = 500 (heavy line). The frequency (f) equals the number of cycles per 10,000 years. threshold is reached (for example, in Figs.3 through 8 the power threshold was 1 percent of the total variance; periodic components with powers less than this value were not resolved). APPLICATION OF PERIODIC REGRESSION TO PALEOCLIMATICDATA

The Milankovitch Theory of climatic change Many theories on the modulation of the Pleistocene "Ice Age" have been advanced and discarded with the notable exception of the Milankovitch Theory {i.e., the Astronomical Theory). What makes this theory so attractive to students of climate is that it alone is amenable to rigorous mathematical testing. According to the Milankovitch Theory (Milankovitch, 1941; Berger, 1977b), the a m o u n t and distribution of solar radiation received by the earth is controlled by the earth's orbital parameters, specifically, precession, obliquity and eccentricity. Of the three parameters, only changes in eccentricity can m o d i f y the total solar energy received. The other two contribute only to a different redistribution of that energy among the latitudes. At m a x i m u m eccentricity (~0.07) the perihelion--aphelion solar insolation difference is 28 percent, whereas at a low value of eccentricity (~0.017), such as exists at present, this difference is 6 percent (Berger, 1977b). The latitudinal redistribution of the total solar energy is principally a function of obliquity. The " t i l t " effect is most prominent at the high latitudes and negligible at the low latitudes. The precessional effect, on the other hand, is more uniform with respect to latitude. The precession cycle controls the timing of the seasons with respect to earth-sun distances. Because these distances are a function of eccentricity, the magnitude of the precessional effect is directly proportional to the eccentricity, being greatest at times of maximum eccentricity. Rather

17 than being viewed as a separate factor, eccentricity is then generally treated as a regulating influence on the precession cycle. It is instructive to note that Milankovitch did not anticipate the 413~000 year signal detected by Briskin and Berggren in 1975 and calculated by Andr~ Berger in 1977. Cyclical variations in the precession, obliquity and eccentricity result in concomitant periodic changes in the global insolation patterns and hence, the climate. Such astronomically induced climatic changes are widely thought to have modulated the ebb and flow of the Pleistocene ice sheets (Broecker and Van Donk, 1970; Imbrie and Kipp, 1971; Kukla, 1972, 1975; Briskin and Berggren, 1975; Hays et al., 1976; Berger, 1977b; Broecker, 1977; Emiliani, 1977). It is widely accepted (Emiliani, 1969; Imbrie and Kipp, 1971; Chappell, 1974; Briskin and Berggren, 1975; Weertman, 1976) that maximum accumulations of land-ice in the high latitudes tend to coincide with times of minimum summer-winter temperature contrasts. Vernekar (1972) and Berger (1977a) have computed the expected periodic secular variations in precession, obliquity and eccentricity for the last two and five million years, respectively. Significant variations are found to exist at numerous frequencies in each of the three parameters. The average wavelengths (in years) for the last two million years of the more important of these frequencies are: 413,000 and 96,000 (eccentricity); 41,000 (obliquity); and 24,000 and 19,000 (precession). These values were determined by applying periodic regression to the data of Berger (1977a, the actual raw data has not been published, but may be obtained from A.L. Berger). For a compilation of other cycles in the earth's orbital geometry see Figs.3 and 7.

Description of data In selecting data for time-series analyses great care must be taken. Of principal concern are those effects which can so distort a deep-sea record as to render the time-series analysis of the paleoclimate parameters meaningless. Both dissolution of the calcium carbonate and nonuniform sedimentation rates are c o m m o n l y recognized as the processes most responsible for the less than ideal deep-sea record and its plethora of chronostratigraphic problems. It is with these concerns in mind that we have endeavored to select the best possible Pleistocene deep-sea records available: (1) a composite of cores R C l 1 - 1 2 0 and E49-18 with the parameters 5 0 is (per mil enrichment of oxygen 18 over oxygen 16 relative to the PDB standard in the planktonic foraminifera Globigerina bulloides), Ts (estimated summer sea surface temperature), and percent of the radiolarian Cycladophora davisiana, a cosmopolitan species sensitive to changes in surface water structure (Hays, Imbrie and Shackleton, 1976), (2) core V16-205 with the parameters 5 0 ~8 (using Globigerinoides sacculifer; Van Donk, 1976) and coarse fraction (percent of the sample retained on a 74 p m sieve, Briskin and Berggren, 1975); core V28-239 with parameter 6 O~S. Core V28-239 (3 ° 15'N, 159°11'E) was raised from the equatorial western

18 Pacific waters of the Solomon Rise from a depth of 3490 m. This core is especially well known for its excellent oxygen isotopic record and constitutes a complete and detailed record of the Pleistocene. Dissolution and variable sedimentation have apparently not significantly affected the chronological control. Tropical North Atlantic core V16-205 (15°24'N, 43°24'W) situated near the gyral perimeter of the central waters was raised from a depth of 4045 m. V16-205 is the only other core with an oxygen isotopic record which spans the entire Pleistocene. The dissolution effects are negligible and the variation in sedimentation rate is within acceptable limits. During the interval of reversed polarity between the Brunhes and Jaramillo magnetic events (Briskin and Berggren, 1975, fig.2), the sedimentation rate appears to undergo a minor change, from an otherwise relatively constant trend of 0.5 cm per 1000 years. We think this apparent change may be due to the distortion introduced by a large radiometric error over a short core interval rather than a drastic change in sediment influx. Hsfi (1978) has reported on the excellent correlation between the Pleistocene chronology of V16-205 and the recognized climatic events of the Black Sea. The good standing of 50 is (an ice volume indicator) as a paleoclimatic parameter is well recognized. However, the relation between coarse fraction and climate is less clear. Broecker et al. (1958) suggested t h a t the average sedimentation rate, at any given locality in the Atlantic Ocean, tends to be the same for both the glacials and interglacials, but that the accumulation of fine sediments, terrigenous as well as carbonate, is proportionally greater during the glacials. Thus, coarse fraction minima should correspond to glacial maxima. The composite of cores RCl1-120 (43 ° 31'S, 79°52'E) and E49-18 (46°03'S, 90°09'E) constitutes a nearly ideal sedimentologic and paleontologic record for the Late Pleistocene. The data is of the highest quality and has been previously subjected to time-series analysis by Hays et al. (1976). Thus the paleoclimatic records available from these cores constitute a good data set for this study. Results and discussion

The reality of the Milankovitch Theory has been convincingly demonstrated by Hays et al. (1976). They found that the 96,000 year eccentricity cycle, and the major precession and obliquity cycles dominated the power spectra of their deep-sea paleoclimatic data. For their time-series analyses, Hays et al. used the more traditional approach of Blackman and Tukey (1958) -the method of the "truncated autocovariance function". The power spectra produced by this method {Hays et al., 1976, fig.5, PATCH ELBOW data) are shown as continuous curves in Fig.2. For purposes of comparison, we applied periodic regression to the raw PATCH ELBOW data provided us by John Imbrie. The resulting line spectrum and other information are presented in Fig.2. Periodic regression was also performed on Berger's (1977a) astronomical data for the same period of

19 TABLE I S u m m a r y of t h e d o m i n a n t cycles in t h e p a l e o c l i m a t i c d a t a o f Hays et al. ( 1 9 7 6 ) a n d t h e a s t r o n o m i c a l d a t a o f Berger ( 1 9 7 7 a ) f o r t h e period: 0 to 4 8 7 , 0 0 0 years B.P. 1 Astronomical data 2

PATCH ELBOW 3 (Hays et al., 1 9 7 6 )

PATCH ELBOW 2 (this study)

Ts

50 ~

%C.d.

Ts

8 0 '~

%C.d.

94

106

1224

104

107

113

40

43

43

42

41

23 --

24 19.5

23--24 20

24 20

Eccentricity 101

Obliquity 41

43

Precession 23--24 19

24 --

24 --

1 F o r details o n t h e o b s e r v e d cycles see Figs.2 a n d 3. T h e cycles are in t h o u s a n d s of years. C o m p u t e d using p e r i o d i c regression. T h e w a v e l e n g t h o f t h e cycles is resolved to t h e n e a r e s t 1 0 0 0 years. 3 C o m p u t e d using the m e t h o d o f the " t r u n c a t e d a u t o c o v a r i a n c e f u n c t i o n " ( B l a c k m a n a n d Tukey, 1958). 4 122 is t h e r e p o r t e d value for t h e d o m i n a n t peak in t h e C. davisiana s p e c t r u m b u t an e x a m i n a t i o n o f Fig.5 of Hays et al. ( 1 9 7 6 ) or Fig.2 o f this p a p e r reveals t h e a c t u a l value is closer t o 145 t h o u s a n d years. This is f u r t h e r c o n f i r m e d b y t h e p e r i o d i c regression results.

time: 0 to 487,000 years B.P. (Fig.3). The results and interpretations of these analyses for the major cycles are summarized in Table I. Many of the minor cycles seen in the periodic regression results for the astronomical data (Fig.3) are seen as well in the results for the paleoclimatic data (Fig.2). The greater correspondence between the cyclic components of the astronomical and paleoclimatic records, as demonstrated by periodic regression, further substantiates the validity of the Milankovitch Theory of climatic change. The results obtained from periodic regression and the method of time-series analysis used by Hays et al. are comparable with periodic regression providing better wavelength resolution in the low to intermediate frequency portion of the power spectrum. Periodic regression with cyclic descent is an excellent time-series technique for resolving the discrete cyclic components actually present in the data. It contrasts sharply with the method used by Hays et al. which resolves only arbitrary frequency bands. In search of broader periods we applied periodic regression to the 6018 (Fig.5) and coarse fraction (Fig.6) records of core V16-205 and the 6 0 ~8 (Fig.4) record of core V28-239. These paleoclimatic records comprise a wellconstituted data set, consisting of two independent parameters from the same core and the same parameter in two different, widely separated cores. Such a data set provides the opportunity to identify the climatic effects of deterministic forcing mechanisms. All three records span the entire Pleistocene. Detection of the 413,000 year eccentricity cycle can only be achieved in long cores with a minimum time span of approximately 2 million years.

20 Presented in Fig.7 are the results of the time-series analysis for Berger's (1977a) astronomical data (precession, obliquity and eccentricity) for the period 0 to 2,100,000 years B.P. In Figs.4, 5 and 6 we are primarily interested in the spectral expression of the 96,000 and 413,000 year eccentricity cycles. Because of the large observation intervals, problematic chronological control and the relatively low sedimentation rates in these cores, resolution of the shorter wavelength cycles predicted by the Milankovitch Theory is not expected. Our purpose at this stage of this research is to report the strong evidence of an eccentricity effect in the paleoclimatic record of deep-sea cores. The possible eccentricity related wavelengths are as follows: in core V16-205 -434 X 103 years and 101--105 X 103 years f o r 5 0 is and 435 × 103 years and 101 X 103 years for coarse fraction; and in core V28-239 -- 412 × 103 years and 101--105 X 103 years for 6 0 is . The PS values in the lower portion of each plot indicate that these periodic components, in almost every case, possess substantial power and are thus not likely to be artifacts of the red climatic noise. The range of 412--435 thousand years for the longer cycle constitutes an excellent agreement between records when one considers the many possible errors in the data. The fact that periodic components with wavelengths of roughly 103,000 and 420,000 years are consistently present in independent parameters and different cores of the Pacific and Atlantic Oceans is strongly suggestive of a deterministic origin. The results of this study affirm the Astronomical Theory of Climate. Work remains to be done on the phase relationships between the paleoclimatic and eccentricity record. Cross-spectral analysis between astronomical parameters and core variables should help to elucidate amplitude and phase relations at various frequencies (T.B. Start, personal communication). Indeed the question as to whether the response model for different paleoclimatic parameters is linear or not remains to be studied. In conclusion, the preliminary results in this study strongly suggest a deep-sea climatic expression of the dominant 413,000 year eccentricity cycle, first detected in tropical North Atlantic core V16-205 by Briskin and Berggren (1975, fig.10, p.187). ACKNOWLEDGMENTS We are indebted to Dr. Robert Kuhn (University of Cincinnati) for his expert guidance and generous contributions in the development of the periodic regression methodology. We are indebted to Dr. John Tukey (Bell Laboratories and Princeton University) for his help and constructive criticism of the manuscript. We are grateful to Dr. Peter Bloomfield (Princeton University) and Dr. Murray Mitchell (National Oceanic and Atmospheric Administration) who generously discussed and helped us to clarify some of the mathematical questions. Their suggestions were very helpful. We thank Dr. Rayner for critically reviewing the manuscript in the early phases of this study.

21

Dr. Andre Berger (University of Louvain) gave his permission to use his astronomical data. Dr. Steve Streeter and Dr. George Kukla (both of the Lamont-Doherty Geological Observatory of Columbia University) provided the astronomical data on magnetic tapes. The senior author, Madeleine Briskin, is indebted to Thomas B. Starr (University of Wisconsin, Madison) for his meticulous review of the final manuscript. His expert advice helped to strengthen the final version of the manuscript. The senior author is grateful to Dr. Alan Hecht (National Science Foundation) for his constant interest, vision and words of encouragement; University Dean Albert Yates and Assistant Dean Frank Tepe were very helpful and provided financial support which made this research possible. An award from the University of Cincinnati Research Council was much appreciated. The Department of Geology Research Fund provided additional funds. Mrs. Jean Carrol kindly typed the manuscript. REFERENCES Berger, A.L., 1976. Obliquity and precession of the last 5,000,000 years. Astron. Astrophys., 51: 127--135. Berger, A.L., 1977a. Support for the astronomical theory of climatic change. Nature, 259 (5623): 44--45. Berger, A.L., 1977b. Power and limitations of energy balance climate model as applied to the astronomical theory of paleoclimates. Palaeogeogr., Palaeoclimatol., Palaeoecol., 21: 227--235. Blackman, R.B. and Tukey, J.W., 1958. The Measurement of Power Spectra from the Point of View of Communications Engineering. Dover, New York, N.Y. Bloomfield, P., 1976. Fourier Analysis of Time-Series: An Introduction. Wiley, New York, N.Y. Brier, G.W., 1961. Some statistical aspects of long-term fluctuations in solar and atmospheric phenomena. In: F.N. Furness (Editor), Solar Variations, Climatic Change and Related Geophysical Problems. New York Acad. Sci., New York, N.Y., pp.173--187. Briskin, M. and Berggren, W.A., 1975. Pleistocene stratigraphy and quantitative paleooceanography of tropical North Atlantic core V16-205, In: T. Saito and L. Burckle (Editors), Late Neogene Epoch Boundaries. Micropaleotology Press, Museum of Natural History, New Y~rk, N.Y., pp. 167--198. Broecker, W.S., 1966. Absolute dating and the astronomical theory of glaciation. Science, 151: 299--304. Broecker, W.S. and Van Donk, J., 1970. Insolation changes, ice volumes and the O 18 record in deep-sea cores. Rev. Geophys. Space Phys., 8: 169. Broecker, W.S., Turekian, K.K. and Heezen, B.C., 1958. The relation of deep-sea sedimentation rates to variations in climate. Am. J. Sci., 256: 503--517. Chappell, J., 1974. Relationship between sea levels, O 18 variations and orbital perturbations during the past 250,000 years. Nature, 252: 199--202. Craddock, J.M., 1957. An analysis of the slower temperature variations at Kew Observatory by means o f mutually exclusive band-pass filters. J, R., Statistical Soc. (ser. A., general), IV, 120: 387--397. Emiliani, C., 1955. Pleistocene temperatures. J. Geol., 63: 538--578. Emiliani, C., 1969. Interglacial high sea levels and the control of Greenland ice by the precession of the equinoxes. Science, 178: 398--401. Gilman, D.L., Fuglister, F.J. and Mitchell, J.M., 1963. On the power spectrum of "red noise". J. Atmos. Sci., 20: 182--184.

22 Hasselman, K., 1976. Stochastic climatic models, I. Theory. Tellus, 2 8 : 4 7 5 - - 4 8 5 . Hays, J.D., Imbrie, J. and Shackleton, N.J., 1976. Astronomical theory of ice ages confirmed. Science, 194 (4270): 1121--1132. Hsi], K.J., 1978. When the Black Sea was drained. Sci. Am., 238: 52--63. Imbrie, J. and Kipp, N.G., 1971. A new micropaleontological method for quantitative paleoclimatology: application to a late Pleistocene Caribbean core. In: K.K. Turekian (Editor), Late Cenozoic Glacial Ages. Yale Univ. Press, New Haven, Conn., pp.71--182. James, W.R., 1966. FORTRAN IV program using double Fourier series for surface fitting of irregularly spaced data. Kansas Geol. Surv. Comput., Contrib. 5. Kukla, G.J., 1975. Missing link between Milankovitch and climate. Nature, 253: 600-603. Kukla, G.J. and Kukla, H.J., 1972. Insolation regime of interglacials. Quat. Res., 2: 412-424. Milankovitch, M., 1941. Canon of Insolation and the Ice-Age Problem. K~niglich Servische Akademie, Beograd. Trans. Israel Program for Scientific Translations. U.S. Dept. Commerce and National Science Foundation, Washington, D.C. (in German). Mitchell, J.M., 1963. Some practical considerations in the analysis of geophysical timeseries (unpublished manuscript). Paper given at the 44th Annual Meeting of the American Geophysical Union, Washington, D.C. Mitchell, J.M., 1966. "Stochastic models of air--sea interaction and climatic fluctuation. In: J.O. Fletcher (Editor), Proceedings of the Symposium on the Arctic Heat Budget and Atmospheric Circulation; The Rand Corp. Memorandum RM-5233-NSF. Mitchell, J.M., 1976. An overview of climatic variability and its causal mechanisms. Quat. Res., 6: 481--493. Schickendanz, P.T. and Bowen, E.G., 1977. The computation of climatological power spectra. J. Appl. Meteorol., 16 (4): 359--367. Shackleton, J.N. and Opdyke, N.D., 1976. Oxygen-isotope and paleomagnetic stratigraphy of Pacific core V28-239 Late Pliocene to Latest Pleistocene. Geol. Soc. Am. Mem., 145: 449--464. Van Donk, J., 1976. O ~s record of the Atlantic Ocean for the entire Pleistocene epoch. Geol. Soc. Am. Mere., 145: 147--163. Vernekar, H.D., 1972. Long-period variations of incoming solar radiation. Meteorol. Monogr., 12(34). Walker, A.M., 1971. On the estimation of a harmonic component in a time-series with stationary independent residuals. Biometrika, 58 (1): 21--36. Weertman, J., 1976. Milankovitch solar radiation variations and ice-age ice-sheet sizes. Nature, 261: 17--20.