BS09 and Mannian Smoothing

Many CA readers have probably noticed Lucia’s recent coverage of Benestad Schmidt (JGR 2009) (with the inevitable acronym BS09). BS09 was an effort by B&S to verify Scafetta and West.

realclimatescientists can also be climateauditors, I guess. It’s too bad that they don’t spend the same amount of energy examining Mann et al 2008 or something like that.

Scafetta replied to BS09 at Pielke Sr here. The rebuttal point that interested Lucia was an issue relating to end points. Scafetta argued that BS09 had failed to understand wavelet decomposition and urged the BS authors to study a wavelet primer. While I’ve used wavelet decomposition in R from time to time and have written algorithms using modwt, my initial reaction was that once they started arguing about wavelets, you’d have to see their scripts to fully understand what the dispute was about.

In fact, the matter is much simpler than that and pertains only to end-point padding, an issue that we’ve covered here on a number of occasions, especially in connection with Mannian smoothing and Rahmstorf smoothing. Lucia referred to the latter discussion, but the discussion has been going longer than that.

As we’ve discussed, there are several methods in play for endpoint padding to create a smooth: padding with the mean (IPCC); reflection; reflection-and-flipping (Mann); linear extrapolation (Rahmstorf).

An option embedded in some algorithms – one which is highly inappropriate for trending series – is “cyclic” padding. Where you pad the end with values from the beginning.

Here is BS09 Figure 4. Note in particular the downslope of the smoothed version of the various series at the end.

Figure 1. Excerpt from BS09.

In his rebuttal at Pielke Sr, Scafetta observes that the BS09 smooth was derived by cyclic padding. . Scafetta observed:

By applying a cyclical periodic mode Benestad and Schmidt are artificially introducing two large and opposite discontinuities in the records in 1900 and 2000, as the above figure shows in 2000. These large and artificial discontinuities at the two extremes of the time sequence disrupt completely the decomposition and force the algorithm to produce very large cycles in proximity of the two borders, as it is clear in their figure 4. This severe error is responsible for the fact that Benestad and Schmidt find unrealistic values for Z22y and Z11y that significantly differ from ours by a factor of three

Here is the figure illustrating the point:

Figure 2. Figure from Pielke Sr blog illustrating the cyclic padding in BS09.

This particular issue is unrelated to wavelet decomposition in the sense that cyclic padding with a gaussian filter or something like that would have yielded a more or less similar smooth. In the latter case, however, you’d have had to make a decision on endpoint padding and even the BS09 authors would have been unlikely to include cyclic padding in their written description. It seems to have crept in as a default mechanism in the algorithm, to the authors’ bad luck.

In a comment at Lucia’s, Gavin blames Scafetta and West for failing to describe their methodology clearly – a criticism that seems a little cheeky for a realclimate colleague of Mann and Steig:

B&S09 clearly stated that we had not been able to fully emulate Scafetta and West’s methodology and so statements that we did something different to them were to be expected. The issue with the periodic vs. reflection boundary conditions in the wavelet decomposition does make a difference – but what they used was never stated in any of their papers (look it up!).

As I noted at Lucia’s, I was unable to find a “clearly statement that we had not been able to fully emulate Scafetta and West’s methodology and so statements that we did something different to them were to be expected”. On the contrary, BS09 seems to clearly say that they had been able to repeat their analysis. BS09:

We repeated the analysis in SW06a and SW06b, and tested the sensitivity of the conclusions to a number of arbitrary choices. The methods of SW06a and SW06b were used (1) …. We reproduced the SW06a study…”.

There is no mention of any concern over how to emulate BS09 endpoint padding. Limited concerns were expressed about getting “exactly the same ratio of amplitudes”, to a “slight” difference in lagged values and to the ACRIM discussion in a separate paper by Scafetta and Willson. Perhaps there’s a “clear statement” elsewhere in the article that I missed. If so, perhaps someone can bring it to my attention.

In any event, the use of cyclic padding is a clanger. I presume that Scafetta is going to pursue this as a reply, in which case, I suggest that it would be useful for him to compare BS09 endpoint padding to IPCC padding, Mannian padding and Rahmstorf padding.

It looks like the realclimatescientists are going to have to start on another Corrigendum. Predictably Gavin declared that the error doesn’t “matter”. Funny, Team errors never seem to “matter”.

Remarkably, Gavin, using terms that I might have used, commented:

Note we still don’t have a perfect emulation, so perhaps you guys could agitate for some ‘code freeing’ to help out. 🙂

Gavin’s got this part exactly right. Verbal descriptions of methods seldom provide enough information; code saves a lot of time and energy. This has obviously been a CA campaign for a long time. I, for one, am quite willing to do my part in “freeing” the code; I’ve already emailed Nicola Scafetta. (In this case, it seems logical to also provide the code for BS09.) Maybe Gavin will join us in resolving long outstanding questions pertaining to MBH that have frustrated us for years (and ignored by Wahl and Ammann): how were the confidence intervals in MBH99 calculated? how were the number of retained principal components calculated?

So Gavin, for MBH,

Note we still don’t have a perfect emulation, so perhaps you guys could agitate for some ‘code freeing’ to help out. 🙂

26 Comments

  1. Posted Aug 7, 2009 at 12:07 PM | Permalink

    I would probably prefer the padding by setting “T” at all future points equal to the last known “T”. At any rate, if some conclusions don’t survive the change of the padding to some of the other algorithms, I would say that they’re not robust. It’s the old observation by Feynman that one can’t build his conclusions on the last point in a graph. If the people were sure that it was reliable, they could add one more. In some sense, it applies to time series, too.

    Is the usage of wavelets and their chosen shape justified? Many of these things look arbitrary to me. They can be OK to get some idea or conclusion if it is robust but if it is not, it is not hard to see that people write papers like BS09.

    • John S.
      Posted Aug 7, 2009 at 1:22 PM | Permalink

      Re: Luboš Motl (#1),

      the old observation by Feynman that one can’t build his conclusions on the last point in a graph

      is certainly a propos here. Wavelet analysis, however, does not rely upon that last point, per se, in looking for nonstationarities, which are not the real issue in Scafetta’s analysis anyway. Although Scafetta could have made his original point perhaps more strongly via other techniques, it is BS09 that makes the bonehead end-point error, while indulging in gratuitous talk about robustness.

      There is nothing “non-robust” about properly applied wavelet analysis despite the different forms of wavelets that may be employed. The onus is entirely on the analyst to make the proper choice (similar to that of choosing a spectral window) and run the analysis intelligently. Those choices affect only the granularity and not the substance of the picture obtained.

      • Steve McIntyre
        Posted Aug 7, 2009 at 1:37 PM | Permalink

        Re: John S. (#10),

        IMO, the issue here has very little to do with wavelet analysis to the extent that it is used as a form of smoothing – and I really don’t want to spend time on the wavelet aspect of this.

        Let’s stick to the endpoint padding issue.

        • John S.
          Posted Aug 7, 2009 at 1:48 PM | Permalink

          Re: Steve McIntyre (#11),

          Dare I point out that it Lubos Motl (#1) who raised the wavelet aspect in the first place? I certainly don’t wish to waste anyone’s time.

        • Mark T
          Posted Aug 7, 2009 at 2:11 PM | Permalink

          Re: Steve McIntyre (#11), Uh, yeah, that’s what I was trying to say… Ahem. 🙂

          Mark

  2. Steve McIntyre
    Posted Aug 7, 2009 at 12:12 PM | Permalink

    Another aspect of this paper is the problem of collinearity between solar and CO2 in a regression. MBH98 did a linear regression of the type criticized. Chefen and I did some useful posts on this a couple of years ago.

  3. Posted Aug 7, 2009 at 12:22 PM | Permalink

    A bit of a ‘clanger’. I like it.

    I should mention that it doesn’t appear from Gavin and Scafetas comments that Gavin or his team have been in contact.

    From Dr. Scafetta on Lucia’s.

    The problem I saw in BS09 is not just the evident mathematical errors it contains but the arrogance of the tone. They never contacted me to ask a clarification, for example.

    But this clanger thing is also clearly stated:

    It is the way how they used the code that is wrong because they do not know how wavelet works and how they must be used to let them to work.

    It was big of Dr. Schmidt to admit error, but the end point treatment they used made so little sense it’s hard to explain.

  4. John S.
    Posted Aug 7, 2009 at 12:35 PM | Permalink

    Dare I say that BS09 is robustly wrong?

  5. Tim
    Posted Aug 7, 2009 at 12:44 PM | Permalink

    In RC’s weekly roundup they mention the recent paper on ENSO and global temps. They end the paragraph with the line, “Nevermore let it be said that you can’t get any old rubbish published in a peer-reviewed journal!”

    I don’t suppose they see the irony in their statement 🙂

  6. stan
    Posted Aug 7, 2009 at 12:45 PM | Permalink

    Perhaps someone should introduce Gavin’s goose to the gander. [Before it’s cooked.]

    Or advise him of the admonition about living in glass houses.

  7. Posted Aug 7, 2009 at 1:00 PM | Permalink

    Note we still don’t have a perfect emulation, so perhaps you guys could agitate for some ‘code freeing’ to help out. 🙂

    Wow. Just wow.

    And wow again…

  8. steven mosher
    Posted Aug 7, 2009 at 1:18 PM | Permalink

    I guess you won’t be handing out brownies. hehe.

  9. Mark T
    Posted Aug 7, 2009 at 2:10 PM | Permalink

    I should point out that the end-point issue with doing a circular transform is not unique to wavelets. Wavelets as implemented here are nothing more than a class of FIR filters that obey certain properties (which are useful for perfect reconstruction). The normal implementation of a FIR filter is a convolution. For some applications, it makes sense to implement the convolution circularly, i.e., by tying the end to the beginning and performing the convolution in a circle. Doing so inevitably forces discontinuities between the end points, be they wavelets or other functions serving as the basis.
    .
    Re: Luboš Motl (#1),

    Is the usage of wavelets and their chosen shape justified?

    I suppose it depends upon your goal. Without a priori knowledge of the “signal” you’re looking for, I’m not sure any one choice is better than another as long as the choice obeys the specific requirements for being called a wavelet.
    .

    Mark

  10. bernie
    Posted Aug 7, 2009 at 2:26 PM | Permalink

    If he needed the code Gavin should have asked:

    They never contacted me to ask a clarification… (Dr. Scafetta on Lucia’s blog)

    BS09 reads like an example of “Ready, Fire, Aim”, which always increases the chances that you will shoot yourself in the foot.

  11. Andrew
    Posted Aug 7, 2009 at 3:09 PM | Permalink

    I love irony! Hehe…

  12. Steve McIntyre
    Posted Aug 7, 2009 at 3:20 PM | Permalink

    I sent the following email to Dr Gavin Schmidt of NASA GISS:

    Dear Dr Schmidt

    In a recent comment at Lucia’s blog, you stated:

    “Note we still don’t have a perfect emulation, so perhaps you guys could agitate for some ‘code freeing’ to help out. ”

    As you may be aware, I’ve been trying for some time to achieve a more “perfect emulation” of MBH methodology. While some code was provided in connection with the House Energy and Commerce Committee hearings a few years ago, the code was unfortunately incomplete. The code did not include the steps in which the confidence intervals in MBH99 were calculated nor how the number of retained principal components for the various tree ring networks was calculated. Existing information is insufficient to permit either to be emulated. Aside from myself, some very able Climate Audit readers (Jean S, UC) have tried hard to figure out these steps and have been frustrated.

    I would appreciate it if you would attempt to persuade your associates to provide the relevant code for these steps. I will do what I can to persuade Scafetta to provide a script for their calculations.

    Resolving small measures like this can go a long way.

    Regards,

    Steve McIntyre

  13. Robinedwards
    Posted Aug 7, 2009 at 4:42 PM | Permalink

    Unfortunately my carefully composed contribution seems to have vanished before being submitted :-(( I’ll try to start again.

    I fear that I must have missed a vital point early in these discussions of smoothing, endpoint treatment and wavelets. I admit immediately that I know /nothing/ about wavelet analysis, so will have to assume that in common with other methods of smoothing it is intended to produce believable (plausible?) estimates of as yet ungathered observations. In other words, it is hoped that the future of the series can be predicted. This demands (I think) that a model must be hypothesised and then tested against known observations so that one can be reasonably convinced that model and observations are mutually “consistent”. Having achieved this the model is projected into the future and its consequences weighed in the balance of politics and science.

    What I would really like to know is the provenance of the data that were used to generate the intriguing plots provided by Steve. The cyclic component is truly amazing! Where can they be found, please? To my untutored eye some of the curves appear to be creations of fantasy – though since I cannot see the individual data points I may be fantasising too.

    Is anyone really going to place credence in these methods? I would like to see the result of applying them to some well-known climate data. I would suggest trying observations from Anchorage, Fairbanks and other Alaskan sites for the thirty year period 1946 to 1976, totally ignoring data from the following years, applying the endpoint manipulation methods of the types described by Steve, and using the fitted model to forecast events from mid 1976 onward. My expectation is that a simple model will do a brilliant job on the known data, 1946 to 1976, producing forecasts for 1977 to say 1987 with impressively narrow confidence bands. This exercise should be repeated for the 30 years after 1976, using the data that emerged during that period. Comparison of the two fitted models and their relationships to the actual data numbers should provide some talking points.

    Underlying all this is my feeling that climate for a given location or region may be essentially resistant to any “rational” attempts to forecast it. A “random” component may in practice outweigh the best-laid schemes of skilled operators in the art of climate forecasting.

    Robin

  14. Posted Aug 7, 2009 at 5:56 PM | Permalink

    Hi –

    To this econemetrician, this really looks like someone desperately thrashing the data. There is no connection between the original data and the final results.

    Of course, that may be the point…

  15. sky
    Posted Aug 8, 2009 at 2:29 PM | Permalink

    Scafetta and West’s conclusions do not stand or fall on the basis of any end-point issue (although it’s not clear how they resolved it)*. It’s the efficacy of wavelet analysis as a bandpass filter throughout the entire record that matters in revealing the similarities between the solar and temperature records. By contrast, SB09’s critique does fall on the basis of inept cyclical padding of trending data. I think that’s the point that John S was driving at in this discussion, which got off on the wrong foot from the start.

    Uncritical resort to cyclical padding has become fairly common practice not only in filtering but in spectral estimation, where the entire record is assumed to repeat indefinetly in FFT-based analyses. When the record is very long and suitable decimation techniques are used, the results can be quite indicative despite such assumptions. That luxury is very seldom available in climate studies. My employer wisely bans all analyses that involve data padding.

    * Maybe Mark T can shed some expert light on this question.

    • Craig Loehle
      Posted Aug 8, 2009 at 6:22 PM | Permalink

      Re: sky (#20), I, personally, have banned myself from using any padding. And I am not nearly clever enough to loop my data in a circle…

      • sky
        Posted Aug 9, 2009 at 5:43 PM | Permalink

        Re: Craig Loehle (#21),

        Good for you. Cleverness is not even required, it’s built in to the periodic nature of Fourier series analysis.

    • Nick Stokes
      Posted Aug 9, 2009 at 5:46 PM | Permalink

      Re: sky (#20), I believe I was initially misled by Scafetta’s framing of the cyclic padding issue, and I think others have too. B&S, as I now see it, did not use cyclic padding for trended data. They fitted a fifth order polynomial to get their initial smooth, reflecting the trend. This does not involve any padding. They then subtracted this to approximate a stationary process, which they band filtered by wavelet analysis, recovering two bands D8 (centred on a 22yr cycle) and D7 (11yr). It was in this band filtering that cyclic padding was used. BS say that S&W used reflection padding.

      Cyclic padding, of course, introduces a discontinuity – reflection makes a discontinuous derivative. Both produce Gibbs effects. In this case, they are not so large, because the process (D8 and D7) is at least approximately stationary, and just a part of the total signal.

      • sky
        Posted Aug 10, 2009 at 4:33 PM | Permalink

        Re: Nick Stokes (#23),

        I don’t think your explanations really fit.

        Had BS09 properly detrended the data, the downturn at the end of their bandpass shouldn’t be as severe as it is. Furthermore, Gavin has effectively admitted to error in using cyclical padding. Lacking any code from him, nobody really knows what was actually done.

        The higher-frequency ripple known as the Gibbs effect is manifest only between cardinal points of a discrete series when insufficient terms are retained in the Fourier series expansion of a sharp-featured signal. This is not an issue here. Both wavelet and F. analysis provide exact results at the cardinal points. If anything, I would expect customarily very smooth finite-duration wavelets to smooth over discontinuities.

        • Nick Stokes
          Posted Aug 11, 2009 at 1:04 AM | Permalink

          Re: sky (#24), First I don’t think Gavin has admitted an error using cyclic padding. He shouldn’t have, because I think his approach using a low-order polynomial to detrend is better than what S&W did, and much reduces the dependence on padding anyway. I think what he did concede was that what B&S did was not a complete emulation of what S&W did.

          I didn’t understand your argument about cardinal points, but it certainly isn’t true that Gibbs oscillations are only manifest when insufficient Fourier terms are retained. Sharp featured signals have many harmonics, extending into high frequencies. They only remain sharp featured if all these are retained with amplitude and phase unchanged. Any filtering changes this, and the harmonics then become evident as oscillations within the data set.

          In the case of these bandpass filters, the artefacts created by the boundary treatment generate oscillations within the frequency band. There would have been higher frequencies as well, but they are attenuated by the band pass.

        • sky
          Posted Aug 11, 2009 at 5:43 PM | Permalink

          Re: Nick Stokes (#25),

          There’s really no solid ground for speculation when only the authors know (???) what they did with the data.

          No matter how sharp-featured or discontinuous the signal, if it satisfies the Dirichlet conditions, the entire stretch of record can be reproduced over continuous time with an infinite Fourier series. That reproduction, however, is not extendable to other stretches, unless the signal itself is strictly periodic. When that signal is sampled at discrete time-intervals, the F. series is necessarily truncated to a finite series at the Nyquist frequency. If there’s no aliasing introduced by sampling, that reproduceability survives. Otherwise, truncation of the series produces the Gibbs effect between the cardinal points at which the underlying signal was sampled. The exactitude of reconstruction at the cardinal points survives in any event. Any post-sampling filtering of the discrete data series can only operate on the frequencies below or at Nyquist. With genuine band-pass filtering, there’s no need to detrend the data; the filter removes that trend by itself.

          I really can’t take any more time to clear up your persistent confusions about the basics of signal analysis/synthesis.