Package 'CovEsts'

Compute the Kernel Regression Estimator.

Description

This function computes the kernel regression estimator of the autocovariance function.

Usage

adjusted_est(
  X,
  x,
  t,
  b,
  kernel_name = c("gaussian", "wave", "rational_quadratic", "bessel_j"),
  kernel_params = c(),
  pd = TRUE,
  type = c("autocovariance", "autocorrelation"),
  meanX = mean(X),
  custom_kernel = FALSE,
  parallel = FALSE,
  ncores = parallel::detectCores() - 1,
  cl_export = NULL,
  cl = NULL
)

Arguments

X	A vector representing observed values of the time series.
x	A vector of lags.
t	The arguments at which the autocovariance function is calculated at.
b	Bandwidth parameter, greater than 0.
kernel_name	The name of the symmetric kernel (see kernel_symm_ec) function to be used. Possible values are: gaussian, wave, rational_quadratic, and bessel_j. Alternatively, a custom kernel function can be provided, see the examples.
kernel_params	A vector of parameters of the kernel function. See kernel_symm_ec for parameters.
pd	Whether a positive-definite estimate should be used. Defaults to TRUE.
type	Compute either the 'autocovariance' or 'autocorrelation'. Defaults to 'autocovariance'.
meanX	The average value of X. Defaults to mean(X).
custom_kernel	If a custom kernel is to be used or not. Defaults to FALSE.
parallel	Whether or not the computations should be done in parallel or not. Defaults to FALSE.
ncores	The number of cores to be used in the parallel computations. Defaults to the number cores - 1 (threads if hyperthreading is available), calculated from parallel::detectCores() - 1.
cl_export	A vector of any additional functions or variables to export for parallel computations. This may be required if estimator is not within the package. Defaults to NULL.
cl	An optional cluster object created by parallel::makeCluster. Defaults to NULL, which creates a temporary PSOCK cluster.

Details

The kernel regression estimator of an autocovariance function is defined as

$\hat{\rho}(t) = \left( \sum_{i=1}^{N} \sum_{j=1}^{N} \check{X}_{ij} K((t - (t_{i} - t_{j})) / b) \right) \left( \sum_{i=1}^{N} \sum_{j=1}^{N} K((t - (t_{i} - t_{j})) / b) \right)^{-1},$

where $\check{X}_{ij} = (X(t_{i}) - \bar{X}) (X(t_{j}) - \bar{X}).$

If pd is TRUE, the estimator will be made positive-definite through the following procedure

1. Take the discrete Fourier cosine transform, $\widehat{\mathcal{F}}(\theta)$, of the estimated autocovariance function

2. Compute a modified spectrum $\widetilde{\mathcal{F}}(\theta) = \max(\widehat{\mathcal{F}}(\theta), 0)$ for all sample frequencies.

3. Perform the Fourier inversion to obtain a new estimator.

Value

A CovEsts S3 object (list) with the following values

acf

A numeric vector containing the autocovariance/autocorrelation estimates.

lags

A numeric vector containing the lag indices used to compute the estimates on.

est_type

The type of estimate, namely 'autocorrelation' or 'autocovariance', this depends on the type parameter.

est_used

The estimator function used, in this case, 'adjusted_est'.

References

Hall, P. & Patil, P. (1994). Properties of nonparametric estimators of autocovariance for stationary random fields. Probability Theory and Related Fields 99(3), 399-424. https://doi.org/10.1007/bf01199899

Hall, P., Fisher, N. I., & Hoffmann, B. (1994). On the nonparametric estimation of covariance functions. The Annals of Statistics 22(4), 2115-2134. https://doi.org/10.1214/aos/1176325774

Examples

X <- c(1, 2, 3, 4)
adjusted_est(X, 1:4, 1:3, 0.1, "gaussian")
my_kernel <- function(x, params) {
  stopifnot(params[1] > 0, length(x) >= 1)
  return(exp(-((abs(x) / params[1])^params[2])) * (2 * params[1] * gamma(1 + 1/params[2])))
}
adjusted_est(X, 1:4, 1:3, 0.1, my_kernel, c(0.25), custom_kernel = TRUE)
## Not run: 
library(parallel)
X <- c(1, 2, 3, 4)
my_cl <- makePSOCKcluster(2)
adjusted_est(X, 1:4, 1:3, 0.1, "gaussian", parallel = TRUE, cl = my_cl)
stopCluster(my_cl)

## End(Not run)

---

Compute Adjusted Splines.

Description

A helper function that is an implementation of the formula from Choi, Li & Wang (2013, p. 616),

$f_{j}^{(l)}(x) = \frac{m + 1}{l} \left( f_{j}^{(l - 1)}(x + 1) - \tau_{j - p} f_{j}^{(l - 1)}(x) + \tau_{j - p + l + 1} f_{j + 1}^{(l - 1)}(x) - f_{j + 1}^{(l - 1)}(x + 1) \right) ,$

where $m$ is the number of nonboundary knots, $p$ is the order of the spline, $l$ is the order of the adjusted spline (the function $f_{j}^{(l)}(\cdot)$) and $j = 1, 2, \dots , m + p.$

Usage

adjusted_spline(x, j, l, p, m, taus)

Arguments

x	Argument of the function.
j	Index of basis function of order $l$.
l	Order of function.
p	The order of the splines.
m	The number of nonboundary knots.
taus	Vector of $\tau$s, see get_taus.

Value

A numeric value of the adjusted spline $f_{j}^{(l)}(x).$

References

Choi, I., Li, B. & Wang, X. (2013). Nonparametric Estimation of Spatial and Space-Time Covariance Function. JABES 18, 611-630. https://doi.org/10.1007/s13253-013-0152-z

Examples

## Not run: 
taus <- get_taus(3, 2)
adjusted_spline(1, 2, 1, 3, 2, taus)

## End(Not run)

---

Area Between Estimated Autocovariance Functions.

Description

This function estimates the area between two estimated autocovariance functions.

Usage

area_between(est1, est2, lags = c(), plot = FALSE)

## S3 method for class 'CovEsts'
area_between(est1, est2, lags = c(), plot = FALSE)

## Default S3 method:
area_between(est1, est2, lags = c(), plot = FALSE)

Arguments

est1	A numeric vector representing the first estimated autocovariance function.
est2	A numeric vector of the same length as est1 representing the second estimated autocovariance function
lags	An optional vector of lags starting from 0 up until some other lag. If empty, a vector of lags is created starting from 0 until len(est1) - 1, by 1.
plot	A boolean determining whether a plot should be created. By default, no plot is created.

Details

This function estimates the area between two estimated autocovariance functions over a set of lags, from 0 up to $h_{n}$ defined by

$\int_{0}^{h_{n}} \left| \hat{C}_{1}(h) - \hat{C}_{2}(h) \right| dh ,$

where $\hat{C}_{1}(\cdot)$ and $\hat{C}_{2}(\cdot)$ are estimated autocovariance functions.

To approximate this integral the trapezoidal rule is used.

If lags is empty, a uniform time grid with a step of 1 will be used which may result in a different area than if lags is specified.

Value

A numeric value representing the estimated area between two estimated autocovariance functions.

Methods (by class)

• area_between(CovEsts): Method for CovEsts objects.

• area_between(default): Method for numeric vectors.

Examples

x <- seq(0, 5, by=0.1)
estCov1 <- exp(-x^2)
estCov2 <- exp(-x^2.1)
area_between(estCov1, estCov2, lags=x)
area_between(estCov1, estCov2, lags=x, plot = TRUE)

---

as.double Method for CovEsts Objects

Description

This function provides a method for as.double for CovEsts objects, allowing for as.numeric to be called. It extracts the acf values from the object.

Usage

## S3 method for class 'CovEsts'
as.double(x, ...)

Arguments

x	A CovEsts S3 object.
...	Refer to base::double.

Value

The acf values from the CovEsts object.

Examples

as.numeric(standard_est(c(1, 2, 3)))

---

Block Bootstrap

Description

This function performs block bootstrap (moving or circular) to obtain a $(1-\alpha)\%$ confidence-interval for the autocovariance function. It will also compute average autocovariance function across all bootstrapped estimates.

Usage

block_bootstrap(
  X,
  maxLag,
  x = 0:length(X),
  n_bootstrap = 100,
  l = ceiling(length(X)^(1/3)),
  estimator = standard_est,
  type = c("autocovariance", "autocorrelation"),
  alpha = 0.05,
  boot_type = c("moving", "circular"),
  parallel = FALSE,
  ncores = parallel::detectCores() - 1,
  cl_export = NULL,
  boot_seed = NULL,
  cl = NULL,
  ...
)

Arguments

X	A vector representing observed values of the time series.
maxLag	The maximum lag to compute the estimated autocovariance function at.
x	A vector of lag indices. Defaults to the sequence 0:length(X). Must be at least as large as maxLag + 1.
n_bootstrap	The number of times to run moving block bootstrap. Defaults to 100.
l	The block length considered for bootstrap. Defaults to $\lceil N \rceil^{1/3}$, where $N$ is the length of the observation window.
estimator	The function name of the estimator to use. Defaults to standard\_est.
type	Compute either the 'autocovariance' or 'autocorrelation'. Defaults to 'autocovariance'.
alpha	The quantile used to compute the $(1 - \alpha)\%$ confidence interval. Defaults to $0.05.$
boot_type	What type of block bootstrap should be used, either 'moving' for moving block bootstrap or 'circular' for circular block bootstrap.
parallel	Whether or not the bootstrap computations should be done in parallel or not. Defaults to FALSE.
ncores	The number of cores to be used in the parallel bootstrap computations. Defaults to the number cores - 1 (threads if hyperthreading is available), calculated from parallel::detectCores() - 1.
cl_export	A vector of any additional functions or variables to export for parallel computations. This may be required if estimator is not within the package. Defaults to NULL.
boot_seed	An integer seed for reproducibility. This is used for
cl	An optional cluster object created by parallel::makeCluster. Defaults to NULL, which creates a temporary PSOCK cluster.
...	Optional named arguments to the chosen estimator. See the examples.

Details

This function performs block bootstrap to obtain a $(1-\alpha)\%$ confidence-interval for the autocovariance function. It will also compute average autocovariance function across all bootstrapped estimates.

Moving block bootstrap can be described as follows. For some times series $X(1), X(2), \dots, X(n),$ construct $k \in N$ overlapping blocks of the form $B_{i} = (X(i), \dots, X(i + \ell - 1)),$ where $\ell \in \{1, \dots , n\}$ is the block length. Randomly sample, with replacement, from the discrete uniform distribution with on $\{1, \dots, n - \ell + 1\}$ to obtain a set of random starting locations $I_{1}, \dots, I_{k}.$ Construct a bootstrapped time series $B_{1}^{\ast}, B_{2}^{\ast}, \dots, B_{k}^{\ast},$ where $B_{i}^{\ast} = B_{I_{i}}.$ The bootstrapped time series is truncated to have length $n,$ and will be of the form $X^{\ast}(1), \dots , X^{\ast}(n).$ The autocovariance function is then computed for the bootstrapped time series.

An alternative to moving block bootstrap is circular block bootstrap. Circular block bootstrap uses the time series like a circle, that is, the observation at $n + i$ is the same as the observation at location $i.$ For example, for the block $B_{n - \ell + 2}$, we obtain $(X(n - \ell + 2) , \dots , X(n), X(n + 1))$ is the same as $(X(n - \ell + 2) , \dots , X(n), X(1)).$ When performing random sampling to obtain starting locations, the set $\{1, \dots, n\}$ is now considered. The procedure for constructing the bootstrap time series after constructing the blocks and selecting the starting indices is the same as moving block bootstrap.

This process is repeated n_bootstrap times to obtain n_boostrap estimates of the autocovariance function using the bootstrapped time series, for which the average autocovariance function and the $(1 - \alpha)\%$ confidence intervals are constructed pointwise for each lag.

The choice of the block length, $\ell ,$ depends on the degree of dependence present in the time series. If the time series has a high degree of dependence, a larger block size should be chosen to ensure the dependency structure is maintained within the block.

Any estimator of the autocovariance function can be used in this function, including a custom estimator not in the package, see the examples.

When parallel = TRUE, L'Ecuyer-CMRG is used to ensure independent random numbers across workers. The serial version uses the Mersenne-Twister algorithm.

Value

A BootEsts S3 object (list) with the following values

acf_avg

A numeric vector containing the average autocovariance/autocorrelation bootstrap estimate.

lags

A numeric vector containing the lag indices used to compute the estimates on.

acf_orig

A numeric vector containing the nonbootstrapped autocovariance/autocorrelation estimate.

acf_mat

A numeric matrix of the a matrix of all of the bootstrap estimates.

conf_lower

A numeric vector containing the lower bounds for the estimated pointwise confidence interval.

conf_upper

A numeric vector containing the upper bounds for the estimated pointwise confidence interval.

est_type

The type of estimate, namely 'autocorrelation' or 'autocovariance', this depends on the argument type.

est_used

The estimator function used, this depends on the argument estimator.

boot_type

The type of estimate, namely 'moving' or 'circular', this depends on the argument boot_type.

alpha

A numeric value, which is the $\alpha$ value used to compute the confidence intervals, this depends on the argument alpha.

References

Chapters 2.5 and 2.7 in Lahiri, S. N. (2003). Resampling Methods for Dependent Data. Springer. https://doi.org/10.1007/978-1-4757-3803-2

Künsch, H. R. (1989). The Jackknife and the Bootstrap for General Stationary Observations. The Annals of Statistics 17(3), 1217-1241. https://doi.org/10.1214/aos/1176347265

Politis, D. N. & Romano, J. P. (1991). A Circular Block-Resampling Procedure for Stationary Data. In R. LePage & L. Billard, eds, Exploring the Limits of Bootstrap, Wiley, 263-270.

Examples

X <- c(1, 2, 3, 3, 2, 1)
block_bootstrap(X, 4, n_bootstrap = 3, l = 2, type = 'autocorrelation')
block_bootstrap(X, 4, n_bootstrap = 3, l = 2, type = 'autocovariance')
block_bootstrap(X, 4, n_bootstrap = 3, l = 2, estimator=tapered_est,
    rho = 0.5, window_name = 'blackman', window_params = c(0.16),
    type='autocorrelation')
my_cov_est <- function(X, maxLag) {
  n <- length(X)
  covVals <- rep(0, maxLag + 1)
  for(h in 0:maxLag) {
    covVals_t <- (X[1:(n-h)] - mean(X)) * (X[(1+h):n] - mean(X))
    covVals[h] <- sum(covVals_t) / (n - h)
  }
  return(covVals)
}
block_bootstrap(X, 4, n_bootstrap = 3, l = 2, estimator=my_cov_est)

plot(LakeHuron, main="Lake Huron Levels", ylab="Feet")
X <- as.vector(LakeHuron)
block_bootstrap(X, 20, n_bootstrap = 100, l = 40, type = 'autocorrelation')
block_bootstrap(X, 20, n_bootstrap = 100, l = 40, type = 'autocorrelation')
block_bootstrap(X, 20, n_bootstrap = 100, l = 40, estimator=tapered_est,
    rho = 0.5, window_name = 'blackman', window_params = c(0.16),
    type='autocorrelation')

my_cov_est <- function(X, maxLag) {
  n <- length(X)
  covVals <- rep(0, maxLag + 1)
  for(h in 0:maxLag) {
    covVals_t <- (X[1:(n-h)] - mean(X)) * (X[(1+h):n] - mean(X))
    covVals[h] <- sum(covVals_t) / (n - h)
  }
  return(covVals)
}
block_bootstrap(X, 20, n_bootstrap = 100, l = 40, estimator = my_cov_est,
    type = 'autocorrelation')
## Not run: 
library(parallel)
X <- c(1, 2, 3, 3, 2, 1)
my_cl <- makePSOCKcluster(2)
block_bootstrap(X, 4, n_bootstrap = 1000, l = 3, parallel = TRUE, cl = my_cl)
stopCluster(my_cl)

## End(Not run)

---

Block Bootstrap Sample

Description

This function generates block bootstrap samples for either moving block bootstrap or circular bootstrap.

Usage

bootstrap_sample(X, l, k, boot_type = c("moving", "circular"))

Arguments

X	A vector representing observed values of the time series.
l	The block length considered for bootstrap.
k	The number of blocks considered for bootstrap.
boot_type	What type of block bootstrap should be used, either 'moving' for moving block bootstrap or 'circular' for circular block bootstrap.

Details

This function generates a block bootstrap sample for a time series $X.$ For the moving block bootstrap and circular bootstrap procedures see block_bootstrap and the included references.

Value

A vector of length length(X) whose values are a bootstrapped time series.

References

Chapters 2.5 and 2.7 in Lahiri, S. N. (2003). Resampling Methods for Dependent Data. Springer. https://doi.org/10.1007/978-1-4757-3803-2

Künsch, H. R. (1989). The Jackknife and the Bootstrap for General Stationary Observations. The Annals of Statistics 17(3), 1217-1241. https://doi.org/10.1214/aos/1176347265

Politis, D. N. & Romano, J. P. (1991). A Circular Block-Resampling Procedure for Stationary Data. In R. LePage & L. Billard, eds, Exploring the Limits of Bootstrap, Wiley, 263-270.

Examples

X <- c(1, 2, 3, 3, 2, 1)
bootstrap_sample(X, 2, 3)

---

Check if an Autocovariance Function Estimate is Positive-Definite or Not.

Description

This function checks if an autocovariance function estimate is positive-definite or not by determining if the eigenvalues of the corresponding matrix (see the Details section) are all positive.

Usage

check_pd(est)

## S3 method for class 'CovEsts'
check_pd(est)

## Default S3 method:
check_pd(est)

Arguments

est	A numeric vector, corresponding cyclic matrix representing an estimated autocovariance function, or a CovEsts S3 object.

Details

For an autocovariance function estimate $\hat{C}(\cdot)$ over a set of lags separated by a constant difference $\{h_{0}, h_{1} , h_{2} , \dots , h_{n} \},$ construct the symmetric matrix

$\left[ {\begin{array}{ccccc} \hat{C}(h_{0}) & \hat{C}(h_{1}) & \cdots & \hat{C}(h_{n - 1}) & \hat{C}(h_{n}) \\ \hat{C}(h_{1}) & \hat{C}(h_{0}) & \cdots & \hat{C}(h_{n - 2}) & \hat{C}(h_{n - 1}) \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ \hat{C}(h_{n - 1}) & \hat{C}(h_{n - 2}) & \cdots & \hat{C}(h_{0}) & \hat{C}(h_{1}) \\ \hat{C}(h_{n}) & \hat{C}(h_{n - 1}) & \cdots & \hat{C}(h_{1}) & \hat{C}(h_{0}) \\ \end{array}} \right] .$

The eigendecomposition of this matrix is computed to determine if all eigenvalues are positive. If so, the estimated autocovariance function is assumed to be positive-definite.

Value

A boolean where TRUE denotes a positive-definite autocovariance function estimate and FALSE for an estimate that is not positive-definite.

Methods (by class)

• check_pd(CovEsts): Method for CovEsts objects.

• check_pd(default): Method for numeric vectors.

Examples

x <- seq(0, 5, by=0.1)
estCov <- exp(-x^2)
check_pd(estCov)

---

Kernel Correction of the Standard Estimator.

Description

This function computes the standard autocovariance estimator and applies kernel correction to it,

$\widehat{C}_{T}^{(a)}(h) = \widehat{C}(h) a_{T}(h),$

where $a_{T}(h) := a(h / N_{T}).$ It uses a kernel $a(\cdot)$ which decays or vanishes to zero (depending on the type of kernel) where $a(0) = 1.$ The rate or value at which the kernel vanishes is $N_{T}$, which is recommended to be of order $0.1 N$, where $N$ is the length of the observation window, however, one may need to play with this value.

Usage

corrected_est(
  X,
  kernel_name = c("gaussian", "exponential", "wave", "rational_quadratic", "spherical",
    "circular", "bessel_j", "matern", "cauchy"),
  kernel_params = c(),
  N_T = 0.1 * length(X),
  pd = TRUE,
  maxLag = length(X) - 1,
  x = 0:length(X),
  type = c("autocovariance", "autocorrelation"),
  meanX = mean(X),
  custom_kernel = FALSE
)

Arguments

X	A vector representing observed values of the time series.
kernel_name	The name of the kernel_ec function to be used. Possible values are: gaussian, exponential, wave, rational_quadratic, spherical, circular, bessel_j, matern, cauchy.
kernel_params	A vector of parameters of the kernel function. See kernel_ec for parameters. In the case of gaussian, wave, rational_quadratic, spherical and circular, N_T takes the place of $\theta$. For kernels that require parameters other than $\theta$, such as the Matern kernel, those parameters are passed.
N_T	The range at which the kernel function vanishes at. Recommended to be $0.1 N$ when considering all lags. This parameter may be large for a lag small estimation lag.
pd	Whether a positive-definite estimate should be used. Defaults to TRUE.
maxLag	An optional parameter that determines the maximum lag to compute the estimated autocovariance function at. Defaults to length(X) - 1.
x	A vector of lag indices. Defaults to the sequence 0:length(X). Must be at least as large as maxLag + 1.
type	Compute either the 'autocovariance' or 'autocorrelation'. Defaults to 'autocovariance'.
meanX	The average value of X. Defaults to mean(X).
custom_kernel	If a custom kernel is to be used or not. Defaults to FALSE. See examples.

Details

The aim of this estimator is gradually bring the estimated values to zero through the use of a kernel multiplier. This can be useful when estimating an autocovariance function that is short-range dependent as estimators can have large fluctuations as the lag increases, or to deal with the wave artefacts for large lags, see Bilchouris and Olenko (2025). This estimator can be positive-definite depending on whether the choice of $\widehat{C}(\cdot)$ and $a$ are chosen to be positive-definite or not.

Value

A CovEsts S3 object (list) with the following values

acf

A numeric vector containing the autocovariance/autocorrelation estimates.

lags

A numeric vector containing the lag indices used to compute the estimates on.

est_type

The type of estimate, namely 'autocorrelation' or 'autocovariance', this depends on the type parameter.

est_used

The estimator function used, in this case, 'corrected_est'.

References

Yaglom, AM (1987). Correlation Theory of Stationary and Related Random Functions. Volume I: Basic Results. Springer New York. https://doi.org/10.1007/978-1-4612-4628-2

Bilchouris, A. & Olenko, A (2025). On Nonparametric Estimation of Covariogram. Austrian Statistical Society (Vol. 54, Issue 1). https://doi.org/10.17713/ajs.v54i1.1975

Examples

X <- c(1, 2, 3)
corrected_est(X, "gaussian")

X <- rnorm(1000)
Y <- c(X[1], X[2])
for(i in 3:length(X)) { Y[i] <- X[i] - 0.3*X[i - 1] - 0.6*X[i - 2] }
plot(Y)
plot(corrected_est(Y, "bessel_j",
     kernel_params=c(0, 1), N_T=0.2*length(Y)))

# Custom kernel
my_kernel <- function(x, theta, params) {
  stopifnot(theta > 0, length(x) >= 1, all(x >= 0))
  return(sapply(x, function(t) ifelse(t == 0, 1,
         ifelse(t == Inf, 0,
         (sin((t^params[1]) / theta) / ((t^params[1]) / theta)) * cos((t^params[2]) / theta)))))
}
plot(corrected_est(Y,
     my_kernel, kernel_params=c(2, 0.25), custom_kernel = TRUE))

---

Create a Cyclic Matrix for a Given Vector.

Description

This helper function creates a symmetric matrix from a given vector $v$.

Usage

cyclic_matrix(v)

Arguments

v	A numeric vector.

Details

This function creates a symmetric matrix for a given vector $v$. If $v = \{v_{0}, v_{1} , \dots , v_{N-1} , v_{N} \},$ then the symmetric matrix will has the form

$\left[ {\begin{array}{ccccc} v_{0} & v_{1} & \cdots & v_{N - 1} & v_{N} \\ v_{1} & v_{0} & \cdots & v_{N - 2} & v_{N - 1} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ v_{N - 1} & v_{N- 2} & \cdots & v_{0} & v_{1} \\ v_{N} & v_{N - 1} & \cdots & v_{1} & v_{0} \\ \end{array}} \right]$

Value

A symmetric matrix.

Examples

## Not run: 
v <- c(1, 2, 3)
cyclic_matrix(v)

## End(Not run)

---

Compute 1D Discrete Cosine Transform

Description

This function computes the Type-II discrete cosine transform.

Usage

dct_1d(X)

Arguments

X	A vector of values for which the discrete cosine transform is being computed.

Details

The Type-II discrete cosine transform is obtained using stats::fft. If $X$ is of length N, construct a new signal $Y$ of length $4N$, with values $Y_{2n} = 0, Y_{2n + 1} = x_{n}$ for $0 \le n < N$, and $Y_{2N} = 0, Y_{4N - n} = y_{n}$ for $0 < n < 2N.$ After this, the Type-II discrete cosine transform is computed by 0.5 * Re(stats::fft(Y))[1:(length(Y) / 4)].

Value

A vector of discrete cosine transform values.

References

Ochoa-Dominguez, H. & Rao, K.R. (2019). Discrete Cosine Transform, Second Edition. CRC Press. https://doi.org/10.1201/9780203729854

Makhoul, J. (1980). A Fast Cosine Transform in One and Two Dimensions. IEEE Transactions on Acoustics, Speech, and Signal Processing 28(1), 27-34. https://doi.org/10.1109/TASSP.1980.1163351

Stasiński, R. (2002). DCT Computation Using Real-Valued DFT Algorithms. Proceedings of the 11th European Signal Processing Conference.

Examples

X <- c(1, 2, 3)
dct_1d(X)

---

Generate Spline Knots.

Description

A helper function that generates $m + 2$ spline knots of the form:

$\kappa_{0} = 0 , \kappa_{1} = 1 / (m + 1) , \dots , \kappa_{m} = m / (m + 1) , \kappa_{m + 1} = 1 .$

The knots are equally spaced with boundary knots $\kappa_{0} = 0$ and $\kappa_{m + 1} = 1 .$

Usage

generate_knots(m)

Arguments

m	The number of nonboundary knots.

Value

A numeric vector representing the knots, including the boundary knots.

References

Choi, I., Li, B. & Wang, X. (2013). Nonparametric Estimation of Spatial and Space-Time Covariance Function. JABES 18, 611-630. https://doi.org/10.1007/s13253-013-0152-z

Examples

## Not run: 
generate_knots(3)

## End(Not run)

---

Get all $\tau$.

Description

A helper function to obtain all $\tau_{i}, i=-p, \dots, m + p + 1,$ where each $\tau_{i}$ is as follows. For $i = -p , -p + 1, \dots , -2, -1 , m + 2, m + 3, \dots , m + p , m + p + 1,$ it is equal to $\tau_{i} = i / (m + 1)$, and for $i = 0, \dots , m + 1,$ it is $\tau_{i} = \kappa_{i}.$ See Choi, Li & Wang (2013, p. 615) for details.

Usage

get_taus(p, m)

Arguments

p	The order of the splines.
m	The number of nonboundary knots.

Value

A numeric vector of all $\tau_{i}, i = -p, \dots, m + p + 1.$

References

Choi, I., Li, B. & Wang, X. (2013). Nonparametric Estimation of Spatial and Space-Time Covariance Function. JABES 18, 611-630. https://doi.org/10.1007/s13253-013-0152-z

Examples

## Not run: 
get_taus(3, 2)

## End(Not run)

---

Compute Normalisation Factor

Description

This helper function is used in the computation of the normalisation factor the function tapered_est,

$H_{2, n}(0) = \sum_{s=1}^{n} a((s - 1/2) / n; \rho)^{2},$

where $a(\cdot; \cdot)$ is a window function.

Usage

H2n(
  n,
  rho,
  window_name = c("tukey", "triangular", "sine", "power_sine", "blackman",
    "hann_poisson", "welch"),
  window_params = c(1),
  custom_window = FALSE
)

Arguments

n	The sample size.
rho	A scale parameter in $(0, 1].$
window_name	The name of the window_ec function to be used. Possible values are: tukey, triangular, power_sine, blackman_window, hann_poisson, welch. Alternatively, a custom window function can be provided, see the example for taper.
window_params	A vector of parameters of the window function.
custom_window	If a custom window is to be used or not. Defaults to FALSE.

Value

A single value being $H_{2, n}(0)$.

Examples

## Not run: 
H2n(3, 0.6, "tukey")

## End(Not run)

---

Hilbert-Schmidt Norm Between Estimated Autocovariance Functions.

Description

This function computes the Hilbert-Schidmt norm between two estimated autocovariance functions.

Usage

hilbert_schmidt(est1, est2)

## S3 method for class 'CovEsts'
hilbert_schmidt(est1, est2)

## Default S3 method:
hilbert_schmidt(est1, est2)

Arguments

est1	A numeric vector representing the first estimated autocovariance function.
est2	A numeric vector of the same length as est1 representing the second estimated autocovariance function

Details

This function computes the Hilbert-Schidmt norm between two estimated autocovariance functions. The Hilbert-Schmidt norm of a matrix

$D = \left[(d_{i,j})_{1 \le i,j \le n}\right] = \left[ {\begin{array}{ccccc} D(h_{0}) & D(h_{1}) & \cdots & D(h_{n - 1}) & D(h_{n}) \\ D(h_{1}) & D(h_{0}) & \cdots & D(h_{n - 2}) & D(h_{n - 1}) \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ D(h_{n - 1}) & D(h_{n - 2}) & \cdots & D(h_{0}) & D(h_{1}) \\ D(h_{n}) & D(h_{n - 1}) & \cdots & D(h_{1}) & D(h_{0}) \\ \end{array}} \right] ,$

over a set of lags $\{h_{0}, h_{1}, \dots , h_{N} \},$ where $D(h) = \hat{C}_{1}(h) - \hat{C}_{2}(h),$ is defined as

${\left\Vert D \right\Vert}_{HS} = \sqrt{\sum_{i,j} d_{i, j}^{2}}.$

Value

A numeric value representing the estimated Hilbert-Schmidt norm between two estimated autocovariance functions.

Methods (by class)

• hilbert_schmidt(CovEsts): Method for CovEsts objects.

• hilbert_schmidt(default): Method for numeric vectors.

Examples

x <- seq(0, 5, by=0.1)
estCov1 <- exp(-x^2)
estCov2 <- exp(-x^2.1)
hilbert_schmidt(estCov1, estCov2)

---

Compute 1D Inverse Discrete Cosine Transform

Description

This function computes the inverse of the Type-II discrete cosine transform.

Usage

idct_1d(X)

Arguments

X	A vector of values for which the discrete cosine transform is being computed.

Details

The Type-II inverse discrete cosine transform is computed using stats::fft. For autocovariance function estimation, the spectrum is given in the input X and then an inverse FFT is applied.

The original spectrum, dct_full, from X is obtained as follows. First, the original sample Fourier spectrum is reconstructed using dct_full <- c(X, 0, -X[-1], 0, rev(X[-1])). After this, an inverse FFT is applied, idct <- Re(stats::fft(dct_full, inverse = TRUE)) * (2 / length(dct_full)), which gives the original function with additional zero-values at even indices. The zeroes are dropped, which gives the untransformed X.

Value

A vector with the inverse transform values.

References

Ochoa-Dominguez, H. & Rao, K.R. (2019). Discrete Cosine Transform, Second Edition. CRC Press. https://doi.org/10.1201/9780203729854

endolith (2013). Fast Cosine Transform via FFT. Signal Processing Stack Exchange. https://dsp.stackexchange.com/a/10606

Examples

X <- dct_1d(c(1, 2, 3))
idct_1d(X)

---

1D Isotropic Kernels.

Description

This function computes one of the isotropic kernels listed below. Note that kernel() is deprecated, please use kernel_ec() instead. Unlike kernel_symm_ec, these kernels are only defined when $x \geq 0.$ They are used as kernel multipliers in estimators corrected_est and kernel_est.

Usage

kernel_ec(
  x,
  name = c("gaussian", "exponential", "wave", "rational_quadratic", "spherical",
    "circular", "bessel_j", "matern", "cauchy"),
  params = c(1)
)

Arguments

x	A vector or matrix of arguments of at least length 1 for which the kernel is computed at.
name	The name of the kernel. Options are: gaussian, exponential, wave, rational_quadratic, spherical, circular, bessel_j, matern, and cauchy.
params	A vector of parameters for the kernel. See the documentation below for the position of the parameters. All kernels have a scale parameter as the first value in the vector.

Details

Gaussian Kernel. The isotropic Gaussian kernel, which is positive-definite for ${R}^{d}, d \in N,$ is defined as

$a(x;\theta) = \exp(-x^{2} / \theta), \theta > 0 .$

The params argument is of the form c($\theta$).

Exponential Kernel. The isotropic exponential kernel, which is positive-definite for ${R}^{d}, d \in N,$ is defined as

$a(x;\theta) = \exp(-x / \theta) , \theta > 0 .$

The params argument is of the form c($\theta$).

Isotropic Wave (Cardinal Sine) Kernel. The isotropic wave (cardinal sine) kernel, which is positive-definite for ${R}^{d}, d \leq 3,$ is given by

$a(x;\theta) = \left\{ \begin{array}{ll} \frac{\theta}{x} \sin\left( \frac{x}{\theta} \right), & x \neq 0 \\ 1, & x = 0 \end{array} \right.$

where $\theta > 0 .$ The params argument is of the form c($\theta$).

Isotropic Rational Quadratic Kernel. The isotropic rational quadratic kernel, which is positive-definite for ${R}^{d}, d \in N,$ is defined as

$a(x;\theta) = \frac{\theta}{x^{2} + \theta} , \theta > 0 .$

The params argument is of the form c($\theta$).

Isotropic Spherical Kernel. The isotropic spherical kernel, which is positive-definite for ${R}^{3}, d \leq 3,$ is given by

$a(x;\theta) = \left\{ \begin{array}{ll} 1 - \frac{3}{2}\frac{x}{\theta} + \frac{1}{2}\left( \frac{x}{\theta} \right)^{3}, & x < \theta \\ 0, & \mbox{otherwise} \end{array} \right.$

where $\theta > 0 .$ The params argument is of the form c($\theta$).

Isotropic Circular Kernel. The isotropic circular kernel, which is positive-definite for ${R}^{d}, d \leq 2,$ is given by

$a(x;\theta) = \left\{ \begin{array}{ll} \frac{2}{\pi}\arccos\left( \frac{x}{\theta} \right) - \frac{2}{\pi}\frac{x}{\theta} \sqrt{ 1 - \left( \frac{x}{\theta} \right)^{2} }, & x < \theta \\ 0, & \mbox{otherwise} \end{array} \right.$

where $\theta > 0 .$ The params argument is of the form c($\theta$).

Isotropic Matérn Kernel. The isotropic Matérn kernel, which is positive-definite for ${R}^{d}, d \in N,$ and when $\nu > 0,$ is defined as

$a(x; \theta, \nu) = \left(\sqrt{2\nu} \frac{x}{\theta} \right)^{\nu} \left(2^{\nu - 1} \Gamma(\nu) \right)^{-1} K_{\nu}\left( \sqrt{2\nu} \frac{x}{\theta} \right) ,$

where $\theta > 0$ and $K_{\nu}(\cdot)$ is the modified Bessel function of the second kind. The params argument is of the form c($\theta, \nu$).

Isotropic Bessel Kernel. The isotropic Bessel kernel, which is positive-definite for ${R}^{d}, d \in N,$ and when $\nu \geq \frac{d}{2} - 1,$ is given by

$a(x; \theta, \nu) = 2^{\nu} \Gamma(\nu + 1) J_{\nu}(x / \theta) (x / \theta)^{-\nu} ,$

where $\theta > 0$ and $J_{\nu}(\cdot)$ is the Bessel function of the first kind. The params argument is of the form c($\theta, \nu, d$ ).

Isotropic Cauchy Kernel. The isotropic Cauchy kernel, which is positive-definite for ${R}^{d}, d \in N,$ and when $0 < \alpha \leq 2$ and $\beta \geq 0,$ is defined by

$a(x ; \theta, \alpha, \beta) = (1 + (x / \theta)^{\alpha})^{-(\beta / \alpha)} , \theta > 0 .$

The params argument is of the form c($\theta, \alpha, \beta$ ).

Value

A vector or matrix of kernel values.

References

Genton, M. (2001). Classes of Kernels for Machine Learning: A Statistics Perspective. Journal of Machine Learning Research. 2, 299-312. https://doi.org/10.1162/15324430260185646

Table 4.2 in Hristopulos, D. T. (2020). Random Fields for Spatial Data Modeling: A Primer for Scientists and Engineers. Springer. https://doi.org/10.1007/978-94-024-1918-4

Examples

x <- c(0.2, 0.4, 0.6)
theta <- 0.9
kernel_ec(x, "gaussian", c(theta))
kernel_ec(x, "exponential", c(theta))
kernel_ec(x, "wave", c(theta))
kernel_ec(x, "rational_quadratic", c(theta))
kernel_ec(x, "spherical", c(theta))
kernel_ec(x, "circular", c(theta))
nu <- 1
kernel_ec(x, "matern", c(theta, nu))
dim <- 1
kernel_ec(x, "bessel_j", c(theta, nu, dim))
alpha <- 1
beta <- 2
kernel_ec(x, "cauchy", c(theta, alpha, beta))
curve(kernel_ec(x, "gaussian", c(theta)), from = 0, to = 5)
curve(kernel_ec(x, "exponential", c(theta)), from = 0, to = 5)
curve(kernel_ec(x, "wave", c(theta)), from = 0, to = 5)
curve(kernel_ec(x, "rational_quadratic", c(theta)), from = 0, to = 5)
curve(kernel_ec(x, "spherical", c(theta)), from = 0, to = 5)
curve(kernel_ec(x, "circular", c(theta)), from = 0, to = 5)
curve(kernel_ec(x, "matern", c(theta, nu)), from = 0, to = 5)
curve(kernel_ec(x, "bessel_j", c(theta, nu, dim)), from = 0, to = 5)
curve(kernel_ec(x, "cauchy", c(theta, alpha, beta)), from = 0, to = 5)

---

Kernel Correction for an Estimated Autocovariance Function.

Description

This function applies kernel correction to an estimated autocovariance function,

$\widehat{C}_{T}^{(a)}(h) = \widehat{C}(h) a_{T}(h),$

where $a_{T}(h) := a(h / N_{T}).$ It uses a kernel $a(\cdot)$ which decays or vanishes to zero (depending on the type of kernel) where $a(0) = 1.$ The rate or value at which the kernel vanishes is $N_{T}$, which is recommended to be of order $0.1 N$, where $N$ is the length of the observation window, however, one may need to play with this value.

Usage

kernel_est(
  estCov,
  kernel_name = c("gaussian", "exponential", "wave", "rational_quadratic", "spherical",
    "circular", "bessel_j", "matern", "cauchy"),
  kernel_params = c(),
  N_T = 0.1 * length(estCov),
  maxLag = length(estCov) - 1,
  x = 0:length(estCov),
  type = c("autocovariance", "autocorrelation"),
  custom_kernel = FALSE
)

## S3 method for class 'CovEsts'
kernel_est(
  estCov,
  kernel_name = c("gaussian", "exponential", "wave", "rational_quadratic", "spherical",
    "circular", "bessel_j", "matern", "cauchy"),
  kernel_params = c(),
  N_T = 0.1 * length(estCov$acf),
  maxLag = length(estCov$acf) - 1,
  x = estCov$lags,
  type = c("autocovariance", "autocorrelation"),
  custom_kernel = FALSE
)

## Default S3 method:
kernel_est(
  estCov,
  kernel_name = c("gaussian", "exponential", "wave", "rational_quadratic", "spherical",
    "circular", "bessel_j", "matern", "cauchy"),
  kernel_params = c(),
  N_T = 0.1 * length(estCov),
  maxLag = length(estCov) - 1,
  x = 0:length(estCov),
  type = c("autocovariance", "autocorrelation"),
  custom_kernel = FALSE
)

Arguments

estCov	A vector whose values are an estimate autocovariance function.
kernel_name	The name of the kernel_ec function to be used. Possible values are: gaussian, exponential, wave, rational_quadratic, spherical, circular, bessel_j, matern, cauchy.
kernel_params	A vector of parameters of the kernel function. See kernel_ec for parameters. In the case of gaussian, wave, rational_quadratic, spherical and circular, N_T takes the place of $\theta$. For kernels that require parameters other than $\theta$, such as the Matern kernel, those parameters are passed.
N_T	The range at which the kernel function vanishes at. Recommended to be $0.1 N$ when considering all lags. This parameter may be large for a lag small estimation lag.
maxLag	An optional parameter that determines the maximum lag to compute the estimated autocovariance function at. Defaults to length(estCov) - 1.
x	A vector of lag indices. Defaults to the sequence 0:length(X). Must be at least as large as maxLag + 1.
type	Compute either the 'autocovariance' or 'autocorrelation'. Defaults to 'autocovariance'.
custom_kernel	If a custom kernel is to be used or not. Defaults to FALSE. See the examples of corrected_est for usage.

Value

A vector whose values are the kernel corrected autocovariance estimates or CovEsts S3 object (list) with the following values

acf

A numeric vector containing the autocovariance/autocorrelation estimates.

lags

A numeric vector containing the lag indices used to compute the estimates on, inherited from the argument estCov.

est_type

The type of estimate, namely 'autocorrelation' or 'autocovariance', this depends on the argument type.

est_used

The estimator function used, in this case, 'kernel_est'.

If a numeric vector is given for the argument estCov, then a numeric vector output is given, and if a CovEsts S3 object is given, a CovEsts object is given as output.

Methods (by class)

• kernel_est(CovEsts): Method for CovEsts objects.

• kernel_est(default): Method for numeric vectors.

Examples

X <- rnorm(1000)
Y <- c(X[1], X[2])
for(i in 3:length(X)) { Y[i] <- X[i] - 0.3*X[i - 1] - 0.6*X[i - 2] }
cov_est <- standard_est(Y)
plot(cov_est)
plot(kernel_est(cov_est,
     "bessel_j", kernel_params=c(0, 1), N_T=0.2*length(Y)))

---

1D Isotropic Symmetric Kernels.

Description

These functions computes values of kernels that have the properties of symmetric probability distributions. Note that kernel_symm() is deprecated, please use kernel_symm_ec() instead. For a kernel $a(x)$, the standardised version is $a(x) / \int_{-\infty}^{\infty} a(x) dx$, so that the integral is 1. The symmetric kernels are different to kernel and are used in the functions adjusted_est and truncated_est.

Usage

kernel_symm_ec(
  x,
  name = c("gaussian", "wave", "rational_quadratic", "bessel_j"),
  params = c(1)
)

Arguments

x	A vector or matrix of arguments of at least length 1 for which the kernel is computed at. Each value can be negative as well as positive.
name	The name of the kernel. Options are: gaussian, wave, rational_quadratic, bessel_j.
params	A vector of parameters for the kernel. See the documentation below for the position of the parameters. All kernels will have a scale parameter as the first value in the vector.

Details

Symmetric Gaussian Kernel. The symmetric Gaussian kernel is defined as

$a(x;\theta) = \sqrt{\pi \theta} \exp(-x^{2} / \theta) , \theta > 0.$

The params argument is of the form c($\theta$).

Symmetric Wave Kernel. The wave (cardinal sine) kernel is given by

$a(x;\theta) = \left\{ \begin{array}{ll} (\sqrt{\theta^{2}} \pi)^{-1} \frac{\theta}{x} \sin\left( \frac{x}{\theta} \right), & x \neq 0 \\ 1, & x = 0 \end{array} \right.$

where $\theta > 0.$ The params argument is of the form c($\theta$)

Symmetric Rational Quadratic Kernel. The symmetric rational quadratic kernel is given by

$a(x;\theta) = (\pi \sqrt{\theta})^{-1} \left(\frac{\theta}{x^{2} + \theta}\right), \theta > 0.$

The params argument is of the form c($\theta$)

Symmetric Besesel Kernel. The symmetric Bessel kernel, which is valid when $\nu \geq \frac{d}{2} - 1,$ is given by

$a(x; \theta, \nu) = \left(\Gamma\left(\frac{1}{2} + \nu\right)/(2 \sqrt{\pi} \theta \Gamma(1 + \nu))\right) ( 2^{\nu} \Gamma(\nu + 1) J_{\nu}(x / \theta) (x / \theta)^{-\nu}), \,\theta > 0, \nu \geq \frac{d}{2} - 1,$

where $J_{\nu}(\cdot)$ is the Bessel function of the first kind and $d$ is the dimension. The params argument is of the form c($\theta, \nu, d$).

Value

A vector or matrix of values.

Examples

x <- c(-2, -1, 0, 1, 2)
theta <- 1
kernel_symm_ec(x, "gaussian", c(theta))
kernel_symm_ec(x, "wave", c(theta))
kernel_symm_ec(x, "rational_quadratic", c(theta))
dim <- 1
nu <- 1
kernel_symm_ec(x, "bessel_j", c(theta, nu, dim))
curve(kernel_symm_ec(x, "gaussian", c(theta)), from = -5, to = 5)
curve(kernel_symm_ec(x, "wave", c(theta)), from = -5, to = 5)
curve(kernel_symm_ec(x, "rational_quadratic", c(theta)), from = -5, to = 5)
curve(kernel_symm_ec(x, "bessel_j", c(theta, nu, dim)), from = -5, to = 5)

---

Lines Method for BootEsts Objects

Description

This function plots a BootEsts object as a line onto another plot.

Usage

## S3 method for class 'BootEsts'
lines(x, type = "l", ...)

Arguments

x	A BootEsts S3 object.
type	Defaults to 'l', see base::plot for the possible values.
...	Additional plotting arguments, refer to graphics::par.

Value

A line of a BootEsts S3 object.

Examples

plot(block_bootstrap(c(1, 2, 3), 2))
lines(block_bootstrap(c(1, 2, 3), 2, pd = FALSE), lwd = 2, lty = 3, col = 'blue')

---

Lines Method for CovEsts Objects

Description

This function plots a CovEsts object as a line onto another plot.

Usage

## S3 method for class 'CovEsts'
lines(x, type = "l", ...)

Arguments

x	A CovEst S3 object.
type	Defaults to 'l', see base::plot for the possible values.
...	Additional plotting arguments, refer to graphics::par.

Value

A line of a CovEsts S3 object.

Examples

plot(standard_est(c(1, 2, 3)))
lines(standard_est(c(1, 2, 3), pd = FALSE))

---

Lines Method for VarioEsts Objects

Description

This function plots a VarioEsts object as a line onto another plot.

Usage

## S3 method for class 'VarioEsts'
lines(x, type = "l", ...)

Arguments

x	A VarioEsts S3 object.
type	Defaults to 'l', see base::plot for the possible values.
...	Additional plotting arguments, refer to graphics::par.

Value

A line of a VarioEsts S3 object.

Examples

plot(to_vario(standard_est(c(1, 2, 3))))
lines(to_vario(standard_est(c(1, 2, 3), pd = FALSE)))

---

Make a Function Positive-Definite

Description

This function can make a function positive-definite using the methods proposed by P. Hall and his coauthors.

Usage

make_pd(x, method.1 = TRUE)

## S3 method for class 'CovEsts'
make_pd(x, method.1 = TRUE)

## Default S3 method:
make_pd(x, method.1 = TRUE)

Arguments

x	A vector of numeric values of an estimated autocovariance function.
method.1	Should method 1 be used (TRUE) or method 2 (FALSE).

Details

This function perform positive-definite adjustments proposed by P. Hall and his coauthors.

Method 1 is as follows:

1. Take the discrete Fourier cosine transform, $\widehat{\mathcal{F}}(\theta)$.

2. Compute a modified spectrum $\widetilde{\mathcal{F}}(\theta) = \max(\widehat{\mathcal{F}}(\theta), 0)$ for all sample frequencies.

3. Perform the Fourier inversion to obtain a new estimator.

Method 2 is as follows:

1. Take the discrete Fourier cosine transform $\widehat{\mathcal{F}}(\theta)$.

2. Find the smallest frequency where its associated value in the spectral domain is negative

   $\hat{\theta} = \inf\{ \theta > 0 : \widehat{\mathcal{F}}(\theta) < 0\}.$

3. Compute a modified spectrum $\widetilde{\mathcal{F}} = \widehat{\mathcal{F}}(\theta)\textbf{1}(\theta < \hat{\theta}),$ where $\textbf{1}(A)$ is the indicator function over the set $A.$

4. Perform the Fourier inversion.

Value

A vector that is the adjusted function or a CovEsts S3 object (list) with the following values

acf

A numeric vector containing the autocovariance/autocorrelation estimates.

lags

A numeric vector containing the lag indices used to compute the estimates on, inherited from the argument x.

est_type

The type of estimate, namely 'autocorrelation' or 'autocovariance', inherited from the argument x.

est_used

The estimator function used, inherited from the argument x.

correction_method

Either 'method.1' or 'method.2' depending on the argument method.1

If a numeric vector is given for the argument x, then a numeric vector output is given, and if a CovEsts S3 object is given, a CovEsts object is given as output.

Methods (by class)

• make_pd(CovEsts): Method for CovEsts objects.

• make_pd(default): Method for numeric vectors.

References

Hall, P. & Patil, P. (1994). Properties of nonparametric estimators of autocovariance for stationary random fields. Probability Theory and Related Fields 99(3), 399-424. https://doi.org/10.1007/bf01199899

Hall, P., Fisher, N. I., & Hoffmann, B. (1994). On the nonparametric estimation of covariance functions. The Annals of Statistics 22(4), 2115-2134. https://doi.org/10.1214/aos/1176325774

Bilchouris, A. & Olenko, A (2025). On Nonparametric Estimation of Covariogram. Austrian Statistical Society 54(1), 112-137. https://doi.org/10.17713/ajs.v54i1.1975

Examples

X <- c(1, 2, 3, 4)
make_pd(X)
check_pd(make_pd(X))
check_pd(make_pd(X, method.1 = FALSE))

---

Maximum Vertical Distance Between Estimated Functions.

Description

This function computes the maximum vertical distance between functions.

Usage

max_distance(est1, est2, lags = c(), plot = FALSE)

## S3 method for class 'CovEsts'
max_distance(est1, est2, lags = c(), plot = FALSE)

## Default S3 method:
max_distance(est1, est2, lags = c(), plot = FALSE)

Arguments

est1	A numeric vector representing the first estimated autocovariance function.
est2	A numeric vector of the same length as est1 representing the second estimated autocovariance function
lags	An optional vector of lags starting from 0 up until some other lag. If empty, a vector of lags is created starting from 0 until len(est1) - 1, by 1.
plot	A boolean as to whether a plot should be created. By default, no plot is created.

Details

This function computes the maximum vertical distance between functions:

$D(\hat{C}_{1}(h), \hat{C}_{2}(h)) = \displaystyle \max_{h} \left| \hat{C}_{1}(h) - \hat{C}_{2}(h) \right| ,$

where $\hat{C}_{1}(\cdot)$ and $\hat{C}_{2}(\cdot)$ are estimated autocovariance functions. It assumes that the estimated values are given for the same set of lags. The vectors of the function values must be of the same length.

Value

A numeric value representing the maximum vertical distance between the two estimated functions.

Methods (by class)

• max_distance(CovEsts): Method for CovEsts objects.

• max_distance(default): Method for numeric vectors.

Examples

x <- seq(0, 5, by=0.1)
estCov1 <- exp(-x^2)
estCov2 <- exp(-x^2.1)
max_distance(estCov1, estCov2, lags=x)
max_distance(estCov1, estCov2, lags=x, plot = TRUE)

---

MSE Between Estimated Autocovariance Functions.

Description

This function computes the mean-square difference/error between two autocovariance functions (estimated or theoretical).

Usage

mse(est1, est2)

## S3 method for class 'CovEsts'
mse(est1, est2)

## Default S3 method:
mse(est1, est2)

Arguments

est1	A numeric vector representing the first estimated autocovariance function.
est2	A numeric vector of the same length as est1 representing the second estimated (or theoretical) autocovariance function

Details

This function computes the mean-square difference/error (MSE) between two estimated autocovariance functions (estimated or theoretical). The MSE is defined as

$\frac{1}{n} \sum_{i=0}^{n} \left(\widehat{C}_{1}(h_{i}) - \widehat{C}_{2}(h_{i})\right)^{2}$

over a set of lags $\{h_{0}, h_{1} , h_{2} , \dots , h_{n} \}.$

Value

A numeric value representing the MSE between two autocovariance functions (estimated or theoretical).

Methods (by class)

• mse(CovEsts): Method for CovEsts object.

• mse(default): Method for numeric vectors.

Examples

x <- seq(0, 5, by=0.1)
estCov1 <- exp(-x^2)
estCov2 <- exp(-x^2.1)
mse(estCov1, estCov2)

---

Compute the Nearest Positive-Definite Matrix.

Description

This function computes the nearest positive-definite matrix to some matrix $A$.

Usage

nearest_pd(X)

## S3 method for class 'CovEsts'
nearest_pd(X)

## Default S3 method:
nearest_pd(X)

Arguments

X	Either a numeric square matrix. If a vector is provided, a matrix will be created of the form found in cyclic_matrix.

Details

This function computes the nearest positive-definite matrix to some matrix $A$.

The procedure to do so is as follows For a matrix $X$, compute the symmetric matrix $B = (A + A^{T}) / 2.$ Let $B = UH$ be the polar decomposition of B. The nearest positive-definite matrix to $X$ is $X_{F} = (B + H) / 2.$

Unlike shrinking, only an autocorrelation matrix can be returned, not an autocovariance function.

The implementation is a translation of https://au.mathworks.com/matlabcentral/fileexchange/42885-nearestspd#functions_tab .

Value

The closest positive-definite autocorrelation matrix.

Methods (by class)

• nearest_pd(CovEsts): Method for CovEsts objects.

• nearest_pd(default): Method for numeric vectors.

References

Higham, N. J. (1988). Computing a nearest symmetric positive semidefinite matrix. Linear Algebra and its Applications, 103, 103-118. https://doi.org/10.1016/0024-3795(88)90223-6

D'Errico, J. (2025). nearestSPD (https://www.mathworks.com/matlabcentral/fileexchange/42885-nearestspd), MATLAB Central File Exchange. Retrieved August 2, 2025.

Examples

X <- c(1, 0, -1.1)
nearest_pd(X)
check_pd(nearest_pd(X))

---

Normalise a CovEsts Object

Description

This function normalises a CovEsts S3 object, that is, turning an autocovariance function into an autocorrelation function. This procedure is a one-way transformation. This function does not support any other type of argument.

Usage

normalise_acf(estCov)

## Default S3 method:
normalise_acf(estCov)

## S3 method for class 'CovEsts'
normalise_acf(estCov)

Arguments

estCov

A CovEst S3 object.

Value

A CovEsts S3 object with a normalised acf attribute.

Methods (by class)

• normalise_acf(default): Method for other objects.

• normalise_acf(CovEsts): Method for 'CovEsts' objects.

Examples

normalise_acf(standard_est(c(1, 2, 3)))

---

Plot Method for BootEsts Objects

Description

This function plots a BootEsts object. It plots the confidence region, average bootstrap estimate, and the original nonbootstrapped estimate.

Usage

## S3 method for class 'BootEsts'
plot(x, type = "l", xlab = "Lag (h)", ylab = "ACF", ...)

Arguments

x	A BootEsts S3 object.
type	Defaults to 'l', see base::plot for the possible values.
xlab	Defaults to 'Lag (h)', see graphics::title for more information.
ylab	Defaults to 'ACF', see graphics::title for more information.
...	Additional plotting arguments, refer to graphics::par.

Value

A plot of a BootEsts S3 object.

Examples

plot(block_bootstrap(c(1, 2, 3), 2))

---

Plot Method for CovEsts Objects

Description

This function plots a CovEsts object, over the lags used to estimate the autocovariance values.

Usage

## S3 method for class 'CovEsts'
plot(x, type = "l", xlab = "Lag (h)", ylab = "ACF", ...)

Arguments

x	A CovEst S3 object.
type	Defaults to 'l', see base::plot for the possible values.
xlab	Defaults to 'Lag (h)', see graphics::title for more information.
ylab	Defaults to 'ACF', see graphics::title for more information.
...	Additional plotting arguments, refer to graphics::par.

Value

A plot of a CovEsts S3 object.

Examples

plot(standard_est(c(1, 2, 3)))

---

Plot Method for VarioEsts Objects

Description

This function plots a VarioEsts object, over the lags used to estimate the variogram values.

Usage

## S3 method for class 'VarioEsts'
plot(x, type = "l", xlab = "Lag (h)", ylab = expression(gamma(h)), ...)

Arguments

x	A VarioEsts S3 object.
type	Defaults to 'l', see base::plot for the possible values.
xlab	Defaults to 'Lag (h)', see graphics::title for more information.
ylab	Defaults to expression(gamma(h)), see graphics::title for more information.
...	Additional plotting arguments, refer to graphics::par.

Value

A plot of a VarioEsts S3 object.

Examples

plot(to_vario(standard_est(c(1, 2, 3))))

---

Print Method for BootEsts Objects

Description

This function plots a BootEsts object, printing estimtor type, bootstrap type, estimator used, the bootstrap average estimate, the original estimate, and the pointwise confidence interval.

Usage

## S3 method for class 'BootEsts'
print(x, n_head = 5, n_tail = 5, digits = NULL, ...)

Arguments

x	A BootEsts S3 object.
n_head	The number of 'head' values to print, defaults to 5.
n_tail	The number of 'tail' values to print, defaults to 5.
digits	The number of digits to print. Defaults to NULL. Note, this value is passed into base::format.
...	Addition printing parameters, see base::print.

Value

Prints a BootEsts object.

Examples

print(block_bootstrap(c(1, 2, 3), 2))

---

Print Method for CovEsts Objects

Description

This function plots a CovEsts object, printing the estimator type, estimator used, lags and values.

Usage

## S3 method for class 'CovEsts'
print(x, n_head = 5, n_tail = 5, digits = NULL, ...)

Arguments

x	A CovEst S3 object.
n_head	The number of 'head' values to print, defaults to 5.
n_tail	The number of 'tail' values to print, defaults to 5.
digits	The number of digits to print. Defaults to NULL. Note, this value is passed into base::format.
...	Addition printing parameters, see base::print.

Value

Prints a CovEsts object.

Examples

print(standard_est(c(1, 2, 3)))

---

Print Method for VarioEsts Objects

Description

This function plots a VarioEsts object, printing the lags and values.

Usage

## S3 method for class 'VarioEsts'
print(x, n_head = 5, n_tail = 5, digits = NULL, ...)

Arguments

x	A VarioEsts S3 object.
n_head	The number of 'head' values to print, defaults to 5.
n_tail	The number of 'tail' values to print, defaults to 5.
digits	The number of digits to print. Defaults to NULL. Note, this value is passed into base::format.
...	Addition printing parameters, see base::print.

Value

Prints a VarioEsts object.

Examples

print(to_vario(standard_est(c(1, 2, 3))))

---