Monthly time-series of water flow at Big Tree water gauge.

Description

Average monthly natural water flow (cubic feet per second) at the Big Tree gauge of the San Lorenzo river in Santa Cruz County, CA from 1937 through 2014.

Usage

BTflow

Format

A time series of length 936.

Source

https://waterdata.usgs.gov/nwis/

References

U.S. Geological Survey (2016). National Water Information System data available on the World Wide Web (USGS Water Data for the Nation). https://waterdata.usgs.gov/nwis/.

Daily time-series of ELI anomalies.

Description

ELI anomalies on the daily time-scale from January 1, 1979 to December 31, 2019 with all February 29ths omitted.

Usage

ELIanoms

Format

A time series of length 14965.

Source

https://portal.nersc.gov/archive/home/projects/cascade/www/ELI

References

Patricola, C.M., O’Brien, J.P., Risser, M.D. et al. Maximizing ENSO as a source of western US hydroclimate predictability. Clim Dyn 54, 351–372 (2020). doi:10.1007/s00382-019-05004-8

Combines state space blocks of an exDQLM

Description

The function combines two models into a single state space model for an exDQLM.

Usage

combineMods(m1, m2)

Arguments

Value

List containing the new combined state space model components.

Examples

trend.comp = polytrendMod(2,rep(0,2),10*diag(2))
seas.comp = seasMod(365,c(1,2,4),C0=10*diag(6))
model = combineMods(trend.comp,seas.comp)
# using dlm package
library(dlm)
model = combineMods(dlmModPoly(order=2,C0=10*diag(2)),dlmModTrig(365,2,C0=10*diag(4)))

Plot a component of an exDQLM

Description

The function plots the dynamic MAP estimates and 95

Usage

compPlot(
  y,
  m1,
  index,
  add = FALSE,
  col = "purple",
  just.theta = FALSE,
  cr.percent = 0.95
)

Arguments

Value

A list of the following is returned:

• 'map.comp' - MAP estimate of the dynamic component (or element of the state vector).

• 'lb.comp' - Lower bound of the 95

• 'ub.comp' - Upper bound of the 95

Examples

y = scIVTmag[1:365]
trend.comp = polytrendMod(2,rep(0,2),10*diag(2))
seas.comp = seasMod(365,c(1,2,4),C0=10*diag(6))
model = combineMods(trend.comp,seas.comp)
M0 = exdqlmISVB(y,p0=0.85,model,df=c(0.98,1),dim.df = c(2,6),
                   gam.init=-3.5,sig.init=15,tol=0.05)
# plot first harmonic component
compPlot(y,M0,index=c(3,4),col="blue")

Create state space model of exDQLM from DLM

Description

The function creates a state space model of an exDQLM from "'dlm'" object.

Usage

dlmMod(m)

Arguments

Value

List containing only the components of 'm' needed for the exDQLM state space model.

Examples

library(dlm)
m = dlmModPoly(order=2,C0=10*diag(2)) + dlmModTrig(365,2,C0=10*diag(4))
model = dlmMod(m)

exDQLM Diagnostics

Description

The function computes the following for the model(s) provided: the posterior predictive loss criterion based off the check loss, the one-step-ahead distribution sequence and its KL divergence from normality. The function also plots the following: the qq-plot and ACF plot corresponding to the one-step-ahead distribution sequence, and a time series plot of the MAP standard forecast errors.

Usage

exdqlmChecks(
  y,
  m1,
  m2 = NULL,
  plot = TRUE,
  cols = c("grey", "grey"),
  ref = NULL
)

Arguments

Value

A list containing the following is returned:

• 'm1.uts' - The one-step-ahead distribution sequence of 'm1'.

• 'm1.KL' - The KL divergence of 'm1.uts' and a standard normal.

• 'm1.pplc' - The posterior predictive loss criterion of 'm1' based off the check loss function.

• 'm1.qq' - The ordered pairs of the qq-plot comparing 'm1.uts' with a standard normal distribution.

• 'm1.acf' - The autocorrelations of 'm1.uts' by lag.

If 'm2' is provided, analogous results for 'm2' are also included in the list.

Examples

y = scIVTmag[1:100]
model = polytrendMod(1,mean(y),10)
M0 = exdqlmISVB(y,p0=0.85,model,df=c(0.95),dim.df = c(1),
                  gam.init=-3.5,sig.init=15)
check.out = exdqlmChecks(y,M0,plot=FALSE)
check.out$m1.KL
check.out$m1

k-step-ahead quantile forecasts

Description

Computes filtered and k-step-ahead forecast quantiles from a fitted dynamic quantile model and optionally adds them to an existing plot.

Usage

exdqlmForecast(
  y,
  start.t,
  k,
  m1,
  fFF = NULL,
  fGG = NULL,
  plot = TRUE,
  add = FALSE,
  cols = c("purple", "magenta"),
  cr.percent = 0.95
)

Arguments

Value

A list with components:

• fa Forecast state mean vectors (p \times k matrix).

• fR Forecast state covariance matrices (p \times p \times k array).

• ff Forecast quantile means (length-k numeric).

• fQ Forecast quantile variances (length-k numeric).

Examples

 # Toy example; keep small and fast
 y <- scIVTmag[1:100]
 model <- polytrendMod(1, stats::quantile(y, 0.85), 10)
 M0 <- exdqlmISVB(y, p0 = 0.85, model, df = c(0.98), dim.df = c(1),
                  gam.init = -3.5, sig.init = 15)
 exdqlmForecast(y, start.t = 90, k = 10, m1 = M0)

exDQLM - ISVB algorithm

Description

The function applies an Importance Sampling Variational Bayes (ISVB) algorithm to estimate the posterior of an exDQLM.

Usage

exdqlmISVB(
  y,
  p0,
  model,
  df,
  dim.df,
  fix.gamma = FALSE,
  gam.init = NA,
  fix.sigma = TRUE,
  sig.init = NA,
  dqlm.ind = FALSE,
  exps0,
  tol = 0.1,
  n.IS = 500,
  n.samp = 200,
  PriorSigma = NULL,
  PriorGamma = NULL,
  verbose = TRUE
)

Arguments

Value

A list of the following is returned:

• 'run.time' - Algorithm run time in seconds.

• 'iter' - Number of iterations until convergence was reached.

• 'dqlm.ind' - Logical value indicating whether gamma was fixed at '0', reducing the exDQLM to the special case of the DQLM.

• 'model' - List of the state-space model including 'GG', 'FF', prior parameters 'm0' and 'C0'.

• 'p0' - The quantile which was estimated.

• 'df' - Discount factors used for each block.

• 'dim.df' - Dimension used for each block of discount factors.

• 'sig.init' - Initial value for sigma, or value at which sigma was fixed if 'fix.sigma=TRUE'.

• 'seq.sigma' - Sequence of sigma estimated by the algorithm until convergence.

• 'samp.theta' - Posterior sample of the state vector variational distribution.

• 'samp.post.pred' - Sample of the posterior predictive distributions.

• 'map.standard.forecast.errors' - MAP standardized one-step-ahead forecast errors.

• 'samp.sigma' - Posterior sample of scale parameter sigma variational distribution.

• 'samp.vts' - Posterior sample of latent parameters, v_t, variational distributions.

• 'theta.out' - List containing the variational distribution of the state vector including filtered distribution parameters ('fm' and 'fC') and smoothed distribution parameters ('sm' and 'sC').

• 'vts.out' - List containing the variational distributions of latent parameters v_t.

If 'dqlm.ind=FALSE', the list also contains:

• 'gam.init' - Initial value for gamma, or value at which gamma was fixed if 'fix.gamma=TRUE'.

• 'seq.gamma' - Sequence of gamma estimated by the algorithm until convergence.

• 'samp.gamma' - Posterior sample of skewness parameter gamma variational distribution.

• 'samp.sts' - Posterior sample of latent parameters, s_t, variational distributions.

• 'gammasig.out' - List containing the IS estimate of the variational distribution of sigma and gamma.

• 'sts.out' - List containing the variational distributions of latent parameters s_t.

Or if 'dqlm.ind=TRUE', the list also contains:

• 'sig.out' - List containing the IS estimate of the variational distribution of sigma.

Examples

y = scIVTmag[1:1095]
trend.comp = polytrendMod(1,mean(y),10)
seas.comp = seasMod(365,c(1,2,4),C0=10*diag(6))
model = combineMods(trend.comp,seas.comp)
M0 = exdqlmISVB(y,p0=0.85,model,df=c(1,1),dim.df = c(1,6),
                 gam.init=-3.5,sig.init=15,tol=0.05)

exDQLM - MCMC algorithm

Description

The function applies a Markov chain Monte Carlo (MCMC) algorithm to sample the posterior of an exDQLM.

Usage

exdqlmMCMC(
  y,
  p0,
  model,
  df,
  dim.df,
  fix.gamma = FALSE,
  gam.init = NA,
  fix.sigma = FALSE,
  sig.init = NA,
  dqlm.ind = FALSE,
  Sig.mh,
  joint.sample = FALSE,
  n.burn = 2000,
  n.mcmc = 1500,
  init.from.isvb = TRUE,
  PriorSigma = NULL,
  PriorGamma = NULL,
  verbose = TRUE
)

Arguments

Value

A list of the following is returned:

• 'run.time' - Algorithm run time in seconds.

• 'model' - List of the state-space model including 'GG', 'FF', prior parameters 'm0' and 'C0'.

• 'p0' - The quantile which was estimated.

• 'df' - Discount factors used for each block.

• 'dim.df' - Dimension used for each block of discount factors.

• 'samp.theta' - Posterior sample of the state vector.

• 'samp.post.pred' - Sample of the posterior predictive distributions.

• 'map.standard.forecast.errors' - MAP standardized one-step-ahead forecast errors.

• 'samp.sigma' - Posterior sample of scale parameter sigma.

• 'samp.vts' - Posterior sample of latent parameters, v_t.

• 'theta.out' - List containing the distributions of the state vector including filtered distribution parameters ('fm' and 'fC') and smoothed distribution parameters ('sm' and 'sC').

If 'dqlm.ind=FALSE', the list also contains the following:

• 'samp.gamma' - Posterior sample of skewness parameter gamma.

• 'samp.sts' - Posterior sample of latent parameters, s_t.

• 'init.log.sigma' - Burned samples of log sigma from the random walk MH joint sampling of sigma and gamma.

• 'init.logit.gamma' - Burned samples of logit gamma from the random walk MH joint sampling of sigma and gamma.

• 'accept.rate' - Acceptance rate of the MH step.

• 'Sig.mh' - Covariance matrix used in MH step to jointly sample sigma and gamma.

Examples

y = scIVTmag[1:100]
trend.comp = polytrendMod(1,mean(y),10)
seas.comp = seasMod(365,c(1,2,4),C0=10*diag(6))
model = combineMods(trend.comp,seas.comp)
M2 = exdqlmMCMC(y,p0=0.85,model,df=c(1,1),dim.df = c(1,6),
                gam.init=-3.5,sig.init=15,
                n.burn=100,n.mcmc=150)

Plot exDQLM

Description

The function plots the MAP estimates and 95

Usage

exdqlmPlot(y, m1, add = FALSE, col = "purple", cr.percent = 0.95)

Arguments

Value

A list of the following is returned:

• 'map.quant' - MAP estimate of the dynamic quantile.

• 'lb.quant' - Lower bound of the 95

• 'ub.quant' - Upper bound of the 95

Examples

y = scIVTmag[1:100]
model = polytrendMod(1,quantile(y,0.85),10)
M0 = exdqlmISVB(y,p0=0.85,model,df=c(0.98),dim.df = c(1),
                   gam.init=-3.5,sig.init=15)
exdqlmPlot(y,M0,col="blue")

Monthly Niño 3.4 Index.

Description

Monthly Niño 3.4 sea surface temperature (SST) Index time series.

Usage

nino34

Format

A time series of length 936.

Source

https://psl.noaa.gov/gcos_wgsp/Timeseries/Nino34/

References

Rayner N. A., D. E. Parker, E. B. Horton, C. K. Folland, L. V. Alexander, D. P. Rowell, E. C. Kent, A. Kaplan. Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century. J. Geophys. Res., 108 (D14), 4407 (2003). doi:10.1029/2002JD002670

Create an n-th order polynomial exDQLM component

Description

The function creates an n-th order polynomial exDQLM component.

Usage

polytrendMod(order, m0, C0)

Arguments

Value

A list of the following:

• FF - Observational vector.

• GG - Evolution matrix.

• m0 - Prior mean of the state vector.

• C0 - Prior covariance of the state vector.

Examples

# create a second order polynomial component
trend.comp = polytrendMod(2,rep(0,2),10*diag(2))

Time series of daily average magnitude IVT in Santa Cruz, CA.

Description

ECMWF Re-Analysis 5 (ERA5) daily average magnitude IVT in Santa Cruz, CA (approximately 22 N, 122 W) from January 1, 1979 to December 31, 2019 with all February 29ths omitted.

Usage

scIVTmag

Format

A time series of length 14965.

Source

https://cds.climate.copernicus.eu

References

Hersbach, H, Bell, B, Berrisford, P, et al. The ERA5 global reanalysis. Q J R Meteorol Soc. 2020; 146: 1999– 2049. doi:10.1002/qj.3803

Create Fourier representation of a periodic exDQLM component

Description

The function creates a Fourier form periodic component for given period and harmonics.

Usage

seasMod(p, h, m0, C0)

Arguments

Value

A list of the following:

• FF - Observational vector.

• GG - Evolution matrix.

• m0 - Prior mean of the state vector.

• C0 - Prior covariance of the state vector.

Examples

# create a seasonal component with first, second and fourth harmonics of a period of 365
seas.comp = seasMod(365,c(1,2,4),C0=10*diag(6))

Transfer Function exDQLM - ISVB algorithm

Description

The function applies an Importance Sampling Variational Bayes (ISVB) algorithm to estimate the posterior of an exDQLM with exponential decay transfer function component.

Usage

transfn_exdqlmISVB(
  y,
  p0,
  model,
  X,
  df,
  dim.df,
  lam,
  tf.df,
  fix.gamma = FALSE,
  gam.init = NA,
  fix.sigma = TRUE,
  sig.init = NA,
  dqlm.ind = FALSE,
  exps0,
  tol = 0.1,
  n.IS = 500,
  n.samp = 200,
  PriorSigma = NULL,
  PriorGamma = NULL,
  tf.m0 = rep(0, 2),
  tf.C0 = diag(1, 2),
  verbose = TRUE
)

Arguments

Value

A list of the following is returned:

• 'run.time' - Algorithm run time in seconds.

• 'iter' - Number of iterations until convergence was reached.

• 'dqlm.ind' - Logical value indicating whether gamma was fixed at '0', reducing the exDQLM to the special case of the DQLM.

• 'model' - List of the augmented state-space model including 'GG', 'FF', prior parameters 'm0' and 'C0'.

• 'p0' - The quantile which was estimated.

• 'df' - Discount factors used for each block, including transfer function component.

• 'dim.df' - Dimension used for each block of discount factors, including transfer function component.

• 'lam' - Transfer function rate parameter lambda.

• 'sig.init' - Initial value for sigma, or value at which sigma was fixed if 'fix.sigma=TRUE'.

• 'seq.sigma' - Sequence of sigma estimated by the algorithm until convergence.

• 'samp.theta' - Posterior sample of the state vector variational distribution.

• 'samp.post.pred' - Sample of the posterior predictive distributions.

• 'map.standard.forecast.errors' - MAP standardized one-step-ahead forecast errors.

• 'samp.sigma' - Posterior sample of scale parameter sigma variational distribution.

• 'samp.vts' - Posterior sample of latent parameters, v_t, variational distributions.

• 'theta.out' - List containing the variational distribution of the state vector including filtered distribution parameters ('fm' and 'fC') and smoothed distribution parameters ('sm' and 'sC').

• 'vts.out' - List containing the variational distributions of latent parameters v_t.

• 'median.kt' - Median number of time steps until the effect of X_t is less than or equal to 1e-3.

If 'dqlm.ind=FALSE', the list also contains:

• 'gam.init' - Initial value for gamma, or value at which gamma was fixed if 'fix.gamma=TRUE'.

• 'seq.gamma' - Sequence of gamma estimated by the algorithm until convergence.

• 'samp.gamma' - Posterior sample of skewness parameter gamma variational distribution.

• 'samp.sts' - Posterior sample of latent parameters, s_t, variational distributions.

• 'gammasig.out' - List containing the IS estimate of the variational distribution of sigma and gamma.

• 'sts.out' - List containing the variational distributions of latent parameters s_t.

Or if 'dqlm.ind=TRUE', the list also contains:

• 'sig.out' - List containing the IS estimate of the variational distribution of sigma.

Examples

y = scIVTmag[1:1095]
X = ELIanoms[1:1095]
trend.comp = polytrendMod(1,mean(y),10)
seas.comp = seasMod(365,c(1,2,4),C0=10*diag(6))
model = combineMods(trend.comp,seas.comp)
M1 = transfn_exdqlmISVB(y,p0=0.85,model=model,
                          X,df=c(1,1),dim.df = c(1,6),
                          gam.init=-3.5,sig.init=15,
                          lam=0.38,tf.df=c(0.97,0.97))