What is this all about?

This is about (irregular, “locally stationary”) time series that can be thought as models for what is seen in either (frequency-time smoothed) spectrogram or (scale-time smoothed) scalogram.

From a different perspective, if a smoothed raw spectrum is associated with stationary series, what are the theoretical time series models for smoothed spectrogram or scalogram?

Motivating examples

Tree rings

# From the R package LSTS (Locally stationary time series)

library(LSTS)
library(rdatamarket)

malleco = dmlist("22tn")
#mammothcreek = dmlist("22r7")
y = malleco$Value
#y = mammothcreek$Value

plot(y,type="l")  

block.smooth.periodogram(y, spar.freq = 0.8, spar.time = 0.8, theta = 45, phi = 0)

library(Rwave)

isize <- length(y)
nvoice <- 25
freqstep <- 1/25
scale <- 20
cgty <-cgt(y,nvoice,freqstep,scale,plot=FALSE)
image(1:isize, seq(0, nvoice * freqstep/2, length = nvoice), Mod(cgty), xlab = "Time", ylab = "Frequency",col = grey(seq(0, 1, length = 256)))
        title("Gabor Transform Modulus")

 # Analysis by blocks of phi and sigma parameters
T. = length(y)
N = 200
S = 100
M = trunc((T. - N)/S + 1)
table = c()
for(j in 1:M){
    x = y[(1 + S*(j-1)):(N + S*(j-1))]
    table = rbind(table,nlminb(start = c(0.65, 0.15), N = N, objective = LS.whittle.loglik,series = x, order = c(p = 1, q = 0))$par)
    }
u = (N/2 + S*(1:M-1))/T.
table = as.data.frame(cbind(u, table))
colnames(table) = c("u", "phi", "sigma")
par(mfrow = c(1,2), cex = 0.8)
spar = 0.6
plot(smooth.spline(phi, spar = spar)$y ~ u, data = table, pch = 20, ylim = c(0,1),
         xlim = c(0,1), las = 1, bty = "n", ylab = expression(phi(u)))

plot(smooth.spline(sigma, spar = spar)$y ~ u, data = table, pch = 20, ylim = c(0,0.2),
       xlim = c(0,1), las = 1, bty = "n", ylab = expression(sigma(u)))

tsl <- 500
phit <- seq(0,.9,.9/tsl)

ts <- rep(0,tsl)

set.seed(1)
ts[1] <- rnorm(1)
for (k in 2:tsl){
  ts[k] <- phit[k]*ts[k-1] + rnorm(1)
}

plot(ts,type="l")

block.smooth.periodogram(ts, spar.freq = 0.8, spar.time = 0.8, theta = 45, phi = 0)

isize <- length(ts)
nvoice <- 25
freqstep <- .04
scale <- 15
cgtts <-cgt(ts,nvoice,freqstep,scale,plot=FALSE)
image(1:isize, seq(0, nvoice * freqstep/2, length = nvoice), Mod(cgtts), xlab = "Time", ylab = "Frequency",col = grey(seq(0, 1, length = 256)))
        title("Gabor Transform Modulus")

library(strucchange)
## Loading required package: sandwich
y0 <- y[2:length(y)]
yl <- y[1:length(y)-1]
y0yl <- data.frame(cbind(y0,yl))
colnames(y0yl) <- c("ts","tsl")


ty.model <- y0 ~ yl

ocusy <- efp(ty.model, type="OLS-CUSUM",data=y0yl)
mey <- efp(ty.model, type="fluctuation", data=y0yl, h=0.2)

bound.ocusy <- boundary(ocusy, alpha=0.05)
plot(ocusy, boundary = FALSE)
lines(bound.ocusy, col = 4)
lines(-bound.ocusy, col = 4)

plot(mey, functional = NULL)

ts0 <- ts[2:length(ts)]
tsl <- ts[1:length(ts)-1]
ts0tsl <- data.frame(cbind(ts0,tsl))
colnames(ts0tsl) <- c("ts","tsl")


t.model <- ts0 ~ tsl

ocus <- efp(t.model, type="OLS-CUSUM",data=ts0tsl)
me <- efp(t.model, type="fluctuation", data=ts0tsl, h=0.2)

bound.ocus <- boundary(ocus, alpha=0.05)
plot(ocus, boundary = FALSE)
lines(bound.ocus, col = 4)
lines(-bound.ocus, col = 4)

plot(me, functional = NULL)

Baby ECG data

# From Nason (2008), Chapter 5
# 1=quiet sleep, 2=between quiet and active sleep, 3=active sleep, 4=awake.

library(wavethresh)
data("BabyECG")
data("BabySS")

myhrs <- c(22, 23, 24, 25, 26, 27, 28, 29, 30)
mylab <- c("22", "23", "00", "01", "02", "03", "04", "05", "06")
initsecs <- 59 + 60 * (17 + 60 * 21)
mysecs <- (myhrs * 3600)
secsat <- (mysecs - initsecs)/16
mxy <- max(BabyECG)
mny <- min(BabyECG)
ro <- range(BabySS)
no <- ((mxy - mny) * (BabySS - ro[1]))/(ro[2] - ro[1]) + mny
rc <- 0:4
nc <- ((mxy - mny) * (rc - ro[1]))/(ro[2] - ro[1]) + mny
plot(1:length(BabyECG), BabyECG, xaxt = "n", type = "l", xlab = "Time (hours)", ylab = "Heart rate (beats per minute)")

lines(1:length(BabyECG), no, lty = 3)
axis(1, at = secsat, labels = mylab)
axis(4, at = nc, labels = as.character(rc))

#library(LSTS)
#block.smooth.periodogram(BabyECG, spar.freq = 0.8, spar.time = 0.8, theta = 30, phi = 30)
isize <- length(BabyECG)
noctave <- 5
nvoice <- 8
pp <- noctave * nvoice
cwtBabyECG <- cwt(BabyECG, noctave, nvoice, plot=FALSE)
image(1:isize, seq(0, pp, length = pp), Mod(cwtBabyECG), xlab = "Time", ylab = "log(scale)",main = "Wavelet Transform Modulus",col = grey(seq(0, 1, length = 256)))

isize <- length(BabyECG)
nvoice <- 100
freqstep <- .001
scale <- 7
cgtBabyECG <-cgt(BabyECG,nvoice,freqstep,scale,plot=FALSE)
image(1:isize, seq(0, nvoice * freqstep/2, length = nvoice), Mod(cgtBabyECG), xlab = "Time", ylab = "Frequency",col = grey(seq(0, 1, length = 256)))
        title("Gabor Transform Modulus")

Data/Theory goals of the lecture

The goal is to describe locally stationary time series (LSTS) that can be used as models for the examples above. Two approaches will be considered, the first in the time/spectral domain and the second in the wavelet domain.

One basic view is that discrete time \(t=1,\ldots,T\) will be rescaled to \(u = \frac{t}{T} \in [0,1]\), and the parameters allowed to depend on \(t/T\). From theoretical perspective, conditions are imposed on relevant “ideal” quantities over \(u\in[0,1]\).

What makes it “locally stationary”?

LSTS: time and spectral domain approaches

For example, an LS-ARMA\((p,q)\) LSTS is defined as \(Y_{t,T}\), \(t=1,\ldots,T\), satisfying \[ \Phi(\frac{t}{T},B) Y_{t,T} = \Theta(\frac{t}{T},B) \sigma(\frac{t}{T})Z_t, \] where \(\{Z_t\}\sim \mbox{WN}(0,1)\), \[ \Phi(u,B) = 1 - \phi_1(u) B - \ldots - \phi_p(u) B^p,\quad \Theta(u,B) = 1 + \theta_1(u) B + \ldots + \theta_q(u) B^q. \] A parametric form could be used for the functions \(\phi_j(u)\) and \(\theta_k(u)\), for example, \(\phi_1(u) = \alpha_0 + \alpha_1 u\).

Estimation can be carried out using Kalman filter (or maximum likelihood estimation in the spectral domain). See Dahlhaus (1997), Palma et al. (2013) (including prediction) for details.

More generally, the LSTS model is written in the spectral domain as \[ X_{t,T} = \mu(\frac{t}{T}) + \int_{-\pi}^\pi e^{i\lambda t} A^0_{t,T}(\lambda) Z(d\lambda). \] Estimation, examples, assumptions and other issues are discussed in Dahlhaus (1997) among others. For example, one basic assumption is to require \[ A^0_{t,T}(\lambda) \approx A(\frac{t}{T},\lambda), \] for a suitable function \(A(u,\lambda)\). The quantity \[ f(u,\lambda) = |A(u,\lambda)|^2 \] is the time varying spectral density of the model, and can be thought as what one would get by locally smoothing the spectrogram. See Dahlhaus (1997).

Example

The autoregressive model with time varying coefficients satisfies \[ \sum_{j=0}^p a_j(\frac{t}{T}) X_{t-j,T} = \sigma(\frac{t}{T}) Z_t\quad \quad (a_0(u)=1) \] has the time varying spectral density \[ f(u,\lambda) = \frac{\sigma^2(u)}{2\pi} \Big| \sum_{j=0}^p a_j(u) e^{-i\lambda j} \Big|^2. \] In the simulated example above, \(p=1\) and \(a_1(u) = 0.8u\).

Back to data example

Tree rings

## start parameters
phi = smooth.spline(table$phi, spar = spar)$y
fit.1 = nls(phi ~ a0+a1*u, start = list(a0 = 0.65, a1 = 0.00))
sigma = smooth.spline(table$sigma, spar = spar)$y
fit.2 = nls(sigma ~ b0+b1*u, start = list(b0 = 0.65, b1 = 0.00))
fit = LS.whittle(series = y, start = c(coef(fit.1), coef(fit.2)), order = c(p = 1, q = 0), ar.order = c(1), sd.order = 1, N = 180, n.ahead = 10)

## Summary
LS.summary(fit)
## $summary
##    Estimate Std. Error z-value Pr(>|z|)
## a0   0.5006     0.0709  7.0606   0.0000
## a1   0.2577     0.1263  2.0410   0.0413
## b0   0.1130     0.0028 40.7288   0.0000
## b1  -0.0122     0.0051 -2.3844   0.0171
## 
## $aic
## [1] -5.313435
## 
## $npar
## [1] 4
## Diagnostic
ts.diag(fit$residuals)

## Fitted Values

op = par(mfrow = c(1,2), cex = 0.8)
spar = 0.6
plot(smooth.spline(phi, spar = spar)$y ~ u, data = table, pch = 20, ylim = c(0,1),
       xlim = c(0,1), las = 1, bty = "n", ylab = expression(phi(u)))
lines(fit$coef[1]+fit$coef[2]*u~u)
plot(smooth.spline(sigma, spar = spar)$y ~ u, data = table, pch = 20, ylim = c(0,0.2),
       xlim = c(0,1), las = 1, bty = "n", ylab = expression(sigma(u)))
lines(fit$coef[3]+fit$coef[4]*u~u)

par(op)
plot(fit$series, type = "l")
lines(fit$fitted.values, col = "red", lty = 1)
for(alpha in seq(0.05,1.00, 0.01)){
  lines(fit$pred + fit$se * qnorm(1-alpha/2), col = gray(1-alpha))
  lines(fit$pred - fit$se * qnorm(1-alpha/2), col = gray(1-alpha))
  }

LSTS: wavelet domain approaches

In the wavelet domain, on the other hand, a (“zero-mean”) LSTS \(X_{t,T}\), \(t=1,\ldots,T=2^J\), can be taken as \[ X_{t,T} = \sum_{j=-J}^{-1} \sum_k w_{j,k;T} \psi_{j,k}(t) \xi_{j,k}, \] where \(\{\xi_{j,k}\}\) is a WN series, \(\{w_{j,k;T}\}\) is a set of weights (amplitudes), and \(\{\psi_{j,k}(t) \}\) is a set of discrete non-decimated wavelets. A basic assumption is that \[ w_{j,k;T} \approx W_j(\frac{t}{T}) \] for suitable functions \(W_j(u)\). See Nason (2008), Section 5.3, who also uses the different convention for scales with \(j\) replaced by \(-j\). The quantity \[ S_j(u) = |W_j(u)|^2 \] is referred to as evoltionary (time varying) wavelet spectrum. (After suitable correction) It can be thought as what one would get by locally smoothing the scalogram of non-decimated DWT. This is all analogous to the time/spectral domain perspective above.

There are no parametric models but the model is useful for simulation, prediction and other tasks.

Example

Take \(T=1024\) and \[ S_j(u) = \left\{ \begin{array}{cl} \sin^2(4\pi u) & \mbox{for}\ j=-6,u \in(0,1),\\ 1 & \mbox{for}\ j=-1, u\in (800/1024, 900/1024),\\ 0 & \mbox{else}. \end{array} \right. \]

myspec <- cns(1024)
myspec <- putD(myspec,level=4,sin(seq(from=0, to=4*pi, length=1024))^2)
burstat800 <- c(rep(0,800), rep(1,100), rep(0,124))
myspec <-putD(myspec, level=9, v=burstat800)
plot(myspec, main="", sub="", ylabchars=-(1:10),ylab="Scale")

##  [1] 1 1 1 1 1 1 1 1 1 1
myproc <- LSWsim(myspec)
plot.ts(myproc, ylab="myproc, X_t")

load("WaveletFiguresJul15.RData")

f.tsa6()

f.tsa7()

f.tsa6 <-
function(cex=1){

op <- par(cex=cex)
plot(tsa6.ansWP, sub="", main="", ylabchars=-(1:10), ylab="Scale", xlab="Time")
par(op)

}

f.tsa7 <-
function(cex=1){

op <- par(cex=cex)
plot(tsa7.ansS, main="", sub="", ylabchars=-(1:10), ylab="Scale", xlab="Time")
par(op)


}

Back to data example

Baby ECG

 # Start time (in hours)
   sttm <- 21+(17+59/60)/60
   # Labels for x axis
   tchar <- c("22", "23", "00", "01", "02", "03", "04", "05", "06")
   # Numerical values for x axis, convert them to x values
   tm2 <- c(22,23, 24, 25, 26, 27, 28, 29, 30)
   tm2 <- tm2 - sttm
   # Convert to seconds
   tm2s <- tm2*60*60
   # Convert to sampling units
   tm2u <- tm2s/16
   # Plot the differenced data
   #

   ts.plot(diff(BabyECG))

   # Now prepare differenced data for analysis (add on an
   # observation at the front so as to make the differenced
   # series an appropriate length).
   dBabyECG <- diff(c(BabyECG[2], BabyECG))
   # Compute corrected smoothed wavelet periodgram with  options
   spec <- ewspec(dBabyECG,smooth.levels=4:10,smooth.policy="universal", smooth.transform=log,smooth.inverse=exp)$S
   # Plot the estimate
   plot(spec, main="", sub="", ylabchars=-(1:11),scaling="by.level", ylab="Scale",xlab="Time (hours)",xlabvals=tm2u, xlabchars=tchar)

##  [1] 7.631794e+01 9.290821e+00 3.637686e+00 2.281380e+00 5.918848e-01
##  [6] 8.920871e-02 1.919850e-01 4.944353e-02 2.407098e-02 4.047056e-03
## [11] 7.316527e-04
FineCoefs <- accessD(spec, lev=10)

myhrs <- c(22, 23, 24, 25, 26, 27, 28, 29, 30)
mylab <- c("22", "23", "00", "01", "02", "03", "04", "05", "06")
initsecs <- 59 + 60 * (17 + 60 * 21)
mysecs <- (myhrs * 3600)
secsat <- (mysecs - initsecs)/16
mxy <- max(FineCoefs)
mny <- min(FineCoefs)
ro <- range(BabySS)
no <- ((mxy - mny) * (BabySS - ro[1]))/(ro[2] - ro[1]) + mny
rc <- 0:4
nc <- ((mxy - mny) * (rc - ro[1]))/(ro[2] - ro[1]) + mny

plot(1:length(BabyECG), FineCoefs, xaxt = "n", type = "l", xlab = "Time (hours)", ylab = "EWS at finest scale")
lines(1:length(BabyECG), no, lty = 3)
axis(1, at = secsat, labels = mylab)
axis(4, at = nc, labels = as.character(rc))

Questions/Homework

Q0: Try some of the considered LSTS methods on real time series.

Q1: Explain how exactly estimation of LSTS models (e.g. the LS-ARMA models) is carried out in the spectral domain. (See Dahlhaus (1997).)

Q2: Explain what correction is made to go from the non-decimated DWT to the estimate of the evolutionary wavelet spectrum. (See Nason (2008).)

References

Dahlhaus, R. 1997. “Fitting Time Series Models to Nonstationary Processes.” The Annals of Statistics 25 (1): 1–37.

Nason, G. P. 2008. Wavelet Methods in Statistics with R. Use R! Springer, New York.

Palma, W., R. Olea, and G. Ferreira. 2013. “Estimation and Forecasting of Locally Stationary Processes.” Journal of Forecasting 32 (1): 86–96.