nearPDToeplitz {dbacf}

Computes the nearest positive definite Toeplitz matrix

Description

Computes the nearest positive definite Toeplitz matrix to an initial approximation, typically a covariance (correlation) matrix of a stationary process. This function implements an alternating projection algorithm that combines Grigoriadis et al. (1994) and Higham (2002). For further details see Section 5 of Tecuapetla-Gómez and Munk (2017).

Usage

nearPDToeplitz(
  matrix,
  type = c("covariance", "correlation"),
  toleranceEigen = 1e-06,
  toleranceConvergence = 1e-06,
  tolerancePosDef = 1e-06,
  maxIterations = 100,
  doEigen = TRUE
)

Arguments

matrix

a symmetric matrix.

type

string indicating whether the elements of the main diagonal must be all equal to 1 ("correlation") or not ("covariance").

toleranceEigen

defines relative positiveness of eigenvalues compared to largest one.

toleranceConvergence

numeric indicating convergence tolerance for alternating projection algorithm.

tolerancePosDef

tolerance for forcing positive definiteness (in the final step) of alternating projection algorithm.

maxIterations

integer giving maximum number of iterations allowed in alternating projection algorithm; when this number is exceeded without convergence a warning is displayed and matrix computed in step maxIterations of algorithm is returned.

doEigen

logical indicating whether finding the closest positive definite matrix -through a eigen step- should be applied to the result of the alternating projection algorithm.

Details

This function is based on an alternating projection algorithm which involves projecting sequentially and iteratively the initial matrix into the set of symmetric positive definite and into the space of Toeplitz matrices, respectively. The iteration process will stop because either a criterion of convergence is met or maxIterations has been exceeded (without convergence). Criterion of convergence: if the Frobenius norm of the difference of the projection matrices computed in the last two iterations of the algorithm is smaller than toleranceConvergence, then the algorithm stops returning the projection matrix computed in the last iteration.

When projecting onto the set of symmetric positive definite matrices, toleranceEigen controls the relative magnitude of any eigenvalue \lambda_k with respect to the largest one \lambda_1 and all eigenvalues \lambda_k are treated as zero if \lambda_k / \lambda_1 \leq \code{toleranceEigen}.

Value

A list containing:

projection

a matrix, the computed symmetric positive definite Toeplitz matrix.

normF

Frobenius norm of the difference between original and projection matrix.

iterations

number of iterations used for alternating projection algorithm.

relativeTolerance

numeric giving relative error (in Frobenius norm) of final approximation with respect to original matrix.

converged

logical indicating if alternating projection algorithm converged.

References

Grigoriadis, K.M., Frazho, A., Skelton, R. (1994). Application of alternating convex projection methods for computation of positive Toeplitz matrices, IEEE Transactions on signal processing 42(7), 1873–1875.

N. Higham (2002). Computing the nearest correlation matrix - a problem from finance, Journal of Numerical Analysis 22, 329–343.

Tecuapetla-Gómez, I and Munk, A. (2017). Autocovariance estimation in regression with a discontinuous signal and m-dependent errors: A difference-based approach. Scandinavian Journal of Statistics, 44(2), 346–368.

See Also

nearPD, projectToeplitz, symBandedToeplitz, posdefify

Examples

# Higham (2002), p. 334
(mat <- matrix(c(1, 1, 0, 1, 1, 1, 0, 1, 1), byrow = TRUE, ncol = 3))
matProj <- matrix(c(1, 0.7607, 0.1573, 0.7607, 1, 0.7607, 0.1573, 0.7607, 1), 
                  byrow = TRUE, ncol = 3)
nrPDT.mat <- nearPDToeplitz(mat, type = "correlation")
stopifnot( identical(unname(matProj), unname(round(as.matrix(nrPDT.mat$projection), 
                     digits=4) ) ))
eigen(nrPDT.mat$projection)$values

# Toeplitz banded matrix near to the covariance matrix of 100 realizations
# of an MA(5) with following parameters:

n <- 1e2
alphas <- c(-2, 0.5, -4, 0, 0.75)
(true.acf <- ARMAacf(ma = alphas))
alphasMat <- symBandedToeplitz(true.acf, n = n)
stopifnot( min(eigen(alphasMat)$values) > 0 ) # alphasMat is a positive definite matrix

(l <- length(true.acf))
(acf.modified <- c(true.acf[-c(l - 1, l)], 0.25)) # modifying original acf
x <- acf.modified
acfMat <- symBandedToeplitz(x, n = n)

# no. of non positive eigenvalues of acfMat (6)
length( eigen(acfMat)$values[eigen(acfMat)$values < 0 ] )
# acfMat is a 100 x 100 symmetric banded Toeplitz matrix
acfMat[1:15, 1:30]

system.time(nrPDT.acfMat <- nearPDToeplitz(acfMat, type = "correlation"))
y <- eigen(nrPDT.acfMat$projection)$values
# no. of non positive eigenvalues of nrPDT.acfMat
length( y[ y < 0 ] ) # none!