Checking a wavelet analysis

Three ways a wavelet scalogram misleads you: a scale bias in the global spectrum, an event that impersonates a cycle, and the resolution knob you never set.

The previous three posts built the machinery: the transform, the significance test, the coherence. Each of those had its own trap and its own fix. This one is about the failures that survive all of them: the map is correct, the null is right, the patch is areawise significant, and the reading is still wrong.

Three of them, each with a number attached.

Check one: the global spectrum favours long periods

The global wavelet spectrum is the time-average of the power map, the wavelet answer to “which period dominates the record”. Feed it two cycles of exactly equal amplitude, one at period 4 and one at period 32, and ask which is bigger.

n <- 1024
t <- 1:n
x <- sin(2 * pi * t / 4) + sin(2 * pi * t / 32)     # identical amplitudes
o <- cwt(x)
P <- Mod(o$W)^2
gws <- rowMeans(P)
i4 <- which.min(abs(o$period - 4)); i32 <- which.min(abs(o$period - 32))
c(gws_at_period_4 = round(gws[i4], 1), gws_at_period_32 = round(gws[i32], 1),
  conventional_ratio = round(gws[i32] / gws[i4], 2), true_ratio = 1)

   gws_at_period_4   gws_at_period_32 conventional_ratio         true_ratio 
              3.30              25.20               7.74               1.00 

The period-32 cycle looks 7.7 times more powerful than the period-4 cycle. They have the same amplitude. The bias is not subtle and it is not noise: it is in the normalisation. The Morlet transform as conventionally normalised returns power proportional to the scale, so any energy at a long period is magnified simply for being at a long period.

Liu et al. (2007) traced it and gave the fix: divide the power by the scale before averaging.

gws_rect <- rowMeans(P / o$scale)
c(rectified_at_period_4 = round(gws_rect[i4], 3),
  rectified_at_period_32 = round(gws_rect[i32], 3),
  rectified_ratio = round(gws_rect[i32] / gws_rect[i4], 2), true_ratio = 1)

 rectified_at_period_4 rectified_at_period_32        rectified_ratio 
                 0.814                  0.787                  0.970 
            true_ratio 
                 1.000 

0.97 against a truth of 1. That is the same division by s that appears inside the coherence formula, and it is needed for the same reason.

Figure 1: Global wavelet spectrum of a signal containing two equal-amplitude cycles, before and after dividing power by scale. Each curve is scaled to its own maximum so the shapes can be compared. The conventional version makes the period-32 cycle look several times larger than the period-4 cycle.

This one matters because the global spectrum is where people read off “the” dominant period, and because the bias always points the same way. Long cycles are the interesting ones. The method hands them to you.

Check two: an event is not a cycle

Here is a series with no cycle in it at all: white noise, plus a single spike at year 128. One year, one shock, nothing periodic.

set.seed(287)
n2 <- 256; t2 <- 1:n2
z <- rnorm(n2)
z[128] <- z[128] + 12                                # one bad winter
oz <- cwt(z)
Pz <- Mod(oz$W)^2
okz <- outer(oz$period, oz$coi, "<=")
az <- as.numeric(acf(z, lag.max = 1, plot = FALSE)$acf[2])
flag <- (Pz > sig95(oz, z, az)) & okz
per_flagged <- range(oz$period[apply(flag, 1, any)])
c(share_of_map_flagged = round(mean(flag[okz]), 3),
  flagged_from_period = round(per_flagged[1], 2),
  flagged_to_period = round(per_flagged[2], 2))

share_of_map_flagged  flagged_from_period    flagged_to_period 
               0.055                2.070               20.350 

The spike lights up 5.5% of the map, across periods from 2.1 to 20.3. Against a correctly specified red-noise null. There is no cycle here, and the map says there is cycling at every timescale it can see.

Worse, look at what the flagged region says about duration:

Worse, the plume has width in time as well. The wavelet at period 8 is several years wide, so a one-year shock comes back looking like a short burst of oscillation rather than a shock:

b8 <- which.min(abs(oz$period - 8))
yrs <- which(flag[b8, ])
col_periods <- oz$period[flag[, 128]]
c(at_period = round(oz$period[b8], 2),
  flagged_from_year = min(yrs), flagged_to_year = max(yrs),
  column_period_span = round(max(col_periods) / min(col_periods), 1))

         at_period  flagged_from_year    flagged_to_year column_period_span 
              8.26             122.00             134.00               9.20 

The last number is the tell. A vertical slice through the plume, at the year of the spike, is significant over a 9.2-fold range of periods. Compare a real transient cycle:

set.seed(3287)
env <- ifelse(t2 >= 100 & t2 <= 156, 1, 0)
r <- rnorm(n2) + 3 * env * sin(2 * pi * t2 / 8)      # a genuine period-8 transient
orr <- cwt(r)
Pr <- Mod(orr$W)^2
okr <- outer(orr$period, orr$coi, "<=")
ar <- as.numeric(acf(r, lag.max = 1, plot = FALSE)$acf[2])
flag_r <- (Pr > sig95(orr, r, ar)) & okr
col_r <- orr$period[flag_r[, 128]]
c(event_column_span = round(max(col_periods) / min(col_periods), 1),
  cycle_column_span = round(max(col_r) / min(col_r), 1))

event_column_span cycle_column_span 
              9.2               1.6 

A vertical slice through the event covers a 9.2-fold span of periods; a slice through the real cycle covers 1.6-fold. That is the diagnostic, and it costs one line of code: a cycle is narrow in period and long in time; an event is narrow in time and broad in period. Look at the shape of the patch, not its colour.

Figure 2: One isolated spike in white noise (upper) against a genuine period-8 transient (lower), both tested against a correctly fitted red noise null. The spike produces a broad vertical plume; the cycle produces a narrow horizontal bar.

Ecological records are full of single events: a hard winter, a fire, a cull, a disease introduction. Every one of them will produce a plume, and the plume will be significant.

Check three: the knob you never set

Every function above takes w0 = 6 and nobody ever changes it. It is the number of oscillations inside the Morlet wavelet’s envelope, and it sets the trade you make between resolution in time and resolution in period. There is no default that is correct; there is only a default that is conventional.

Make it visible on a signal with a known answer: period 10 for the first half, period 5 for the second, no noise at all.

half <- n2 / 2
sig2 <- c(sin(2 * pi * (1:half) / 10), sin(2 * pi * ((half + 1):n2) / 5))
fwhm <- function(g, i) {                             # width of the peak at i, at half its height
  half_h <- g[i] / 2
  lo <- i; while (lo > 1 && g[lo - 1] >= half_h) lo <- lo - 1
  hi <- i; while (hi < length(g) && g[hi + 1] >= half_h) hi <- hi + 1
  c(lo, hi)
}
probe <- function(w0) {
  ow <- cwt(sig2, w0 = w0)
  Pw <- Mod(ow$W)^2
  g <- rowMeans(Pw / ow$scale)                       # rectified, per check one
  b10 <- which.min(abs(ow$period - 10))
  w <- fwhm(g, b10)
  band <- Pw[b10, ] / max(Pw[b10, ])                 # how long to notice the cycle stopped?
  hi <- max(which(band[1:half] > 0.9))
  lo <- min(which(band[half:n2] < 0.1)) + half - 1
  c(w0 = w0,
    fourier_factor = round(fourier_factor(w0), 4),
    period_halfwidth = round(ow$period[w[2]] - ow$period[w[1]], 2),
    transition_years = lo - hi,
    coi_cost_at_period_10 = round(10 / (fourier_factor(w0) / sqrt(2)), 1))
}
tab <- t(sapply(c(4, 6, 12), probe))
tab

     w0 fourier_factor period_halfwidth transition_years coi_cost_at_period_10
[1,]  4         1.5246             3.33               14                   9.3
[2,]  6         1.0330             2.05               21                  13.7
[3,] 12         0.5218             0.69               41                  27.1

Read the columns. period_halfwidth is how blurred the period-10 peak is, and w0 = 12 resolves it 4.8 times more sharply than w0 = 4. transition_years is how long the transform takes to notice the period-10 cycle stopped, and there w0 = 12 needs 2.9 times longer. coi_cost_at_period_10 is how many years at each end of the record you lose at that band, and w0 = 12 costs you 3 times as many.

That is the trade-off, made numeric. A sharper answer about which period costs you a vaguer answer about when, and it costs you record length as well. The convention w0 = 6 is a compromise, not a discovery, and if your question is “did the period change in 1987” rather than “is the period 9 or 11” then the compromise is not the one you want.

Figure 3: The period-10 band through the same regime shift, at three settings of the Morlet parameter. All three see the cycle stop halfway through the record; they disagree about how sharply.

The checklist

Everything the four posts have established, in the order you would actually use it:

1. Is the null red rather than white? If it is white, the map is decoration. (post two)
2. Is the significant patch bigger than the patches the null makes on its own? Pointwise significance says nothing. (post two)
3. Is the patch inside the cone? Do not read the first or last 0.73 * period years of any band.
4. Is the patch narrow in period and long in time, or is it a plume from an event?
5. If you read a dominant period off the global spectrum, did you divide power by scale first?
6. If you report a phase lag, is it the only lag consistent with that angle, or one of several? (post three)
7. Could a common driver produce the same coherence? It almost always could. (post three)
8. Did you choose w0, or did it choose you?

Honest limits

None of these checks makes wavelet analysis confirmatory. The method is descriptive: it says where in the time-period plane the variance sits, and the tests say only that a particular background model does not account for it. A significant, narrow, well-placed, correctly rectified patch is a reason to go and look at what happened in those years. It is not a cycle, and the map cannot tell you what caused it.

Related tutorials

• Wavelet analysis of population cycles
• Wavelet significance and the red noise null
• Wavelet coherence and phase
• Spectral analysis of population cycles

References

Torrence and Compo 1998 Bulletin of the American Meteorological Society 79(1):61-78

Liu, Liang and Weisberg 2007 Journal of Atmospheric and Oceanic Technology 24(12):2093-2102 (10.1175/2007JTECHO511.1)

Maraun, Kurths and Holschneider 2007 Physical Review E 75:016707 (10.1103/PhysRevE.75.016707)

Cazelles, Chavez, Berteaux, Menard, Vik, Jenouvrier and Stenseth 2008 Oecologia 156(2):287-304 (10.1007/s00442-008-0993-2)

Percival and Walden 2000 Wavelet Methods for Time Series Analysis, Cambridge University Press, ISBN 978-0-521-64068-8