inst/doc/multivariate.R

In psd: Adaptive, Sine-Multitaper Power Spectral Density and Cross Spectrum Estimation

## ----eval=TRUE, echo=FALSE, label='prereq'----------------------------
library(knitr)
options(width=72)
knitr::opts_chunk$set(tidy = TRUE, size="small", split=TRUE, par=TRUE, fig.width=6)
knit_hooks$set(par = function(before, options, envir){
  if (before && options$fig.show!='none'){
    par(cex.lab=.95,cex.axis=.9,mgp=c(2.7,.6,0),tcl=-.4)
    par(las=1, lend='square', ljoin='mitre')
  }
})

rm(list=ls())
rng_et <- c(0.7, 2.2)
rng_ba <- c(0.0, 10)

## ----eval=TRUE, echo=TRUE, label='loadwipp'---------------------------
library(psd)
data(wipp30)
str(wipp30)
dat <- wipp30[, c('baro','et','wl')] # output as last column, inputs as first two
dat <- window(ts(dat), 100, 2400)
head(dat)
freqs <- 1:(nrow(dat)/2) * (24) / nrow(dat) # frequency in cpd

## ----eval=TRUE, echo=TRUE, label="calcmv"-----------------------------
mv <- pspectrum(dat, riedsid_column= 0L, verbose = FALSE, plot = FALSE)

## ----eval=TRUE, echo=TRUE, label="names"------------------------------
names(mv)

## ----eval=TRUE, echo=TRUE, label="plotnames"--------------------------
pspectrum(dat, riedsid_column= 0L, verbose = FALSE, plot = TRUE)

## ----eval=TRUE, echo=TRUE, label="Tapers"-----------------------------

mn <- pspectrum(dat, riedsid_column= 0L, verbose = FALSE)   # minimum number 
mx <- pspectrum(dat, riedsid_column=-1L, verbose = FALSE)   # maximum number
c1 <- pspectrum(dat, riedsid_column= 1L, verbose = FALSE)   # column 1
c2 <- pspectrum(dat, riedsid_column= 2L, verbose = FALSE)   # column 2
c3 <- pspectrum(dat, riedsid_column= 3L, verbose = FALSE)   # column 3


## ----eval=TRUE, echo=FALSE, label="Taperplot", fig.height = 4.8, fig.width=5.6----

ylim <- c(0, max(mx$taper))
layout(matrix(1:6, ncol = 1, nrow = 6), heights = c(rep(0.19, 5), 0.05))
par(mar = c(0.5, 4, 0, 0), mgp = c(2, 0.5, 0))
plot(mn$taper, ylim = ylim, xaxt="n")
legend('topleft', 'minimum')
plot(mx$taper, ylim = ylim, xaxt="n")
legend('topleft', 'maximum')
plot(c1$taper, ylim = ylim, xaxt="n")
legend('topleft', 'Water level column')
plot(c2$taper, ylim = ylim, xaxt="n")
legend('topleft', 'Barometric column')
plot(c3$taper, ylim = ylim)
legend('topleft', 'Earth tide column')
mtext('taper index', 1, line = 1.2, cex = 0.8)


## ----eval=TRUE, echo=TRUE, label="taperspecification", fig.height =5.5----
par(mar = c(3, 4, 1, 2), mgp = c(2, 0.5, 0))
layout(matrix(1:2), heights=c(3,1.2))
nti <- 11
one_val <- pspectrum(dat, riedsid_column= 0L, verbose = FALSE, plot = TRUE,
                     niter = 0, 
                     ntap.init = nti)
with(one_val,{plot(taper, xi=freq, xlab='frequency')})
axis(4, at=nti)

## ----eval=TRUE, echo=TRUE, label="taperspecification2", fig.height = 5.5----
par(mar = c(3, 4, 1, 2), mgp = c(2, 0.5, 0), xaxs='i')
layout(matrix(1:2), heights=c(3,1.2))
vec_val <- pspectrum(dat, riedsid_column= 0L, verbose = FALSE, plot = TRUE,
                 niter = 0, 
                 ntap = as.tapers(seq_len(nrow(dat) %/% 2) %/% 7 + 7))

with(vec_val,{plot(taper, xi=freq, xlab='frequency')})
axis(4, at=nti)

## ----eval=TRUE, echo=TRUE, label="pspec"------------------------------

spec <- pspectrum(dat, riedsid_column= 0L, verbose = FALSE)$pspec
dim(spec)


## ----eval=TRUE, echo=TRUE, label="coh"--------------------------------
scoh <- pspectrum(dat, riedsid_column= 0L, 
                 verbose = FALSE,
                 niter = 0, ntap.init = 11)
coh <- scoh[['coh']]
scoh_base <- stats::spec.pgram(dat[, c(3,1,2)], spans = c(11), 
                              fast = FALSE, plot = FALSE)
coh_base <- scoh_base[['coh']]

## ----eval=TRUE, echo=FALSE, label="cohbaplot", fig.height=2.3---------
par(mar = c(3, 3, 2, 0), mgp = c(2, 1, 0))
plot(coh[,1]~freqs, type='l', 
     ylab = 'WL and Barometric',
     xlab = 'Frequency (cpd)',
     xlim = rng_ba,
     main = 'Coherence')
points(y = coh_base[,1], x = freqs, type='l',  col = 'red')
legend('topright', 
       legend = c('psd::pspectrum', 'stats::spec.pgram'), 
       col = c('black', 'red'), lty = 1, cex=0.8)

## ----eval=TRUE, echo=FALSE, label="cohetplot", fig.height=2.3---------
par(mar = c(3, 3, 2, 0), mgp = c(2, 1, 0))

plot(coh[,2]~freqs, type='l', 
     ylab = 'WL and Earth tide',
     xlab = 'Frequency (cpd)',
     xlim = rng_et,
     main = 'Coherence')
points(y = coh_base[,2], x = freqs, type='l',  col = 'red')



## ----eval=TRUE, echo=TRUE, label="phase"------------------------------
phase <- pspectrum(dat, riedsid_column= 0L, verbose = FALSE,
                 niter = 0, ntap.init = 11)$phase
phase_base <- spec.pgram(dat[, c(3,1,2)], spans = c(11), fast = FALSE, plot = FALSE)$phase

## ----eval=TRUE, echo=FALSE, label="phasebaplot", fig.height=2.3-------
par(mar = c(3, 3, 2, 0), mgp = c(2, 1, 0))
plot(phase[,1]~freqs, type='l', 
     ylab = 'WL and Barometric',
     xlab = 'Frequency (cpd)',
     xlim = rng_ba,
     main = 'Phase Difference')
points(y = phase_base[,1], x = freqs, type='l',  col = 'red')
legend('bottomleft', 
       legend = c('psd', 'base'), inset=c(0.1,0),
       col = c('black', 'red'), lty = 1, cex=0.8)

## ----eval=TRUE, echo=FALSE, label="phaseetplot", fig.height=2.3-------
par(mar = c(3, 3, 2, 0), mgp = c(2, 1, 0))

plot(phase[,2]~freqs, type='l', 
     ylab = 'WL and Earth tide',
     xlab = 'Frequency (cpd)',
     xlim = rng_et,
     main = 'Phase')
points(y = phase_base[,2], x = freqs, type='l',  col = 'red')



## ----eval=TRUE, echo=TRUE, label="transfergain", fig.height=3.3-------
transfer <- pspectrum(dat, riedsid_column= -1L, verbose = FALSE)$transfer
gain <- Mod(transfer)

## ----eval=TRUE, echo=FALSE, label="transfergainplot", fig.height=2.3----


par(mar = c(3, 3, 2, 0), mgp = c(2, 1, 0))
plot(gain[, 1]~freqs, type='l', 
     xlim = rng_ba, 
     ylim = c(0, 1),
     ylab = 'WL and Barometric',
     xlab = 'Frequency (cpd)',
     main = 'Gain')
ysc <- 1e5
plot(ysc*gain[, 2]~freqs, type='l', 
     xlim = rng_et,
     ylim = ysc*c(0, 0.00005),
     ylab = 'WL and Earth tide',
     xlab = 'Frequency (cpd)',
     main = parse(text=sprintf('bold("Gain") ~ "(times" ~ 10^%s * ")"',log10(ysc))))

## ----eval=TRUE, echo=TRUE, label="transferphase"----------------------
phase <- Arg(transfer)

## ----eval=TRUE, echo=FALSE, label="transferphaseplot", fig.height=2.3----
par(mar = c(3, 3, 2, 0), mgp = c(2, 1, 0))

plot(phase[, 1]~freqs, type='l', xlim = rng_ba,
     ylab = 'WL and Barometric',
     xlab = 'Frequency (cpd)',
     main = 'Phase')
plot(phase[, 2]~freqs, type='l', xlim = rng_et, 
     ylab = 'WL and Earth tide',
     xlab = 'Frequency (cpd)',
     main = 'Phase')

## ----eval=TRUE, echo=TRUE, label=SI-----------------------------------
utils::sessionInfo()