Milankovic Pseudo-cycles Recorded in Sediments and Ice Cores Extracted by Singular Spectrum Analysis

Milankovic cycles describe the changes in the Earth's orbit and rotation axis and their impact on its climate over thousands of years. Singular Spectrum Analysis (SSA) is a signal processing method that is best known for its ability to find and extract pseudo-cycles in complex signals. In this short paper, we propose to apply it to three time series that have been proposed as geological reference time scales, in order to retrieve, compare and identify their Milankovic periodicities: (1) LR04, a stack of Plio15 Pleistocene benthic microfossil records (Lisiecki and Raymo, 2005), (2) the CO2 and CH4 records from the Vostok ice core (Petit et al, 1999) and (3) the long-term orbital solution La04 for the insolation of Laskar et al (2004). The Vostok CO2 and CH4 series share the first 7 SSA components, three main ones (98, 104, 39 kyr), and four smaller ones (18, 22, 65, 180 kyr). CO2 displays a component at 28kyr and a doublet at 61 and 62 kyr. CH4 displays a doublet near 50 kyr. 18/22 ky is a precession doublet, 62 kyr an insolation component, and 95/105 kyr an insolation/eccentricity doublet. The 49/50 kyr doublet in CH4 is 20 not found in the orbital model. The SSA results for the La04 orbital solution are in excellent agreement with the values obtained by Laskar et al (2004). Four SSA components of obliquity are almost identical (rounded figures are 41, 54, 29 and 39 kyr). As far as eccentricity is concerned, the first five components are 404, 95, 124, 99, and 132 kyr. The next components are not found in our list of components for eccentricity, but they are in the SSA of insolation, at 2338, 970, 488 and 684 kyr. With more than 20 components, the LR04 stack is the richest series. In order of decreasing amplitude, one 25 encounters 41, 95 and 75 kyr components. Next are smaller 39.5 and 53.6 kyr components, and a 22.4 kyr component. One recognizes one of the two main precession components, the doublet of obliquity components, a line at 47.4 kyr that is not found in any of the other spectra, and a doublet at 53.6 and 55.7 kyr, corresponding to the line at 54 kyr found in all four orbital quantities. Next comes a line at 63.6 kyr that may correspond to a line in insolation, CH4 and CO2. Then come components from eccentricity variations at 75.2, 94.5, 107.2, 132.1, 198.6 and 400.9 kyr. The remaining components of 30 LR04 show up in La04. The “elusive ~200 kyr eccentricity cycle” of Hilgen et al (2020) is actually present in all three series, in the La04 orbital model as a 195±6 kyr component of eccentricity and in LR04 as a 198.6±5.6 kyr component. Finding not only the main expected Milankovic periodicities but also many “secondary” components with much smaller amplitudes https://doi.org/10.5194/cp-2021-126 Preprint. Discussion started: 1 November 2021 c © Author(s) 2021. CC BY 4.0 License.

Detailed responses to the reviewer's criticisms, including two more extensive notes, one on Milankovitch theory, the second on the SSA method.
1: I have difficulties to understand the main topic of this manuscript: Either the main message is that paleoclimatic records contain astronomical periodicities, something which is well-known for almost 50 years. Or, more probably, the aim is to present a "new" spectral analysis method that could detect extremely faint periodicities in geological records. In the first case, I would reject the paper on the ground that it brings nothing new.
Of course the paleoclimatic records contain astronomical periodicities; that has been known at least since Milankovitch (1920). But the question here is to use SSA to determine how many of these periodicities (SSA components) are present and what are their respective amplitudes (contributions to the original signal). And we succeed in doing it and our results are indeed in part new (number, periodicity content and amplitude (contribution to the total variance) of main contributing (quasi-)periodicities (components) determined using SSA, more complete than found in previous publications to our knowledge). See also Note I.

2:
But the authors are not even citing the famous foundational Hays et al (1976) paper in their reference list, which is, to say the least, very peculiar and which suggest that they may not be familiar with this topic and its extremely abundant literature.
We do know and cite Hays et al (1976) in our text, but we have indeed forgotten it in the reference list. We will correct this upon revision. We note that that paper only finds 7 periodicities (in ) when we identify 5 more with periods longer than 40 kyr (in both CH4 and CO2 in the same Vostok core): 50, 66, 100, 180 kyr. Our reference list and previous published papers show that we have some familiarity with the topics involved in our paper.

3:
Besides (see below), it appears that the authors have even missed the main point of Milankovitch theory, which is based on summer insolation, and not global mean insolation...
Wrong. Note I below answers this and other comments at some length.

NOTE I (ON MILANKOVIC THEORY)
Milankovich's thesis (1920, in  We refer to Milankovitch's original papers (that are in French; f.i. Part II, chapter 1, page 169, entitled « Distribution de la radiation solaire à la surface du globe terrestre et son climat mathématique »). In short, Milankovitch does not mention summer insolation but global insolation.
We now follow in detail Milankovitch's reasoning.
The last sentence of the reviewer is indeed encountered in many if not most papers on present and past climate. We claim that it is both wrong and rather disconcerting.
Why disconcerting? Let us use the case of marine sediments from Lisiecki and Raymo (2005); their time sampling, which is not regular, has a mean of 2500 years. What is then the physical meaning of considering only summer, i.e. 3 months of the year, when the time series under analysis has a sampling rate of several thousand years and the goal is to determine astronomical parameters the smallest of which is 20,000 years ?
Why wrong? Let us go back to the writings of Milankovitch (1920Milankovitch ( , 1948. We give a link to the full texts, the former in French, the latter in English and German. The latter is not scientifically innovating with respect to the former as he himself states (left column, page 414, second paragraph); « In dem Buche Théorie mathématique des phénomènes thermiques produits par la radiation solaire (Paris 1920)  regions where the Sun doesn't set, Milankovitch shows that the daily variation is given by: Milankovitch represents these two equations for the mean radiation as an axonometric figure ( Figure 13 page 182) that we reproduce below.
It is therefore quite clear that Milankovitch does take into account all latitudes. It is this very same Figure 13 that is referred to as Zackenkurven in Milankovitch (1948).
We next come to the topic that Milankovitch calls "Variations séculaires des éléments 6/28 astronomiques et le problème paléoclimatique" (page 221, Chapter 2). We attempt to follow Milankovitch's reasoning in his thesis as closely as possible. At the beginning of page 221, the author recalls the 4 parameters that influence the Earth's insolation: the eccentricity of the planet's orbit, the longitude of the perihelion, the general precession and the inclination of the planet's axis of rotation.
Thanks to formulas derived by Stockwell and to his own equations (12) and (13)  Milankovitch's study of the secular variation of insolation leads to Figure  This means that if the rotation axis "straightens up", then it will be colder in the Summer (the Earth being farther from the Sun) and warmer in Winter (the Earth being closer to the Sun), as is indeed observed currently.
We now come to chapter 56, part 2 (page 246), the most relevant to our study and to this response to the reviewer, entitled "Théories astronomiques des époques glaciaires". Milankovitch first presents his master Adhémar's theory, that he writes is "sous l'empire des idées de Fourier" (page 247) and that he finally rejects. So does he reject Croll's theory, of which he acknowledges only the idea of an astronomical forcing. In this rejection, there is a remarkably simple development, grounded on observations. For instance, he writes on page 251 that beyond the 66th parallel ( ), a zone that in Milankovitch's theory obeys a second set of equations, "Par conséquent à ces latitudes, on peut obtenir presque les mêmes effets par la variation de l'obliquité 7/28 de l'écliptique que par la variation de la durée des saisons".
This is an important remark, as soon as one is concerned by geochronology at these latitudes. This is seen in the data as he clearly spells out on page 252: "Ceci ressort également de la Figure   21

END OF NOTE I
________________________________________________________________________________ 4: In the second case, I would also reject the paper since it does not show any information on the key question of spectral analysis: statistical relevance.
We have provided the key reference for all practical and mathematical aspects of SSA (Golyandina, N., and Zhigljavsky, A.,: Singular Spectrum Analysis for time series, Vol. 120, Berlin, Springer, 2013), in particular statistical relevance. We do give uncertainties for all component periods in Table 1. We devote Note II to SSA and statistics to refute the reviewer's criticism.
This note recalls two methods, one of which (SSA) was used in this paper. Its aim is to recall the bases of frequency analysis, and to compare them to SSA. This note is a summary, but a rather extensive one, given the fact that SSA does not seem to have enjoyed the success one could have imagined. Also, this is to make sure that the reviewer shares our understanding and competence in using the method.

NOTE II (ON SPECTRAL ANALYSIS , STATISTICAL ANALYSIS AND SSA)
There are several definitions of spectral analysis (e.g. Claerbout, 1976 or Bracewell andBracewell, 1986). We choose the simplest and most basic one (cf. Stoica and Moses, 2005, chapter 1.1 page 1): "From a finite record of stationary data sequence, estimate how the total power is distributed over frequency". Spectral analysis, the most ancient branch of harmonic analysis, of whom Fourier series (Prony, 1795, Fourier, 1822 are the most ancient representative, consists in projecting any stationary (in a strict sense, statistically) signal on an orthogonal basis of sine functions. Its limitations are well known (e.g. Shannon and Burks, 1951;Kay and Marple, 1981).
The discrete Fourier transform ( ) of any time series sampled every and represented by N samples, is given by the following equation : The sampling rate  in the frequency domain must satisfy Shannon's dual law ( ) where T is the duration of the signal. k is an integer describing the position of a given frequency in the dual space, whose maximum value must be at least equal to N (second Shannon criterion). n is an integer describing the position of a sample in and W is defined as: The minus sign in (02) is a convention and must be replaced by + in the reverse Fourier transform.
One sees that projection (01) 12/28 (03) can be written , which is the traditional expression of a direct linear problem (Menke, 1989;Tarantola, 2005), in which G is a Vandermonde matrix with ascending diagonal. G possesses remarkable symmetry properties imposed by the definition of W (all complex exponentials being equal to each other modulo 2. Let us, at this point, emphasize that a spectral analysis and a SSA are two very different things, since the two matrices underlying the two methods are completely different. We next discuss the link between spectral analysis and statistics, or rather probabilities. The convolution between two functions f(t) and g(t) is expressed by: This is also the general formulation of the continuous Fourier transform ( ).
Convolution appears in probability theory in the following way. Let  and  be two independent random variables with respective probability densities a() and b(). The probability density g() of the variable  is given by the convolution: This property, associated with the central limit theorem, explains why the normal law occupies a special place in statistics.
We now recall that the autocovariance (or selfcovariance) of a stochastic process allows one to characterize the linear tendencies that exist within this process. Given a signal x(t), the self covariance at time  is the mean product: In spectral analysis, one hopes to write x(t) as a sum of sine functions (see equations (01) and/or (03)), and it can then be written: In that case, (06) becomes: (08) contains the same frequencies as the signal x(t). Amplitudes are squared and phases ( ) have disappeared. The Fourier transform of (08) is called the spectral density S(f).
This is a very short summary of the links between spectral analysis and statistical or probability theory. If the reviewer wants to know more about the statistics of our analysis, he/she can find the values of amplitudes between parentheses in our Table (since energies are conserved from data space to dual space, whichever decomposition is involved).
We now come to the algorithm of Singular Spectral Analysis. Paradoxically, this method is not a spectral analysis in a strict sense. It has been popularized in the geophysical community by Step 1 (embedding step) is divided into K segments of length L in order to build a matrix X with dimension with 1 will condition our decomposition. This is the first "tuning knob". Integrating in X yields a Hankel matrix: whereas matrix (03) is a Vandermonde matrix. Comparing the two methods is meaningless. 14/28 Embedding, the first step in a SSA, consists in projecting the one-dimensional time series in a multidimensional space of series (X1, … , Xk), such that vectors belong to , where . This definition was introduced by Mané (1981) and Takens (1981), with the aim to build a space in which strange attractors could be described correctly, often a Banach space. A strange attractor is an object whose dynamical properties may change and evolve into chaos (that is by nature non-linear). The parameter that controls the embedding is , the size of the analyzing window, an integer between 2 and . The Hankel matrix has a number of symmetry properties: its transpose , called the trajectory matrix, has dimension . Embedding is a compulsory step in the analysis of non linear series. It consists formally in the empirical evaluation of all pairs of distances between two offset vectors, delayed (lagged) in order to calculate the correlation dimension of the series. This dimension is rather close to the fractal dimension of strange attractors that could generate that type of series. In this case, it is advised to select for the size of window L very small values (resp. very large K). On the contrary, for SSA, L must be sufficiently large, so that each vector contain an important part of the information contained in the initial time series . From a mathematical point of view, one must work in the frame of Structural Total Least Squares for a Hankel matrix (e.g. Lemmerling and Van Huffel, 2001). This is quite different from the fractal dimension mentioned above. Another benefit from considering very large L is the possibility of considering sub-vectors Xi as independent sub-series with distinct dynamics, and thus to be able to identify common characteristics in collections of these sub-series.
Step 2 (decomposition in singular values -SVD) SVD (Golub et Reinsch, 1971) of non-zero trajectory matrix (dimensions ) takes the shape: where the eigenvalues of matrix are arranged in order of decreasing amplitudes. Eigenvectors Ui and Vi are given by : The Vi form an orthonormal basis and are arranged in the same order as the i. Let Xi be a part of matrix X such that :

15/28
Embedding matrix X can then be represented as a simple linear sum of elementary matrices Xi. If all eigenvalues are equal to 1, then decomposition of X into a sum of unitary matrices is : d being the rank of X ( ). SVD allows one to write X as a sum of d unitary matrices, defined in a univocal way.
Let us now discuss the nature and the characteristics of the embedding matrix. Its rows and columns are sub-series of the original time series (or signal). The eigenvectors Ui and Vi have a time structure, and they can be considered as a representation of temporal data. Let X be a suite of L lagged parts of and ( ) the linear basis of its eigenvectors. If we let: with , then (13) can be written: that is for the j th elementary matrix: where is a component of vector . This means that vector is composed of the ith components of vector . In the same way, if we let: we obtain for the transposed trajectory matrix : that corresponds to a representation of the K lagged vectors in the orthogonal basis ( ).
One sees why SVD is a very good choice for the analysis of the embedding matrix, since it provides two different geometrical descriptions.
Remark 1: there are very strong common features between performing a SVD of the trajectory matrix, as is done in SSA, and multivariate analyses in principal components (PCA) or Karhunen-

16/28
Loève (KL) decompositions, that are commonly found in time series analysis. SSA differs through the nature of its trajectory matrix, a Hankel matrix with a particular structure : its rows and columns are sub-parts of the signal to be analyzed. They both have a meaning, obviously a temporal physical meaning. Such is not the case for PCA or KL.
Remark 2 : In general, the orthonormal basis (Ui) associated to the trajectory matrix and obtained through SVD can be replaced by any orthonormal basis (Pi). In that case, equation (14) becomes with , for those who want to perform a Fourier analysis with an orthonormal basis made of sine functions. A classical example of an alternate basis is that of the eigenvectors of a self-covariant matrix (Toeplitz SSA).
Step 3 (reconstruction) As we have seen, Xi matrices are unit matrices, and (as in the classical approach) one can "regroup" these matrices into a physically homogeneous quantity (or energetically homogeneous, etc…). This is the second "tuning knob" of SSA. In order to regroup the unit matrices, one divides the set of indices into m disjoint subsets of indices .
Let I be the grouping of p indices of ; because (14) is linear, then the resulting matrix that regroups indices I can be written: This step is called regrouping the eigen-triplets (, U and V). In the limit case , becomes exactly (14), and we find again the unit matrices.
Next, how can one associate pairs of eigen-triplets? This means separating the additive components of a time series. One must first consider the concept of separability.
Let be the sum of two time series  and  , such that for any .
Let L be the analyzing window (with fixed length), and the embedding matrices of series ,  and  . These two sub-series are separable (even weakly) in equation (14) if there is a collection of indices such that , respectively if there is a collection of indices such that In the case when separability does exist, the contribution of for instance corresponds to 17/28 the ratio of associated eigenvalues ( ) to total eigenvalues ( ) .
Still in the case of (14), let be the set of indices corresponding to the first time series, with corresponding matrix . If both this matrix and that corresponding to the second time series ( ) are close or identical to a Hankel matrix, then the time series are approximately or perfectly separable. So, regrouping SVD components can be summarized by the decomposition into several elementary matrices, whose structure must be as close as possible to a Hankel matrix of the initial trajectory matrix (this is true on paper only, in reality things are much more difficult).

END OF NOTE II
________________________________________________________________________________ 5: Indeed, the authors seem to be unaware that white noise contains "all frequencies".
Ironic or somewhat insulting. What is the use of this remark in the review? We do not understand. Either the reviewer considers that the data (which are not ours) are actually white (or pink, or…) noise, or he/she has not seen the Table with the values (between parentheses) of absolute amplitudes of the extracted components that are clearly not identical (hence one cannot speak of white noise).
6: Therefore, finding frequencies in a signal is NOT (has never been) the key point of spectral analysis.

20/28
Did we understand this comment? Spectral analysis is not about the determination of the frequency contents of a time series? It may not always be "the key point" but it certainly is "a key point". "Spectral analysis" is often called "frequency analysis".

7:
The only interesting question is to estimate their statistical significance.
The "only" interesting question? When projecting a signal on an orthogonal basis, such as in Fourier spectral analysis, the statistics are certainly not the main goal.
___________________________________________________________________________ 8: But there is simply no mention at all of statistics in this paper.
Indeed, in this paper we do not need a statistical analysis, except for the rough estimate of uncertainties of the central frequency of individual SSA components. If there were large numbers of spectral peaks/SSA components and only some could be attributed to planetary frequencies, we could not decide which ones were significant. But when most observed frequencies are recognized as naturally present in the series that is analyzed, it is unreasonable to assign it to chance.
As we have recalled above in Note II, classical spectral analysis and SSA are not statistical analyses, even if one can link these two different branches of mathematics. If the reviewer is asking about the meaning of periodicities we identify, then he/she must have the same concern about the work of Hays et al (1976) or Vautard and Ghil (1989). We actually share part of their lists of periods; not only do we confirm but we extend their lists.

9:
SA is now a very classical spectral analysis method. It was largely developed in the 1980s (Broomhead and King, 1986) and was used for paleoclimatic studies since then (Vautard and Ghil, 1989).
We know and cite this useful paper by V and G (the reviewer?). Indeed SSA has been used in climate studies but not in a number of fields and series we are now analyzing. If the method is so classical, why the questions on the statistics of our results? Could the reviewer help us by giving us the references to papers that perform statistical tests on the components of a SSA analysis? In any case, we have discussed this at length in Note II above. Now, if the question is about the small 21/28 absolute percentages we indicate in our Table, asking about the reality of the components we extract, there are two answers. The first is that this percentage is the fraction of the total signal. But the largest (first) component is in general the trend, which is non linear. It would be legitimate to de-trend the series, since we are interested in periodic and quasi-periodic components, and the amplitudes of all the other components would be correspondingly increased. The second answer is that indeed several of these components are very small compared to the general trend: to which extent could FA as well as SSA (i.e. SVD) be able to separate information in these signals? This is not anymore a question of algorithm (stricto sensu), but of numerical resolving power and precision of computation. This means of course computer, operating system and language (eg. (which is on the order of 10 -16 ) and the fact that our results are exact, since they rely on mathematical functions (SVD, Hankel, FFT, etc … ) implemented by the 5000 mathematicians, numericists and computer scientists that work in the company (Mathworks) that sells Matlab.
10: It was even applied to the Vostok record as early as 1994 (Yiou et al., 1994).
True and we do not contradict the results, we complement them.

11:
Again, it is strange not to find these references in a scientific paper on SSA applied to paleoclimate records in general and to the Vostok data in particular. Yiou et al. (1994) find with SSA the same components found by Petit and Jouzel (1977) using Fourier analysis. The papers mentioned by the reviewer, e.g. Vautard and Ghil (1989), do not describe the full mathematical apparatus of SSA, but are a very good introduction to the method.
We can make a comparison with the discovery of wavelets by Jean Morlet (Morlet et al., 1982). A rudimentary tool first developed by a geophysicist, then developed much more fully by mathematicians who were specialists in time series analysis (eg. Grossmann and Morlet, 1984), up to the pyramidal algorithm of Mallat (2009) that is now (as is SVD) in the top ten of the most used algorithms in research (cf. Cipra, 2000). In the same way, the work of Vautard and Ghil barely touches on the concept of separability of information, as we defined it in Note II above. This is the reason why we systematically cite Golyandina and Zhigljavsky's (2013) book. We have detected more periods because, as stated in our paper, we have used a variant of SSA, called iterative SSA 22/28 (or iSSA). This implies changing the size of the analysing window at each step of detection and extraction of a period (component) so that equation (14) remains as close as possible to a Hankel matrix. Indeed, there are a large number of variants of SSA: iterative SSA (e.g Golyandina and Zhigljavsky, 2013), varimax Rotation SSA (e.g. Portes et Aguirre, 2016), multitaper SSA (Dettinger et al., 1995), circulating SSA (e.g. Bogalo et al., 2021), etc … All aim to ensure optimal separability as a function of the information carried by a given signal (time series). Vautard and Ghil (1992) had performed the simplest, most basic kind of SSA, so the differences in results (improvements) are not surprising.
difficult to find any novelty in this paper.
The many papers that have all found the same set of 4 or 5 periodicities need not be all cited!
The point of our paper is that we find more components and can identify them. Also, this is a way of testing our results on well-known series before applying them to new complementary and/or longer, better sampled series. We find more "Milankovitch"-type periods.
13: Spectral analysis is a branch of statistics.
No! It can use tools from the field of statistics, it is not an (ancillary) branch of this field. One can use a number of statistical tools with Fourier, not with wavelets or SSA, or PCA. One can make statistical inferences with Fourier (autocorrelation, DSP). SSA is a decomposition on an ad-hoc basis, not sinus functions (and therefore not Laplacian). Same for wavelets. See Note II.
14: Finding many frequencies (real or spurious) in a signal is known to be extremely easy, since the discovery of the Fourier transform. In particular, finding astronomical frequencies and their numerous harmonics in paleoclimatic signals is certainly not surprising.
Yes but whereas Fourier creates harmonics (Sinc function), SSA a priori cannot: it is difficult to create a harmonic of a signal which is not periodical stricto sensu. See Note II.

15:
Actually, the opposite would be much more unexpected. This is even more true for the insolation itself !… Why for insolation, that is calculated from solar irradiance (a constant), eccentricity, obliquity and precession that are all cyclical?
16: Not finding the main and secondary astronomical periodicities that are used as input of the insolation computation would be a sure signature of some methodological error. I therefore do not understand the value of applying a spectral analysis method (in particular on insolation), without discussing the "detection level" problem.
The secondary astronomical periodicities are not input but actual observations! We fail to see the link with the previous comment. As already stated above, the detection level is a concept that is computer-dependant (32 vs 64 bit, Linux vs Windows, MKL vs Matlab, etc … ) and not method-24/28 dependant. An example is given in the link below, that shows how the actual realization of a decomposition depends on the computer, the library, the software and the nature of the processor.
https://www.boost.org/doc/libs/1_47_0/libs/math/doc/sf_and_dist/html/math_toolkit/special/sf _poly/sph_harm.html 17: Indeed, the key mathematical difficulty is to decide whether it is "noise" or "signal" using some statistical test. There is no mention of any statistics in this paper, which certainly disqualifies it entirely. This is strange since standard SSA tools used in the climate community (eg. SSA-MTM toolkit) do provide some statistical tools and safeguards.
Again on statistics. See above and Note II. As soon as one exits from the (restricted) climate community, people use it in the same way we do (see references above). Statistics is not what qualifies a signal as being noise or not. By definition, a noise is a signal (stochastic or not) that is superimposed on a signal of interest. For instance, when one performs electrical measurements in urban zones, one often adds a 50 Hz signal, that is an almost perfect sine function. Which statistics will tell us that the 50 Hz signal is a noise?

18:
The harmonics and combination tones found by the authors have often a very small amplitude (below 1% or 0.1% of the variance) and to prove their significance is certainly not an easy task, given that the records are rather "short" and "noisy". This is precisely the main difficulty of spectral analysis, and a problem that climate scientists are usually facing. But the authors of this paper are not even mentioning or discussing this central point, as if they were unaware of the main problem of "spectral analysis estimation".
These small amplitudes must be discussed in parallel with the trend (see point 9). Which as we argue could be a "short" segment of a component with a quasi-period much longer than the interval of time over which the time series is known. And if we let aside the trend, then most of the components (quasi-periods) fit the characteristic astronomical quasi-periods. This analysis is not necessarily "difficult"; see the critical sampling intervals in Papoulis (1977) or Claerbout (1984).
Moreover, we have already answered this point (Note II).

19:
In addition, it is worth mentioning that detecting several significant frequencies (let alone NUMEROUS significant frequencies) is a much more difficult problem than detecting a single one.
Point already made, already answered. This might have been problematic in the '60s, but not 25/28 any more. One can now reconstruct a hushed voice with a mobile phone (using piezo electricity thus frequencies).

20:
Indeed, the null-hypothesis cannot be a simple (red or white) noise hypothesis.
What is the point?

21:
A classical work-around is to perform Monte-Carlo simulations of synthetic signals with and without each single periodicity and assess whether the results are statistically distinguishable.
But this would probably be a very different paper.
If there are synthetic signals, then there is a model behind. In the present paper, we do not introduce any model but we perform a signal analysis of cores that have already been analyzed in several papers and in which NO model is necessary.
We do not use Monte-Carlo since the work by Kirkpatrick (1977) on Simulated annealing. It is much more efficient (and elegant) than a simple distribution of models. In signal analysis, statistics underlie the matrices and decompositions we use (e.g. Mañé, 1981 ;Takens, 1981 ;or Lemmerling and Van Huffel, 2001 components with dubious relevance. The same phenomenon applies to all spectral methods: the more degrees of freedom you are allowing, the more frequencies you get, but the less significant they are. What is the embedding dimension of the analysis here? This is not even mentioned ! But according to line 62, "SSA was able to extract up to 75 cycles". My guess is that the embedding dimension is a few hundreds, which represents quite a lot of freedom for the method to generate a 26/28 lot of "spurious" frequencies. We agree with the first part of this comment. But in our case the embedding dimensions (given the size of the time windows over which the series are defined) are small! In iterative SSA (iSSA) they change at each step. Should we add a glacial cycles while a more recent one (Dome C), covering twice as many cycles, has been available for more than a decade now (see eg: Guo et al (2012) for the SSA analysis of CH4 at Dome C over 8 cycles, not 4 as here).
It is not strange: it is one of the most quoted sequences. We wanted to check our method and thought the reader might be interested by this comparison. Of course (if decided by the editor), we could do the analysis for the Guo et al (2012) series. But in our view, it is not needed in the frame of the present paper.

24:
Obviously, the results are likely to be statistically more robust with a longer time series, if statistics matters of course. "Finding not only the main expected Milankovitch periodicities but also many "secondary" components with much smaller amplitudes gives confidence in our iterative SSA method". But I would expect standard methods to do the same, with attributing to these periodicities a very low confidence. In any case again, finding many frequencies is easy and NOT the main problem of "spectral analysis". Besides, what is the subject of the paper? Is it about climate or mathematics?
What does this mean? Why separate both? We use mathematical techniques to test and check their power in analyzing actual observational time series commonly associated with climate and its changes. Otherwise, see Note II.

25:
The authors appear unaware that the word "insolation" is not sufficient to describe the astronomical forcing of climate. They are presenting the spectral analysis of "the insolation of Laskar et al (2004)". This is clearly a very insufficient description. When looking more closely at the spectrum (Figure 3), it seems that they are in fact discussing global averaged insolation (linked 27/28 to eccentricity only)! Otherwise the spectrum would be dominated by the 41k and the 23k periodicities. Since there is no power in the obliquity and precession band, the only explanation is that the authors are using global averaged insolation. But this is a major misunderstanding on the Milankovitch forcing. Indeed, it is well known since the 19 th century that global insolation changes are much too weak to affect climate and scientists have discussed local seasonal insolation, not mean global one. When climate scientists talk about "insolation" in the context of glacial cycles, they are talking about summer insolation at high northern latitude according to Milankovitch theory.
THIS IS NOT what is used in this analysis.
See Note II.

26:
In other words, this paper has missed entirely its claimed goal of discussing the Milankovitch forcing, simply by looking at the wrong insolation… To conclude, its seems to me that the authors have misunderstood both the aim of "spectral analysis" and the foundations of the "Milankovitch theory". I cannot be favorable to publication.
On these two points, we have shown that the reviewer is wrong. We have specifically devoted the two longer Notes I and II above to explain where and how. We appeal to the Editor and confirm our submission. 28/28