Bayesian spectral inference for time series

Description

Bayesian parametric, nonparametric and semiparametric procedures for spectral density inference of univariate (locally) stationary time series and multivariate stationary time series

Details

The package contains several methods (parametric, nonparametric and semiparametric) for Bayesian spectral density inference. The main algorithms to fit the models for univariate stationary time series are:

• gibbs_ar: Parametric, autoregressive (AR) model

• gibbs_np: Nonparametric model with Whittle's likelihood and Bernstein-Dirichlet prior from Choudhuri et al (2007)

• gibbs_npc: Semiparametric model with corrected AR likelihood and Bernstein-Dirichlet prior from Kirch et al (2018)

The package also contains the following models for multivariate stationary time series:

• gibbs_var: Parametric, vector autoregressive (VAR) model

• gibbs_vnp: Nonparametric model with Whittle's likelihood and Bernstein-Hpd-Gamma prior from Meier (2018)

The main function for univariate locally stationary time series is:

• gibbs_bdp_dw: Nonparametric model with BDP-DW approach from Tang et al (2023)

as well as some useful utility functions. To get started, it is recommended to consider the examples and documentation of the functions listed above. The work was supported by DFG grants KI 1443/3-1 and KI 1443/3-2.

Author(s)

Claudia Kirch, Renate Meyer, Matthew C. Edwards, Alexander Meier, Yifu Tang

Maintainer: Renate Meyer <renate.meyer@auckland.ac.nz>

References

N. Choudhuri, S. Ghosal and A. Roy (2004) Bayesian estimation of the spectral density of a time series JASA <doi:10.1198/016214504000000557>

C. Kirch, M. C. Edwards, A. Meier and R. Meyer (2018) Beyond Whittle: Nonparametric Correction of a Parametric Likelihood with a Focus on Bayesian Time Series Analysis Bayesian Analysis <doi:10.1214/18-BA1126>

A. Meier (2018) A matrix Gamma process and applications to Bayesian analysis of multivariate time series PhD thesis, OvGU Magdeburg <doi:10.25673/13407>

Y. Tang, C. Kirch, J. E. Lee and R. Meyer (2023) Bayesian nonparametric spectral analysis of locally stationary processes ArXiv preprint <arXiv:2303.11561>

adjoint of complex matrix

Description

adjoint of complex matrix

Usage

Adj(m)

Psi-weight calculation for a VARMA model. NOTE: This is an exact copy of the MTS::PSIwgt function (only with the plot functionality removed, as not needed). This has to be done because the MTS package has been removed from CRAN in April 2022.

Description

Psi-weight calculation for a VARMA model. NOTE: This is an exact copy of the MTS::PSIwgt function (only with the plot functionality removed, as not needed). This has to be done because the MTS package has been removed from CRAN in April 2022.

Usage

PSIwgt(Phi = NULL, Theta = NULL, lag = 12, plot = TRUE, output = FALSE)

This is a nearly exact copy of the MTS::VARMAcov function, where the output commands at the end are removed. This has to be done because the function is called repeatedly within the MCMC algorithm. For future versions of the package, a better solution is intended.

Description

This is a nearly exact copy of the MTS::VARMAcov function, where the output commands at the end are removed. This has to be done because the function is called repeatedly within the MCMC algorithm. For future versions of the package, a better solution is intended.

Usage

VARMAcov_muted(Phi = NULL, Theta = NULL, Sigma = NULL, lag = 12, trun = 120)

VAR regressor matrix, see Section 2.2.3 in Koop and Korobilis (2010)

Description

VAR regressor matrix, see Section 2.2.3 in Koop and Korobilis (2010)

Usage

VAR_regressor_matrix(data, var.order)

Computing acceptance rate based on trace Note: Only use for traces from continous distributions!

Description

Computing acceptance rate based on trace Note: Only use for traces from continous distributions!

Usage

acceptanceRate(trace)

Build an nd times nd Block Toeplitz matrix from the (d times d) autocovariances gamma(0),...,gamma(n-1)

Description

Build an nd times nd Block Toeplitz matrix from the (d times d) autocovariances gamma(0),...,gamma(n-1)

Usage

acvBlockMatrix(acv)

Build an n times n Toeplitz matrix from the autocovariance values gamma(0),...,gamma(n-1)

Description

Build an n times n Toeplitz matrix from the autocovariance values gamma(0),...,gamma(n-1)

Usage

acvMatrix(acv)

Negative ARMA(p, q) log likelihood

Description

Negative ARMA(p, q) log likelihood

Usage

arma_conditional(zt, ar, ma, nll_fun, full_lik, ...)

a generic method for bdp_dw_result class

Description

a generic method for bdp_dw_result class

Usage

bayes_factor(obj)

Arguments

Extracting the Bayes factor of k1=1 from bdp_dw_result class

Description

Extracting the Bayes factor of k1=1 from bdp_dw_result class

Usage

## S3 method for class 'bdp_dw_result'
bayes_factor(obj)

Arguments

Value

the estimated Bayes factor of k1=1

Estimating the Bayes factor of hypothesis "k1 = 1".

Description

Estimating the Bayes factor of hypothesis "k1 = 1".

Usage

bdp_dw_bayes_factor_k1(post_sample, precision = 1000)

Arguments

Value

The Savage-Dickey estimate of the Bayes factor and its theoretical upper bound. c.f. section 3.3 of Tang et al. (2023).

References

Tang et al. (2023) Bayesian nonparametric spectral analysis of locally stationary processes ArXiv preprint <arXiv:2303.11561>

Calculating the estimated posterior mean, median and credible region (tv-PSD)

Description

Calculating the estimated posterior mean, median and credible region (tv-PSD)

Usage

bdp_dw_est_post_stats(post_sample, rescaled_time, freq, unif_CR = FALSE)

Arguments

Value

list containing the following fields:

MH sampler for BDP-DW method

Description

Obtain samples of the posterior of the Whittle likelihood in conjunction with a Bernstein-Dirichlet prior on the spectral density.

Usage

bdp_dw_mcmc(
  data,
  m,
  likelihood_thinning = 1,
  mcmc_params,
  prior_params,
  monitor = FALSE,
  print_interval = 100
)

Arguments

Details

Further detail can be found in the simulation study section in the references papers.

Value

list containing the following fields:

References

Y. Tang et al. (2023) Bayesian nonparametric spectral analysis of locally stationary processes ArXiv preprint <arXiv:2303.11561>

Generate a list of values for MCMC algorithm

Description

Generate a list of values for MCMC algorithm

Usage

bdp_dw_mcmc_params_gen(
  Ntotal = 110000,
  burnin = 60000,
  thin = 10,
  adaptive.batchSize = 50,
  adaptive.targetAcceptanceRate = 0.44
)

Arguments

Value

A list of MCMC parameter values

Generate a list of parameter values in prior elicitation

Description

Generate a list of parameter values in prior elicitation

Usage

bdp_dw_prior_params_gen(
  M = 1,
  g0.alpha = 1,
  g0.beta = 1,
  k1.theta = 0.01,
  k2.theta = 0.01,
  tau.alpha = 0.001,
  tau.beta = 0.001,
  k1max = 100,
  k2max = 100,
  L = 10,
  bernstein1_l = 0.1,
  bernstein1_r = 0.9,
  bernstein2_l = 0.1,
  bernstein2_r = 0.9
)

Arguments

Value

A list of prior parameter values

Construct Bernstein polynomial basis of degree k on omega

Description

Construct Bernstein polynomial basis of degree k on omega

Usage

betaBasis_k(omega, k, coarsened)

Arguments

Construct Bernstein polynomial basis of degree k on omega

Description

Construct Bernstein polynomial basis of degree k on omega

Usage

betaBasis_k_dw(omega, k, bernstein_l = 0, bernstein_r = 1)

Arguments

Value

A matrix

Examples

## Not run: 
omega <- seq(0, 1, by = 0.01)
betaBasis_k_dw(omega, 100, 0.1, 0.9)

## End(Not run)

mean center a numerical vector

Description

mean center a numerical vector

Usage

center(x, ...)

Construct coarsened Bernstein polynomial basis of degree l on omega

Description

Construct coarsened Bernstein polynomial basis of degree l on omega

Usage

coarsened_bernstein(omega, l)

Arguments

Details

See Appendix E.1 in Ghosal/Van der Vaart, Fundamentals, (2017)

Helping function for coarsened_bernstein

Description

Helping function for coarsened_bernstein

Usage

coarsened_bernstein_i(omega, l, i)

Inverse function to realValuedPsd

Description

Inverse function to realValuedPsd

Usage

complexValuedPsd(f_)

I/O: Only used within Rcpp Note: Same workaround as cx_cube_from_ComplexVector

Description

I/O: Only used within Rcpp Note: Same workaround as cx_cube_from_ComplexVector

Usage

cube_from_NumericVector(x)

I/O: Only used within Rcpp

Description

This workaround for parsing cubes was neccessary at development time because there was a (presumable) bug in RcppArmadillo that sometimes caused the parsing of arma::cx_cube objects to fail, such that the function received an un-initialized object instead of the parsed one.

Usage

cx_cube_from_ComplexVector(x)

Details

The workaround parses an Rcpp vector instead, and manually copies the data in an arma::cx_cube object. Besides being redundant, it also makes the code less readable and it is hoped that this workaround can be removed in future revisions.

Construct Bernstein polynomial basises of degree up to kmax on omega

Description

Construct Bernstein polynomial basises of degree up to kmax on omega

Usage

dbList(n, kmax, bernstein_l = 0, bernstein_r = 1, coarsened = F, verbose = F)

Arguments

Value

A list of length kmax, where the k-th list element is a matrix containing the polynomial basis of degree k

Construct Bernstein polynomial bases of degree up to kmax on omega

Description

Construct Bernstein polynomial bases of degree up to kmax on omega

Usage

dbList_dw_Bern(omega, kmax, bernstein_l = 0, bernstein_r = 1)

Arguments

Value

A list of length kmax, where the k-th list element is a matrix containing the polynomial basis of degree k

Examples

## Not run: 
omega <- seq(0, 1, by = 0.01)
dbList_dw_Bern(omega, 100, 0.1, 0.9)

## End(Not run)

Construct Bernstein polynomial bases of degree up to kmax on omega for frequency parameter lambda

Description

Construct Bernstein polynomial bases of degree up to kmax on omega for frequency parameter lambda

Usage

dbList_dw_Bern_for_lambda(
  omega,
  kmax,
  bernstein_l = 0,
  bernstein_r = 1,
  m,
  time_grid
)

Arguments

Value

A list of length kmax, where the k-th list element is a matrix containing the polynomial basis of degree k

Examples

## Not run: 
omega <- time_grid <- seq(0, 1, by = 0.01)
dbList_dw_Bern_for_lambda(omega, 100, 0.1, 0.9)

## End(Not run)

Construct a density mixture from mixture weights and density functions.

Description

Construct a density mixture from mixture weights and density functions.

Usage

densityMixture(weights, densities)

epsilon process (residuals) of VAR model

Description

epsilon process (residuals) of VAR model

Usage

epsilon_var(zt, ar)

Fast Fourier Transform

Description

Fast Fourier Transform

Usage

fast_ft(x, real = T)

Details

If real: computes F_n X_n with the real-valued Fourier transformation matrix F_n (see Section 2.1 in Kirch et al (2018)). If !real: computes the complex-valued Fourier coefficients (see (4.5) in Meier (2018)).

Fast Inverse Fourier Transform

Description

Fast Inverse Fourier Transform

Usage

fast_ift(x, real = T, TOL = 1e-15)

Details

inverse function of fast_ft

Help function to compute the mean.

Description

Help function to compute the mean.

Usage

fast_mean(x)

Fourier frequencies

Description

Fourier frequencies on [0,pi], as defined by 2*pi*j/n for j=0,...,floor(n/2).

Usage

fourier_freq(n)

Arguments

Value

numeric vector of length floor(n/2)+1

Get epsilon process (i.e. model residuals) for ARMA(p,q)

Description

Get epsilon process (i.e. model residuals) for ARMA(p,q)

Usage

genEpsARMAC(zt, ar, ma)

Get U from phi, vectorized, cpp internal only

Description

Get U from phi, vectorized, cpp internal only

Usage

get_U_cpp(u_phi)

Construct psd mixture

Description

Construct psd mixture

Usage

get_f_matrix(U_phi, r, Z, k, dbList)

Gibbs sampler for Bayesian AR model in PACF parametrization, including support for TS to be a nuisance parameter

Description

Gibbs sampler for Bayesian AR model in PACF parametrization, including support for TS to be a nuisance parameter

Usage

gibbs_AR_nuisance_intern(
  data,
  mcmc_params,
  prior_params,
  model_params,
  full_lik = F
)

Gibbs sampling algorithm for VAR model

Description

Gibbs sampling algorithm for VAR model

Usage

gibbs_VAR_nuisance_intern(
  data,
  mcmc_params,
  prior_params,
  model_params,
  full_lik = F
)

Gibbs sampler for an autoregressive model with PACF parametrization.

Description

Obtain samples of the posterior of a Bayesian autoregressive model of fixed order.

Usage

gibbs_ar(
  data,
  ar.order,
  Ntotal,
  burnin,
  thin = 1,
  print_interval = 500,
  numerical_thresh = 1e-07,
  adaption.N = burnin,
  adaption.batchSize = 50,
  adaption.tar = 0.44,
  full_lik = F,
  rho.alpha = rep(1, ar.order),
  rho.beta = rep(1, ar.order),
  sigma2.alpha = 0.001,
  sigma2.beta = 0.001
)

Arguments

Details

Partial Autocorrelation Structure (PACF, uniform prior) and the residual variance sigma2 (inverse gamma prior) is used as model parametrization. The DIC is computed with two times the posterior variance of the deviance as effective number of parameters, see (7.10) in the referenced book by Gelman et al. Further details can be found in the simulation study section in the referenced paper by C. Kirch et al. For more information on the PACF parametrization, see the referenced paper by Barndorff-Nielsen and Schou.

Value

list containing the following fields:

References

C. Kirch et al. (2018) Beyond Whittle: Nonparametric Correction of a Parametric Likelihood With a Focus on Bayesian Time Series Analysis Bayesian Analysis <doi:10.1214/18-BA1126>

A. Gelman et al. (2013) Bayesian Data Analysis, Third Edition

O. Barndorff-Nielsen and G. Schou On the parametrization of autoregressive models by partial autocorrelations Journal of Multivariate Analysis (3),408-419 <doi:10.1016/0047-259X(73)90030-4>

Examples

## Not run: 

##
## Example 1: Fit an AR(p) model to sunspot data:
##

# Use this variable to set the AR model order
p <- 2

data <- sqrt(as.numeric(sunspot.year))
data <- data - mean(data)

# If you run the example be aware that this may take several minutes
print("example may take some time to run")
mcmc <- gibbs_ar(data=data, ar.order=p, Ntotal=10000, burnin=4000, thin=2)

# Plot spectral estimate, credible regions and periodogram on log-scale
plot(mcmc, log=T)


##
## Example 2: Fit an AR(p) model to high-peaked AR(1) data
##

# Use this variable to set the AR model order
p <- 1

n <- 256
data <- arima.sim(n=n, model=list(ar=0.95)) 
data <- data - mean(data)
omega <- fourier_freq(n)
psd_true <- psd_arma(omega, ar=0.95, ma=numeric(0), sigma2=1)

# If you run the example be aware that this may take several minutes
print("example may take some time to run")
mcmc <- gibbs_ar(data=data, ar.order=p, Ntotal=10000, burnin=4000, thin=2)

# Compare estimate with true function (green)
plot(mcmc, log=F, pdgrm=F, credib="uniform")
lines(x=omega, y=psd_true, col=3, lwd=2)

# Compute the Integrated Absolute Error (IAE) of posterior median
cat("IAE=", mean(abs(mcmc$psd.median-psd_true)[-1]) , sep="")

## End(Not run)

BDP-DW method: performing posterior sampling and calculating statistics based on the posterior samples

Description

BDP-DW method: performing posterior sampling and calculating statistics based on the posterior samples

Usage

gibbs_bdp_dw(
  data,
  m,
  likelihood_thinning = 1,
  monitor = TRUE,
  print_interval = 100,
  unif_CR = FALSE,
  rescaled_time,
  freq,
  Ntotal = 110000,
  burnin = 60000,
  thin = 10,
  adaptive.batchSize = 50,
  adaptive.targetAcceptanceRate = 0.44,
  M = 1,
  g0.alpha = 1,
  g0.beta = 1,
  k1.theta = 0.01,
  k2.theta = 0.01,
  tau.alpha = 0.001,
  tau.beta = 0.001,
  k1max = 100,
  k2max = 100,
  L = 10,
  bernstein1_l = 0.1,
  bernstein1_r = 0.9,
  bernstein2_l = 0.1,
  bernstein2_r = 0.9
)

Arguments

Value

list containing the following fields:

References

Tang et al. (2023) Bayesian nonparametric spectral analysis of locally stationary processes ArXiv preprint <arXiv:2303.11561>

Examples

## Not run: 

##
## Example: Applying BDP-DW method to a multi-peaked heavy-tailed tvMA(1) process
##

# set seed
set.seed(2)
# set the length of time series
len_d <- 1500
# generate data from DGP LS2b defined in Section 4.2.2 of Tang et al. (2023). 
# see also ?sim_tvarma12
sim_data <- sim_tvarma12(len_d = 1500, dgp = "LS2", innov_distribution = "b")
set.seed(NULL)
# specify grid-points at which the tv-PSD is evaluated
res_time <- seq(0, 1, by = 0.005); freq <- pi * seq(0, 1, by = 0.01)
# calculate the true tv-PSD of DGP LS2b at the pre-specified grid
true_tvPSD <- psd_tvarma12(rescaled_time = res_time, freq = freq, dgp = "LS2")
# plot the true tv-PSD
# type ?plot.bdp_dw_tv_psd for more info
plot(true_tvPSD)

# If you run the example be aware that this may take several minutes
print("This example may take some time to run")
result <- gibbs_bdp_dw(data = sim_data, 
m = 50, 
likelihood_thinning = 2, 
rescaled_time = res_time,
 freq = freq)

# extract bayes factor and examine posterior summary
bayes_factor(result)
summary(result)

# compare estimate with true function
# type ?plot.bdp_dw_result for more info
par(mfrow = c(1,2))

plot(result, which = 1,
zlim = range(result$tvpsd.mean, true_tvPSD$tv_psd)
)
plot(true_tvPSD,
zlim = range(result$tvpsd.mean, true_tvPSD$tv_psd),
main = "true tv-PSD")

par(mfrow = c(1,1))


## End(Not run)

Gibbs sampler for corrected parametric likelihood + Bernstein-Dirichlet mixture, including possibility of using time series as mere nuisance parameter Multivariate case

Description

Gibbs sampler for corrected parametric likelihood + Bernstein-Dirichlet mixture, including possibility of using time series as mere nuisance parameter Multivariate case

Usage

gibbs_multivariate_nuisance(
  data,
  mcmc_params,
  corrected,
  prior_params,
  model_params
)

Gibbs sampler in Cpp

Description

Gibbs sampler in Cpp

Usage

gibbs_multivariate_nuisance_cpp(
  data,
  NA_pos,
  FZ,
  eps_r,
  eps_Z,
  eps_U,
  k_0,
  r_0,
  Z_0,
  U_phi_0,
  phi_def,
  eta,
  omega,
  Sigma,
  Ntotal,
  print_interval,
  numerical_thresh,
  verbose,
  L,
  k_theta,
  dbList
)

Gibbs sampler for Bayesian nonparametric inference with Whittle likelihood

Description

Obtain samples of the posterior of the Whittle likelihood in conjunction with a Bernstein-Dirichlet prior on the spectral density.

Usage

gibbs_np(
  data,
  Ntotal,
  burnin,
  thin = 1,
  print_interval = 100,
  numerical_thresh = 1e-07,
  M = 1,
  g0.alpha = 1,
  g0.beta = 1,
  k.theta = 0.01,
  tau.alpha = 0.001,
  tau.beta = 0.001,
  kmax = 100 * coars + 500 * (!coars),
  trunc_l = 0.1,
  trunc_r = 0.9,
  coars = F,
  L = max(20, length(data)^(1/3))
)

Arguments

Details

Further details can be found in the simulation study section in the references papers.

Value

list containing the following fields:

References

C. Kirch et al. (2018) Beyond Whittle: Nonparametric Correction of a Parametric Likelihood With a Focus on Bayesian Time Series Analysis Bayesian Analysis <doi:10.1214/18-BA1126>

N. Choudhuri et al. (2004) Bayesian Estimation of the Spectral Density of a Time Series JASA <doi:10.1198/016214504000000557>

S. Ghosal and A. van der Vaart (2017) Fundamentals of Nonparametric Bayesian Inference <doi:10.1017/9781139029834>

Examples

## Not run: 

##
## Example 1: Fit the NP model to sunspot data:
##

data <- sqrt(as.numeric(sunspot.year))
data <- data - mean(data)

# If you run the example be aware that this may take several minutes
print("example may take some time to run")
mcmc <- gibbs_np(data=data, Ntotal=10000, burnin=4000, thin=2)

# Plot spectral estimate, credible regions and periodogram on log-scale
plot(mcmc, log=T)


##
## Example 2: Fit the NP model to high-peaked AR(1) data
##

n <- 256
data <- arima.sim(n=n, model=list(ar=0.95)) 
data <- data - mean(data)
omega <- fourier_freq(n)
psd_true <- psd_arma(omega, ar=0.95, ma=numeric(0), sigma2=1)

# If you run the example be aware that this may take several minutes
print("example may take some time to run")
mcmc <- gibbs_np(data=data, Ntotal=10000, burnin=4000, thin=2)

# Compare estimate with true function (green)
plot(mcmc, log=F, pdgrm=F, credib="uniform")
lines(x=omega, y=psd_true, col=3, lwd=2)

# Compute the Integrated Absolute Error (IAE) of posterior median
cat("IAE=", mean(abs(mcmc$psd.median-psd_true)[-1]) , sep="")

## End(Not run)

Gibbs sampler for Bayesian semiparametric inference with the corrected AR likelihood

Description

Obtain samples of the posterior of the corrected autoregressive likelihood in conjunction with a Bernstein-Dirichlet prior on the correction.

Usage

gibbs_npc(
  data,
  ar.order,
  Ntotal,
  burnin,
  thin = 1,
  print_interval = 100,
  numerical_thresh = 1e-07,
  adaption.N = burnin,
  adaption.batchSize = 50,
  adaption.tar = 0.44,
  full_lik = F,
  rho.alpha = rep(1, ar.order),
  rho.beta = rep(1, ar.order),
  eta = T,
  M = 1,
  g0.alpha = 1,
  g0.beta = 1,
  k.theta = 0.01,
  tau.alpha = 0.001,
  tau.beta = 0.001,
  trunc_l = 0.1,
  trunc_r = 0.9,
  coars = F,
  kmax = 100 * coars + 500 * (!coars),
  L = max(20, length(data)^(1/3))
)

Arguments

Details

Partial Autocorrelation Structure (PACF, uniform prior) and the residual variance sigma2 (inverse gamma prior) is used as model parametrization. A Bernstein-Dirichlet prior for c_eta with base measure Beta(g0.alpha, g0.beta) is used. Further details can be found in the simulation study section in the referenced paper by Kirch et al. For more information on the PACF parametrization, see the referenced paper by Barndorff-Nielsen and Schou.

Value

list containing the following fields:

References

C. Kirch et al. (2018) Beyond Whittle: Nonparametric Correction of a Parametric Likelihood With a Focus on Bayesian Time Series Analysis Bayesian Analysis <doi:10.1214/18-BA1126>

S. Ghosal and A. van der Vaart (2017) Fundamentals of Nonparametric Bayesian Inference <doi:10.1017/9781139029834>

O. Barndorff-Nielsen and G. Schou On the parametrization of autoregressive models by partial autocorrelations Journal of Multivariate Analysis (3),408-419 <doi:10.1016/0047-259X(73)90030-4>

Examples

## Not run: 

##
## Example 1: Fit a nonparametrically corrected AR(p) model to sunspot data:
##

# Use this variable to set the AR model order
p <- 2

data <- sqrt(as.numeric(sunspot.year))
data <- data - mean(data)

# If you run the example be aware that this may take several minutes
print("example may take some time to run")
mcmc <- gibbs_npc(data=data, ar.order=p, Ntotal=10000, burnin=4000, thin=2)

# Plot spectral estimate, credible regions and periodogram on log-scale
plot(mcmc, log=T)


##
## Example 2: Fit a nonparametrically corrected AR(p) model to high-peaked AR(1) data
##

# Use this variable to set the autoregressive model order
p <- 1

n <- 256
data <- arima.sim(n=n, model=list(ar=0.95)) 
data <- data - mean(data)
omega <- fourier_freq(n)
psd_true <- psd_arma(omega, ar=0.95, ma=numeric(0), sigma2=1)

# If you run the example be aware that this may take several minutes
print("example may take some time to run")
mcmc <- gibbs_npc(data=data, ar.order=p, Ntotal=10000, burnin=4000, thin=2)

# Compare estimate with true function (green)
plot(mcmc, log=F, pdgrm=F, credib="uniform")
lines(x=omega, y=psd_true, col=3, lwd=2)

# Compute the Integrated Absolute Error (IAE) of posterior median
cat("IAE=", mean(abs(mcmc$psd.median-psd_true)[-1]) , sep="")

## End(Not run)

Gibbs sampler for corrected parametric likelihood + Bernstein-Dirichlet mixture, including possibility of using time series as mere nuisance parameter

Description

Gibbs sampler for corrected parametric likelihood + Bernstein-Dirichlet mixture, including possibility of using time series as mere nuisance parameter

Usage

gibbs_nuisance(
  data,
  mcmc_params,
  corrected,
  prior_params,
  model_params,
  full_lik = NULL
)

Gibbs sampler for vector autoregressive model.

Description

Obtain samples of the posterior of a Bayesian VAR model of fixed order. An independent Normal-Inverse-Wishart prior is employed.

Usage

gibbs_var(
  data,
  ar.order,
  Ntotal,
  burnin,
  thin = 1,
  print_interval = 500,
  full_lik = F,
  beta.mu = rep(0, ar.order * ncol(data)^2),
  beta.Sigma = 10000 * diag(ar.order * ncol(data)^2),
  Sigma.S = 1e-04 * diag(ncol(data)),
  Sigma.nu = 1e-04
)

Arguments

Details

See Section 2.2.3 in Koop and Korobilis (2010) or Section 6.2 in Meier (2018) for further details

Value

list containing the following fields:

References

G. Koop and D. Korobilis (2010) Bayesian Multivariate Time Series Methods for Empirical Macroeconomics Foundations and Trends in Econometrics <doi:10.1561/0800000013>

A. Meier (2018) A Matrix Gamma Process and Applications to Bayesian Analysis of Multivariate Time Series PhD thesis, OvGU Magdeburg <https://opendata.uni-halle.de//handle/1981185920/13470>

Examples

## Not run: 

##
## Example 1: Fit a VAR(p) model to SOI/Recruitment series:
##

# Use this variable to set the VAR model order
p <- 5

data <- cbind(as.numeric(astsa::soi-mean(astsa::soi)), 
              as.numeric(astsa::rec-mean(astsa::rec)) / 50)
data <- apply(data, 2, function(x) x-mean(x))

# If you run the example be aware that this may take several minutes
print("example may take some time to run")
mcmc <- gibbs_var(data=data, ar.order=p, Ntotal=10000, burnin=4000, thin=2)

# Plot spectral estimate, credible regions and periodogram on log-scale
plot(mcmc, log=T)



##
## Example 2: Fit a VAR(p) model to VMA(1) data
##

# Use this variable to set the VAR model order
p <- 5

n <- 256
ma <- rbind(c(-0.75, 0.5), c(0.5, 0.75))
Sigma <- rbind(c(1, 0.5), c(0.5, 1))
data <- sim_varma(model=list(ma=ma), n=n, d=2)
data <- apply(data, 2, function(x) x-mean(x))

# If you run the example be aware that this may take several minutes
print("example may take some time to run")
mcmc <- gibbs_var(data=data, ar.order=p, Ntotal=10000, burnin=4000, thin=2)

# Plot spectral estimate, credible regions and periodogram on log-scale
plot(mcmc, log=T)

## End(Not run)

Gibbs sampler for multivaiate Bayesian nonparametric inference with Whittle likelihood

Description

Obtain samples of the posterior of the multivariate Whittle likelihood in conjunction with an Hpd AGamma process prior on the spectral density matrix.

Usage

gibbs_vnp(
  data,
  Ntotal,
  burnin,
  thin = 1,
  print_interval = 100,
  numerical_thresh = 1e-07,
  adaption.N = burnin,
  adaption.batchSize = 50,
  adaption.tar = 0.44,
  eta = ncol(data),
  omega = ncol(data),
  Sigma = 10000 * diag(ncol(data)),
  k.theta = 0.01,
  kmax = 100 * coars + 500 * (!coars),
  trunc_l = 0.1,
  trunc_r = 0.9,
  coars = F,
  L = max(20, length(data)^(1/3))
)

Arguments

Details

A detailed description of the method can be found in Section 5 in Meier (2018).

Value

list containing the following fields:

References

A. Meier (2018) A Matrix Gamma Process and Applications to Bayesian Analysis of Multivariate Time Series PhD thesis, OvGU Magdeburg <https://opendata.uni-halle.de//handle/1981185920/13470>

Examples

## Not run: 

##
## Example: Fit multivariate NP model to SOI/Recruitment series:
##

data <- cbind(as.numeric(astsa::soi-mean(astsa::soi)), 
              as.numeric(astsa::rec-mean(astsa::rec)) / 50)
data <- apply(data, 2, function(x) x-mean(x))

# If you run the example be aware that this may take several minutes
print("example may take some time to run")
mcmc <- gibbs_vnp(data=data, Ntotal=10000, burnin=4000, thin=2)

# Visualize results
plot(mcmc, log=T)


##
## Example 2: Fit multivariate NP model to VMA(1) data
##

n <- 256
ma <- rbind(c(-0.75, 0.5), c(0.5, 0.75))
Sigma <- rbind(c(1, 0.5), c(0.5, 1))
data <- sim_varma(model=list(ma=ma), n=n, d=2)
data <- apply(data, 2, function(x) x-mean(x))

# If you run the example be aware that this may take several minutes
print("example may take some time to run")
mcmc <- gibbs_vnp(data=data, Ntotal=10000, burnin=4000, thin=2)

# Plot spectral estimate, credible regions and periodogram on log-scale
plot(mcmc, log=T)

## End(Not run)

Does a matrix have an eigenvalue smaller than 0?

Description

Does a matrix have an eigenvalue smaller than 0?

Usage

hasEigenValueSmallerZero(A, TOL = 0)

Check if a matrix is Hermitian positive definite

Description

Check if a matrix is Hermitian positive definite

Usage

is_hpd(A, tol = 1e-15)

Is l quadratic?

Description

Is l quadratic?

Usage

is_quadratic(l, thresh = 1e-15)

Check if a matrix is symmetric positive definite

Description

Check if a matrix is symmetric positive definite

Usage

is_spd(A, tol = 1e-05)

Likelihood of an autoregressive time series model with i.i.d. normal innovations

Description

Likelihood of an autoregressive time series model with i.i.d. normal innovations

Usage

lik_ar(x, ar, mean = 0, sd = 1, log = F, full = T)

Arguments

Value

numeric value for the likelihood or log-likelihood

Log corrected parametric AR likelihood (Gaussian)

Description

Log corrected parametric AR likelihood (Gaussian)

Usage

llike(
  omega,
  FZ,
  ar,
  ma,
  v,
  w,
  k,
  tau,
  corrected,
  toggle.q,
  pdgrm,
  beta_basis_k,
  nll_fun,
  f.alpha,
  excludeBoundary,
  full_lik
)

Details

See (5) in Kirch et al (2018)

Time domain AR(p) likelihood for nuisance/noise time series

Description

Time domain AR(p) likelihood for nuisance/noise time series

Usage

llike_AR(x_noise, rho, sigma2, full_lik)

Calculating log likelihood

Description

Calculating log likelihood

Usage

llike_dw(FZ, npsd, tau)

VAR(p) likelihood

Description

VAR(p) likelihood

Usage

llike_var(zt, ar, sigma, full_lik)

VAR(p) full likelihood

Description

VAR(p) full likelihood

Usage

llike_var_full(zt, ar, sigma)

VAR(p) partial likelihood (unnormalized) Note: Fine for fixed p, but not suited for model comparison

Description

VAR(p) partial likelihood (unnormalized) Note: Fine for fixed p, but not suited for model comparison

Usage

llike_var_partial(zt, ar, sigma)

Calculate the moving Fourier transform ordinates

Description

Calculate the moving Fourier transform ordinates

Usage

local_moving_FT_zigzag(x, m, thinning_factor)

Arguments

Value

A list containing the moving Fourier transform and corresponding time grid.

References

Y. Tang et al. (2023) Bayesian nonparametric spectral analysis of locally stationary processes ArXiv preprint <arXiv:2303.11561>

Examples

set.seed(1); x <- rnorm(1500)
local_moving_FT_zigzag(x, 50, 1)

Log determinant of stick breaking transformation V -> p

Description

Log determinant of stick breaking transformation V -> p

Usage

logDet_stickBreaking(v)

Details

See Section 3 in Choudhuri et al (2004)

Fuller Logarithm

Description

Fuller Logarithm

Usage

logfuller(x, xi = 1e-08)

Details

see Fuller (1996), page 496

References

W. Fuller (1996) Introduction to Statistical Time Series Wiley Series in Probability and Statistics

Log posterior = log prior + log corrected parametric likelihood

Description

Log posterior = log prior + log corrected parametric likelihood

Usage

lpost(
  omega,
  FZ,
  ar,
  ma,
  v,
  w,
  k,
  tau,
  M,
  g0.alpha,
  g0.beta,
  k.theta,
  tau.alpha,
  tau.beta,
  corrected,
  toggle.q,
  pdgrm,
  beta_basis_k,
  nll_fun,
  f.alpha,
  rho,
  rho.alpha,
  rho.beta,
  excludeBoundary,
  full_lik
)

Log Posterior = Log Prior + (conditional) Log Likelihood

Description

Log Posterior = Log Prior + (conditional) Log Likelihood

Usage

lpost_AR(
  x_noise,
  rho,
  rho.alpha,
  rho.beta,
  sigma2,
  sigma2.alpha,
  sigma2.beta,
  full_lik
)

Log prior of Bernstein-Dirichlet mixture and parametric working model – all unnormalized

Description

Log prior of Bernstein-Dirichlet mixture and parametric working model – all unnormalized

Usage

lprior(
  v,
  w,
  k,
  tau,
  M,
  g0.alpha,
  g0.beta,
  k.theta,
  tau.alpha,
  tau.beta,
  f.alpha,
  corrected,
  rho,
  rho.alpha,
  rho.beta
)

Details

See Section 3 in Kirch et al (2018). Hyperparameters are M, g0.a, g0.b, k.theta, tau.alpha, tau.beta. Note: Flat prior on f.alpha.

Log prior for PACF (~Beta) and sigma2 (~InverseGamma), unnormalized

Description

Log prior for PACF (~Beta) and sigma2 (~InverseGamma), unnormalized

Usage

lprior_AR(rho, rho.alpha, rho.beta, sigma2, sigma2.alpha, sigma2.beta)

Calculation of log prior

Description

Calculation of log prior

Usage

lprior_dw(
  tilde.E,
  w1,
  w2,
  k1,
  k2,
  tau,
  g0.alpha,
  g0.beta,
  k1.theta,
  k2.theta,
  tau.alpha,
  tau.beta
)

Multivariate discrete (fast) Fourier Transform

Description

Multivariate discrete (fast) Fourier Transform

Usage

mdft(Z, real = F)

Multivariate inverse discrete (fast) Fourier Transform

Description

Multivariate inverse discrete (fast) Fourier Transform

Usage

midft(FZ, real = F)

Get string representation for missing values position from vector index

Description

Get string representation for missing values position from vector index

Usage

missingValues_str_help(ind, n)

Get mixture weights of Bernstein-Dirchlet-Mixtures

Description

Get mixture weights of Bernstein-Dirchlet-Mixtures

Usage

mixtureWeight(p, w, k)

Compute Periodgram matrix from (complex-valued) Fourier coefficients

Description

Compute Periodgram matrix from (complex-valued) Fourier coefficients

Usage

mpdgrm(FZ)

Details

see (4.7) in Meier (2018)

Generate a random samples from a Dirichlet distribution

Description

Generate a random samples from a Dirichlet distribution

Usage

my_rdirichlet(alpha)

Arguments

Value

a vector of the same length as alpha

Negative log likelihood of iid standard normal observations [unit variance] Note: deprecated

Description

Negative log likelihood of iid standard normal observations [unit variance] Note: deprecated

Usage

nll_norm(epsilon_t, ...)

Fourier frequencies rescaled on the unit interval

Description

Fourier frequencies rescaled on the unit interval

Usage

omegaFreq(n)

Get p from v in Stick Breaking DP representation

Description

Get p from v in Stick Breaking DP representation

Usage

pFromV(v)

C++ function for computing AR coefficients from PACF. See Section III in Barndorff-Nielsen and Schou (1973)

Description

C++ function for computing AR coefficients from PACF. See Section III in Barndorff-Nielsen and Schou (1973)

Usage

pacf2AR(pacf)

References

O. Barndorff-Nielsen and G. Schou (1973) On the Parametrization of Autoregressive Models by Partial Autocorrelations Journal of Multivariate Analysis (3),408-419 <doi:10.1016/0047-259X(73)90030-4>

Convert partial autocorrelation coefficients to AR coefficients.

Description

Convert partial autocorrelation coefficients to AR coefficients.

Usage

pacf_to_ar(pacf)

Arguments

Details

See Section 2 in Kirch et al (2018) or Section III in Barndorff-Nielsen and Schou (1973) for further details

Value

numeric vector of autoregressive model coefficients

References

C. Kirch et al Supplemental material of Beyond Whittle: Nonparametric Correction of a Parametric Likelihood With a Focus on Bayesian Time Series Analysis Bayesian Analysis <doi:10.1214/18-BA1126SUPP>

O. Barndorff-Nielsen and G. Schou On the parametrization of autoregressive models by partial autocorrelations Journal of Multivariate Analysis (3),408-419 <doi:10.1016/0047-259X(73)90030-4>

See Also

acf2AR, ARMAacf

Convert vector parametrization (beta) to matrix-parametrization (phi), the latter as e.g. used in MTS::VAR()$ar

Description

Convert vector parametrization (beta) to matrix-parametrization (phi), the latter as e.g. used in MTS::VAR()$ar

Usage

phiFromBeta_normalInverseWishart(beta, K, p)

Arguments

Value

K times K*p coefficient matrix

Plot method for bdp_dw_result class

Description

Plot method for bdp_dw_result class

Usage

## S3 method for class 'bdp_dw_result'
plot(
  x,
  which = 1:4,
  ask = prod(par("mfcol")) < length(which) && dev.interactive(),
  col = hcl.colors(200, "Blue-Red 3"),
  ...
)

Arguments

Plot method for bdp_dw_tv_psd class

Description

Plot method for bdp_dw_tv_psd class

Usage

## S3 method for class 'bdp_dw_tv_psd'
plot(x, col = hcl.colors(200, "Blue-Red 3"), ...)

Arguments

Details

Visualizes the spectral density function of time-varying

Plot method for gibbs_psd class

Description

Plot method for gibbs_psd class

Usage

## S3 method for class 'gibbs_psd'
plot(x, pdgrm = T, credib = "both", log = T, ...)

Arguments

Details

Visualizes the spectral density estimate (pointwise posterior median), along with the periodogram and credibility regions. If the data has missing values, the periodogram is computed with a linearly interpolated version of the data using na.interp.

Visualization of multivariate PSDs Used in plot.gibbs_psd

Description

Visualization of multivariate PSDs Used in plot.gibbs_psd

Usage

plotMPsd(
  f,
  g = NULL,
  lty = rep(1, 1 + length(g)),
  col = rep(1, 1 + length(g)),
  type = rep("l", 1 + length(g)),
  pch = rep("0", 1 + length(g)),
  lwd = rep(1, 1 + length(g)),
  log = F,
  ylim.compound = T,
  mar = c(4, 4, 4, 3),
  mfrow = c(dim(f)[1], dim(f)[1]),
  lambda_scaling = pi,
  ylab_prefix = "f",
  xlab = "Frequency",
  excludeBoundary = T,
  ...
)

Print method for bdp_dw_result class

Description

Print method for bdp_dw_result class

Usage

## S3 method for class 'bdp_dw_result'
print(x, ...)

Arguments

Print method for gibbs_psd class

Description

Print method for gibbs_psd class

Usage

## S3 method for class 'gibbs_psd'
print(x, ...)

Arguments

Help function to print MCMC state

Description

Help function to print MCMC state

Usage

print_mcmc_state(i, Ntotal, tim, tim0)

Helping function for print and summary (both are quite similar)

Description

Helping function for print and summary (both are quite similar)

Usage

print_summary_gibbs_psd_help(mcmc, flag = T)

Arguments

Help function to print debugging messages

Description

Help function to print debugging messages

Usage

print_warn(msg)

ARMA(p,q) spectral density function

Description

Evaluate the ARMA(p,q) spectral density at some frequencies freq in [0,pi), Note that no test for model stationarity is performed.

Usage

psd_arma(freq, ar, ma, sigma2 = 1)

Arguments

Details

See section 4.4 in the referenced book

Value

numeric vector of the (real-valued) spectral density values

References

P. J. Brockwell and R. Davis (1996) Time Series: Theory and Methods (Second Edition)

Convert psd vector to array (compatibility: to use plotMPsd for univariate functions as well)

Description

Convert psd vector to array (compatibility: to use plotMPsd for univariate functions as well)

Usage

psd_array(f)

Time series model X_t=e_t, E[e_t]=0

Description

Time series model X_t=e_t, E[e_t]=0

Usage

psd_dummy_model()

time-varying spectral density function of the tvARMA(1,2) processes for illustrations

Description

time-varying spectral density function of the tvARMA(1,2) processes for illustrations

Usage

psd_tvarma12(
  rescaled_time,
  freq,
  dgp = NULL,
  a1 = function(u) {
     rep(0, length(u))
 },
  b1 = function(u) {
     rep(0, length(u))
 },
  b2 = function(u) {
     rep(0, length(u))
 }
)

Arguments

Details

See sim_tvarma12 for the precise definition of a tvARMA(1,2) process. The time-varying spectral density function of this process is defined as

where u is called rescaled time and \lambda is called frequency.

For dgp = "LS1", it is a tvMA(2) process (MA order is 2) with

For dgp = "LS2", it is a tvMA(1) process (MA order is 1) with

For dgp = "LS3", it is a tvAR(1) process (MA order is 0) with

Value

a matrix of dimension length(rescaled_time) by length(freq).

References

Tang et al. (2023) Bayesian nonparametric spectral analysis of locally stationary processes ArXiv preprint <arXiv:2303.11561>

Examples

## Not run: 
res_time <- seq(0, 1, by = 0.005); freq <- pi*seq(0, 1, by = 0.01)
true_tvPSD <- psd_tvarma12(rescaled_time = res_time, freq = freq, dgp = "LS2")
plot(true_tvPSD)

## End(Not run)

VARMA(p,q) spectral density function

Description

Evaluate the VARMA(p,q) spectral density at some frequencies freq in [0,pi). Note that no test for model stationarity is performed.

Usage

psd_varma(
  freq,
  ar = matrix(nrow = nrow(Sigma), ncol = 0),
  ma = matrix(nrow = nrow(Sigma), ncol = 0),
  Sigma
)

Arguments

Details

See section 11.5 in the referenced book

Value

an array containing the values of the varma psd matrix at freq

References

P. J. Brockwell and R. Davis (1996) Time Series: Theory and Methods (Second Edition)

helping function for psd_varma

Description

helping function for psd_varma

Usage

psd_varma_help(
  freq,
  ar = matrix(nrow = nrow(sigma), ncol = 0),
  ma = matrix(nrow = nrow(sigma), ncol = 0),
  sigma
)

Compute normalized PSD in the Bernstein-Dirichlet parametrization.

Description

Compute normalized PSD in the Bernstein-Dirichlet parametrization.

Usage

qpsd(omega, v, w, k, beta_basis_k, epsilon = 1e-20)

Details

See (5) in Choudhuri et al (2004)

Evaluation of normalized time-varying spectral density function (based on posterior samples)

Description

Evaluation of normalized time-varying spectral density function (based on posterior samples)

Usage

qpsd_dw(v, w1, w2, k1, k2, beta_basis_1_k, beta_basis_2_k)

Evaluation of normalized time-varying spectral density function (for MCMC algorithm)

Description

Evaluation of normalized time-varying spectral density function (for MCMC algorithm)

Usage

qpsd_dw.tilde_zigzag_cpp_expedited(
  tilde.v,
  tilde.w1,
  tilde.w2,
  k1,
  k2,
  beta_basis_1_k,
  beta_basis_2_k
)

Store imaginary parts above and real parts below the diagonal

Description

Store imaginary parts above and real parts below the diagonal

Usage

realValuedPsd(f_)

Simulate from a Multivariate Normal Distribution

Description

Produces one or more samples from the specified multivariate normal distribution.

Usage

rmvnorm(n, d, mu = rep(0, d), Sigma = diag(d), ...)

Arguments

Details

This is a simple wrapper function based on mvrnorm, to be used within sim_varma

Value

If n=1 a vector of length d, otherwise an n by d matrix with one sample in each row.

Negative log AR likelihood values for scree-type plots

Description

(Approximate) negative maximum log-likelihood for for different autoregressive orders to produce scree-type plots.

Usage

scree_type_ar(data, order.max, method = "yw")

Arguments

Details

By default, the maximum likelihood is approximated by the Yule-Walker method, due to numerical stabililty and computational speed. Further details can be found in the simulation study section in the referenced paper.

Value

a data frame containing the autoregressive orders p and the corresponding negative log likelihood values nll

References

C. Kirch et al. (2018) Beyond Whittle: Nonparametric Correction of a Parametric Likelihood With a Focus on Bayesian Time Series Analysis Bayesian Analysis <doi:10.1214/18-BA1126>

Examples

## Not run: 

###
### Interactive visual inspection for the sunspot data
###

data <- sqrt(as.numeric(sunspot.year))
data <- data <- data - mean(data)

screeType <- scree_type_ar(data, order.max=15)

# Determine the autoregressive order by an interactive visual inspection of the scree-type plot
plot(x=screeType$p, y=screeType$nll, type="b")
p_ind <- identify(x=screeType$p, y=screeType$nll, n=1, labels=screeType$p)
print(screeType$p[p_ind])

## End(Not run)

simulate from the tvARMA(1,2) process for illustration

Description

simulate from the tvARMA(1,2) process for illustration

Usage

sim_tvarma12(
  len_d,
  dgp = NULL,
  ar_order = 1,
  ma_order = 2,
  a1 = NULL,
  b1 = NULL,
  b2 = NULL,
  innov_distribution = NULL,
  wn = NULL
)

Arguments

Details

This function simulates from the following time-varying Autoregressive Moving Average model with order (1,2):

where T is the length specified and {w_t} are a sequence of i.i.d. random variables with mean 0 and standard deviation 1.

For dgp = "LS1", it is a tvMA(2) process (MA order is 2) with

For dgp = "LS2", it is a tvMA(1) process (MA order is 1) with

For dgp = "LS3", it is a tvAR(1) process (MA order is 0) with

Value

a numeric vector of length len_d simulated from the given process.

References

Tang et al. (2023) Bayesian nonparametric spectral analysis of locally stationary processes ArXiv preprint <arXiv:2303.11561>

Examples

## Not run: 
sim_tvarma12(len_d = 1500, 
dgp = "LS2", 
innov_distribution = "a") # generate from LS2a

sim_tvarma12(len_d = 1500, 
dgp = "LS2", 
wn = rnorm(1502)) # again, generate from LS2a

sim_tvarma12(len_d = 1500, 
ar_order = 0, 
ma_order = 1, 
b1 = function(u){1.1*cos(1.5 - cos(4*pi*u))}, 
innov_distribution = "a") # again, generate from LS2a

sim_tvarma12(len_d = 1500, 
ar_order = 0, 
ma_order = 1, 
b1 = function(u){1.1*cos(1.5 - cos(4*pi*u))}, 
wn = rnorm(1502)) # again, generate from LS2a

## End(Not run)

Simulate from a VARMA model

Description

Simulate from a Vector Autoregressive Moving Average (VARMA) model. Note that no test for model stationarity is performed.

Usage

sim_varma(model, n, d, rand.gen = rmvnorm, burnin = 10000, ...)

Arguments

Value

If n=1 a vector of length d, otherwise an n by d matrix with one sample in each row.

See Also

arima.sim to simulate from univariate ARMA models

Examples

## Not run: 
# Example: Draw from bivariate normal VAR(2) model
ar <- rbind(c(.5, 0, 0, 0), c(0, -.3, 0, -.5))
Sigma <- matrix(data=c(1, .9, .9, 1), nrow=2, ncol=2)
x <- sim_varma(n=256, d=2, model=list(ar=ar))
plot.ts(x)

## End(Not run)

sum of multivariate normal log densities with mean 0 and covariance Sigma, unnormalized

Description

sum of multivariate normal log densities with mean 0 and covariance Sigma, unnormalized

Usage

sldmvnorm(z_t, Sigma)

Summary method for bdp_dw_result class

Description

Summary method for bdp_dw_result class

Usage

## S3 method for class 'bdp_dw_result'
summary(object, ...)

Arguments

Summary method for gibbs_psd class

Description

Summary method for gibbs_psd class

Usage

## S3 method for class 'gibbs_psd'
summary(object, ...)

Arguments

VARMA transfer polynomials

Description

VARMA transfer polynomials

Usage

transfer_polynomial(lambda, coef)

Helping function for uci_matrix

Description

Helping function for uci_matrix

Usage

uci_help(fpsd.sample, alpha, log = F)

Uniform credible intervals in matrix-valued case

Description

Uniform credible intervals in matrix-valued case

Usage

uci_matrix(fpsd.sample, alpha, uniform_among_components = F)

Details

see (6.5) in Meier (2018)

Uniform maximum, as needed for uniform credible intervals

Description

Uniform maximum, as needed for uniform credible intervals

Usage

uniformmax(sample)

Details

see Section 4.1 in Kirch et al (2018)

Helping function for uci_matrix

Description

Helping function for uci_matrix

Usage

uniformmax_help(sample)

Helping function for uci_matrix

Description

Helping function for uci_matrix

Usage

uniformmax_multi(mSample)

Range intervals I_l, see (63) in Mittelbach et al.

Description

Range intervals I_l, see (63) in Mittelbach et al.

Usage

unit_trace_I_l(l)

Get L (lower triangular Cholesky) from x Called U^* in Mittelbach et al, see (60) there

Description

Get L (lower triangular Cholesky) from x Called U^* in Mittelbach et al, see (60) there

Usage

unit_trace_L_from_x(x)

Get U (Hpd with unit trace) matrix from phi (hyperspherical coordinates) vector.

Description

Get U (Hpd with unit trace) matrix from phi (hyperspherical coordinates) vector.

Usage

unit_trace_U_from_phi(phi)

Get log(c) vector, see (70) in Mittelbach et al. Helping function for unit_trace_runif

Description

Get log(c) vector, see (70) in Mittelbach et al. Helping function for unit_trace_runif

Usage

unit_trace_log_c(p, q)

Get log(d) vector, see (39) in Mittelbach et al, adjusted to complex case Helping function for unit_trace_runif

Description

Get log(d) vector, see (39) in Mittelbach et al, adjusted to complex case Helping function for unit_trace_runif

Usage

unit_trace_log_d(p, q)

Get log(f_l), see (66) in Mittelbach et al. Helping function for unit_trace_runif

Description

Get log(f_l), see (66) in Mittelbach et al. Helping function for unit_trace_runif

Usage

unit_trace_log_f_l(phi, p, q, log_c, l)

Get mu vector, see (36) in Mittelbach et al. Helping function for unit_trace_runif

Description

Get mu vector, see (36) in Mittelbach et al. Helping function for unit_trace_runif

Usage

unit_trace_mu(p, q)

Get log(nu) vector, see (38) in Mittelbach et al. Helping function for unit_trace_runif

Description

Get log(nu) vector, see (38) in Mittelbach et al. Helping function for unit_trace_runif

Usage

unit_trace_nu(sigma2_vec, log_c_vec, log_d_vec)

Get p vector, see (67) in Mittelbach et al.

Description

Get p vector, see (67) in Mittelbach et al.

Usage

unit_trace_p(d)

Get q vector, see (68) in Mittelbach et al.

Description

Get q vector, see (68) in Mittelbach et al.

Usage

unit_trace_q(d)

Draw uniformly from Hpd matrices with unit trace

Description

Draw uniformly from Hpd matrices with unit trace

Usage

unit_trace_runif(n, d, verbose = F)

Obtain one uniform draw from d times d Hpd matrices with unit trace See Algorithm 2 in Mittelbach et al. (adjusted to complex case)

Description

Obtain one uniform draw from d times d Hpd matrices with unit trace See Algorithm 2 in Mittelbach et al. (adjusted to complex case)

Usage

unit_trace_runif_single(d)

Get sigma2 vector, see (70) in Mittelbach et al. Helping function for unit_trace_runif

Description

Get sigma2 vector, see (70) in Mittelbach et al. Helping function for unit_trace_runif

Usage

unit_trace_sigma2(p, q)

Get x from phi, see (62) in Mittelbach et al.

Description

Get x from phi, see (62) in Mittelbach et al.

Usage

unit_trace_x_from_phi(phi)

Redundantly roll out a PSD from length N=floor(n/2) to length n

Description

Redundantly roll out a PSD from length N=floor(n/2) to length n

Usage

unrollPsd(qPsd, n)

Get v from p (DP inverse stick breaking) Note: p is assumed to have length L, i.e. it does NOT contain p_0

Description

Get v from p (DP inverse stick breaking) Note: p is assumed to have length L, i.e. it does NOT contain p_0

Usage

vFromP(p, eps = 1e-08)

Get VARMA PSD from transfer polynomials Helping function for psd_varma

Description

Get VARMA PSD from transfer polynomials Helping function for psd_varma

Usage

varma_transfer2psd(transfer_ar_, transfer_ma_, sigma)