Nothing
#' Maximum Likelihood Estimation of an R-Vine Copula Model
#'
#' This function calculates the maximum likelihood estimate (MLE) of the
#' R-vine copula model parameters using sequential estimates as initial values
#' (if not provided).
#'
#'
#' @param data An N x d data matrix (with uniform margins).
#' @param RVM An [RVineMatrix()] object including the structure and
#' the pair-copula families and parameters (if known).
#' @param start Lower triangular d x d matrix with zero diagonal entries with
#' starting values for the pair-copula parameters (optional; otherwise they are
#' calculated via \cr [RVineSeqEst()]; default: `start =
#' RVM$par`).
#' @param start2 Lower triangular d x d matrix with zero diagonal entries with
#' starting values for the second parameters of pair-copula families with two
#' parameters (optional; otherwise they are calculated via
#' [RVineSeqEst()]; default: `start2 = RVM$par2`).
#' @param maxit The maximum number of iteration steps (optional; default:
#' `maxit = 200`).
#' @param max.df Numeric; upper bound for the estimation of the degrees of
#' freedom parameter of the t-copula (default: `max.df = 30`; for more
#' details see [BiCopEst()]).
#' @param max.BB List; upper bounds for the estimation of the two parameters
#' (in absolute values) of the BB1, BB6, BB7 and BB8 copulas \cr (default:
#' `max.BB = list(BB1=c(5,6),BB6=c(6,6),BB7=c(5,6),BB8=c(6,1))`).
#' @param grad If RVM$family only contains one parameter copula families or the
#' t-copula the analytical gradient can be used for maximization of the
#' log-likelihood (see [RVineGrad()]; default: `grad = FALSE`).
#' @param hessian Logical; whether the Hessian matrix of parameter estimates is
#' estimated (default: `hessian = FALSE`). Note that this is not the
#' Hessian Matrix calculated via [RVineHessian()] but via finite
#' differences.
#' @param se Logical; whether standard errors of parameter estimates are
#' estimated on the basis of the Hessian matrix (see above; default: `se =
#' FALSE`).
#' @param ... Further arguments for `optim` (e.g. `factr` controls
#' the convergence of the "L-BFGS-B" method, or `trace`, a non-negative
#' integer, determines if tracing information on the progress of the
#' optimization is produced.) \cr For more details see the documentation of
#' [optim()].
#'
#' @return \item{RVM}{[RVineMatrix()] object with the calculated
#' parameters stored in `RVM$par` and `RVM$par2`. Additional
#' information about the fit is added (e.g., log-likelihood, AIC, BIC).}
#' \item{value}{Optimized log-likelihood value corresponding to the estimated
#' pair-copula parameters.} \item{convergence}{An integer code indicating
#' either successful convergence (`convergence = 0`) or an error:\cr
#' `1` = the iteration limit `maxit` has been reached \cr `51` =
#' a warning from the "L-BFGS-B" method; see component `message` for
#' further details \cr `52` = an error from the "L-BFGS-B" method; see
#' component `message` for further details} \item{message}{A character
#' string giving any additional information returned by [optim()], or
#' `NULL`.} \item{counts}{A two-element integer vector giving the number
#' of calls to `fn` and `gr` respectively. This excludes those calls
#' needed to compute the Hessian, if requested, and any calls to `fn` to
#' compute a finite-difference approximation to the gradient.}
#' \item{hessian}{If `hessian = TRUE`, the Hessian matrix is returned. Its
#' calculation is on the basis of finite differences (output of `optim`).}
#'
#' @note `RVineMLE` uses the L-BFGS-B method for optimization. \cr If the
#' analytical gradient is used for maximization, computations may be up to 10
#' times faster than using finite differences.
#'
#' @author Ulf Schepsmeier, Jeffrey Dissmann
#'
#' @seealso [RVineSeqEst()],
#' [RVineStructureSelect()],
#' [RVineMatrix()],
#' [RVineGrad()],
#' [RVineHessian()]
#'
#' @references Dissmann, J. F., E. C. Brechmann, C. Czado, and D. Kurowicka
#' (2013). Selecting and estimating regular vine copulae and application to
#' financial returns. Computational Statistics & Data Analysis, 59 (1), 52-69.
#'
#' Stoeber, J. and U. Schepsmeier (2013). Estimating standard errors in regular
#' vine copula models. Computational Statistics, 1-29
#' <https://link.springer.com/article/10.1007/s00180-013-0423-8#>.
#'
#' @examples
#'
#' # define 5-dimensional R-vine tree structure matrix
#' Matrix <- c(5, 2, 3, 1, 4,
#' 0, 2, 3, 4, 1,
#' 0, 0, 3, 4, 1,
#' 0, 0, 0, 4, 1,
#' 0, 0, 0, 0, 1)
#' Matrix <- matrix(Matrix, 5, 5)
#'
#' # define R-vine pair-copula family matrix
#' family <- c(0, 1, 3, 4, 4,
#' 0, 0, 3, 4, 1,
#' 0, 0, 0, 4, 1,
#' 0, 0, 0, 0, 3,
#' 0, 0, 0, 0, 0)
#' family <- matrix(family, 5, 5)
#'
#' # define R-vine pair-copula parameter matrix
#' par <- c(0, 0.2, 0.9, 1.5, 3.9,
#' 0, 0, 1.1, 1.6, 0.9,
#' 0, 0, 0, 1.9, 0.5,
#' 0, 0, 0, 0, 4.8,
#' 0, 0, 0, 0, 0)
#' par <- matrix(par, 5, 5)
#'
#' # define second R-vine pair-copula parameter matrix
#' par2 <- matrix(0, 5, 5)
#'
#' # define RVineMatrix object
#' RVM <- RVineMatrix(Matrix = Matrix, family = family,
#' par = par, par2 = par2,
#' names = c("V1", "V2", "V3", "V4", "V5"))
#'
#' # simulate a sample of size 300 from the R-vine copula model
#' set.seed(123)
#' simdata <- RVineSim(300, RVM)
#'
#' # compute the MLE
#' mle <- RVineMLE(simdata, RVM, grad = TRUE, trace = 0)
#'
#' # compare parameters
#' round(mle$RVM$par - RVM$par, 2)
#'
RVineMLE <- function(data, RVM, start = RVM$par, start2 = RVM$par2, maxit = 200, max.df = 30,
max.BB = list(BB1 = c(5, 6), BB6 = c(6, 6), BB7 = c(5, 6), BB8 = c(6, 1)),
grad = FALSE, hessian = FALSE, se = FALSE, ...) {
## preprocessing of arguments
args <- preproc(c(as.list(environment()), call = match.call()),
check_data,
remove_nas,
check_nobs,
check_if_01,
check_est_pars,
check_RVMs,
prep_RVMs,
na.txt = " Only complete observations are used.")
list2env(args, environment())
if (all(RVM$family == 0)) {
warning("RVM contains only Independence copulas, RVineMLE does nothing")
return(RVM)
}
## sanity checks
if (maxit <= 0)
stop("'maxit' has to be greater than zero.")
## sanity checks for start parameters
Matrix <- RVM$Matrix
sel <- lower.tri(Matrix)
if (!all(start %in% c(0, NA))) {
BiCopCheck(RVM$family[sel], start[sel], start2[sel])
if (grad == TRUE) {
if (any(RVM$family %in% c(7:10, 17:20, 27:30, 37:40, 104, 114, 124, 134, 204, 214, 224, 234))) {
message("The combination 'grad=TRUE' and a copula family ",
"of the vector (7:10,17:20,27:30,37:40,104,114,124,134,204,214,224,234) ",
"is not possible. The algorithm will continue with 'grad=FALSE'.")
grad <- FALSE
}
}
}
d <- dim(data)[2]
T <- dim(data)[1]
n <- d
N <- T
## normalization of R-vine matrix
o <- diag(RVM$Matrix)
oldRVM <- RVM
RVM <- normalizeRVineMatrix(RVM)
if (is.null(colnames(data)))
colnames(data) <- paste0("V", 1:ncol(data))
cnms <- colnames(data)
data <- data[, o[length(o):1]]
n <- dim(RVM)
N <- dim(data)[1]
## sequential estimation of start parameters if not provided
if (all(start == 0)) {
est_start <- RVineSeqEst(data, RVM, max.df = max.df, max.BB = max.BB)
start <- est_start$par
start2 <- est_start$par2
}
## Position of parameters in R-vine matrix
posParams <- (RVM$family > 0)
posParams2 <- (RVM$family %in% c(2, 7, 8, 9, 10, 17, 18, 19, 20,
27, 28, 29, 30, 37, 38, 39, 40,
104, 114, 124, 134, 204, 214, 224, 234))
posParams[is.na(posParams)] <- FALSE
posParams2[is.na(posParams2)] <- FALSE
## number of parameters
nParams <- sum(posParams, na.rm = TRUE)
nParams2 <- sum(posParams2, na.rm = TRUE)
## vectors of start parameters and corresponding pair-copula families
startpar <- double(nParams + nParams2)
Copula.Types <- RVM$family[posParams]
startpar[1:nParams] <- start[posParams]
startpar[nParams + seq_len(nParams2)] <- start2[posParams2]
## lower and upper bounds
lb <- double(nParams + nParams2)
ub <- double(nParams + nParams2)
for (i in seq_len(nParams)) {
if (Copula.Types[i] %in% c(1, 2)) {
# Normal
lb[i] <- -0.99
ub[i] <- 0.99
} else if (Copula.Types[i] %in% c(3, 13)) {
# clayton
lb[i] <- 1e-04
ub[i] <- 28
} else if (Copula.Types[i] %in% c(23, 33)) {
# rotated clayton
lb[i] <- -28
ub[i] <- -1e-04
} else if (Copula.Types[i] %in% c(4, 14)) {
# gumbel
lb[i] <- 1.0001
ub[i] <- 17
} else if (Copula.Types[i] %in% c(24, 34)) {
# rotated gumbel
lb[i] <- -17
ub[i] <- -1.0001
} else if (Copula.Types[i] == 5) {
# frank
lb[i] <- -100
ub[i] <- 100
} else if (Copula.Types[i] %in% c(6, 16)) {
# joe
lb[i] <- 1.0001
ub[i] <- 30
} else if (Copula.Types[i] %in% c(26, 36)) {
# rotated joe
lb[i] <- -30
ub[i] <- -1.0001
} else if (Copula.Types[i] %in% c(7, 17)) {
# bb1
lb[i] <- 0.001
ub[i] <- max.BB$BB1[1]
} else if (Copula.Types[i] %in% c(8, 18)) {
# bb6
lb[i] <- 1.001
ub[i] <- max.BB$BB6[1]
} else if (Copula.Types[i] %in% c(9, 19)) {
# bb7
lb[i] <- 1.001
ub[i] <- max.BB$BB7[1]
} else if (Copula.Types[i] %in% c(10, 20)) {
# bb8
lb[i] <- 1.001
ub[i] <- max.BB$BB8[1]
} else if (Copula.Types[i] %in% c(27, 37)) {
# rotated bb1
lb[i] <- -max.BB$BB1[1]
ub[i] <- -0.001
} else if (Copula.Types[i] %in% c(28, 38)) {
# rotated bb6
lb[i] <- -max.BB$BB6[1]
ub[i] <- -1.001
} else if (Copula.Types[i] %in% c(29, 39)) {
# rotated bb7
lb[i] <- -max.BB$BB7[1]
ub[i] <- -1.001
} else if (Copula.Types[i] %in% c(30, 40)) {
# rotated bb8
lb[i] <- -max.BB$BB8[1]
ub[i] <- -1.001
} else if (Copula.Types[i] %in% c(43, 44)) {
# double Clayton, Gumbel
lb[i] <- -0.9999
ub[i] <- 0.9999
} else if (Copula.Types[i] %in% c(104, 114, 204, 214)) {
# Tawn
lb[i] <- 1.001
ub[i] <- 20
} else if (Copula.Types[i] %in% c(124, 134, 224, 234)) {
# Tawn
lb[i] <- -20
ub[i] <- -1.001
}
}
todo <- which(Copula.Types %in% c(2, 7:10, 17:20, 27:30, 37:40, 104, 114, 124, 134, 204, 214, 224, 234))
for (i in seq_len(nParams2)) {
if (Copula.Types[todo[i]] == 2) {
# t
lb[nParams + i] <- 2.001
ub[nParams + i] <- max.df
} else if (Copula.Types[todo[i]] %in% c(7, 17)) {
# bb1
lb[nParams + i] <- 1.001
ub[nParams + i] <- max.BB$BB1[2]
} else if (Copula.Types[todo[i]] %in% c(8, 18)) {
# bb6
lb[nParams + i] <- 1.001
ub[nParams + i] <- max.BB$BB6[2]
} else if (Copula.Types[todo[i]] %in% c(9, 19)) {
# bb7
lb[nParams + i] <- 0.001
ub[nParams + i] <- max.BB$BB7[2]
} else if (Copula.Types[todo[i]] %in% c(10, 20)) {
# bb8
lb[nParams + i] <- 0.001
ub[nParams + i] <- max.BB$BB8[2]
} else if (Copula.Types[todo[i]] %in% c(27, 37)) {
# rotated bb1
lb[nParams + i] <- -max.BB$BB1[2]
ub[nParams + i] <- -1.001
} else if (Copula.Types[todo[i]] %in% c(28, 38)) {
# rotated bb6
lb[nParams + i] <- -max.BB$BB6[2]
ub[nParams + i] <- -1.001
} else if (Copula.Types[todo[i]] %in% c(29, 39)) {
# rotated bb7
lb[nParams + i] <- -max.BB$BB7[2]
ub[nParams + i] <- -0.001
} else if (Copula.Types[todo[i]] %in% c(30, 40)) {
# rotated bb8
lb[nParams + i] <- -max.BB$BB8[2]
ub[nParams + i] <- -0.001
} else if (Copula.Types[todo[i]] %in% c(104, 114, 124, 134, 204, 214, 224, 234)) {
# Tawn
lb[nParams + i] <- 0.001
ub[nParams + i] <- 0.99
}
}
## log-likelihood function to be maximized
optim_LL <- function(parm, data, posParams, posParams2, Copula.Types, start, start2, RVM, calcupdate = NA) {
nParams <- sum(posParams, na.rm = TRUE)
nParams2 <- sum(posParams2, na.rm = TRUE)
matrixParams <- start
matrixParams2 <- start2
matrixParams[posParams] <- parm[1:nParams]
matrixParams2[posParams2] <- parm[nParams + seq_len(nParams2)]
ll <- RVineLogLik(data, RVM, par = matrixParams, par2 = matrixParams2)
if (is.finite(ll$loglik)) {
return(ll$loglik)
} else {
if (is.na(ll$loglik)) {
message(parm)
}
message(ll$loglik)
return(-10^305)
}
}
## gradient
gradient <- function(parm, data, posParams, posParams2, Copula.Types, start, start2, RVM, calcupdate) {
nParams <- sum(posParams, na.rm = TRUE)
nParams2 <- sum(posParams2, na.rm = TRUE)
matrixParams <- start
matrixParams2 <- start2
matrixParams[posParams] <- parm[1:nParams]
if (nParams2 > 0) {
matrixParams2[posParams2] <- parm[nParams + seq_len(nParams2)]
}
grad <- RVineGrad(data = data,
RVM = RVM,
par = matrixParams,
par2 = matrixParams2,
posParams = posParams)$gradient
return(grad)
}
## default values for parscale (see optim)
pscale <- numeric()
for (i in seq_len(nParams)) {
pscale[i] <- ifelse(Copula.Types[i] %in% c(1, 2, 43, 44), 0.01, 1)
}
pscale2 <- numeric()
for (i in seq_len(nParams2)) {
pscale2[i] <- ifelse(Copula.Types[i] %in% c(104, 114, 124, 134, 204, 214, 224, 234), 0.05, 1)
}
pscale <- c(pscale, pscale2)
## (default) values for control parameters of optim
ctrl <- list(fnscale = -1,
maxit = maxit,
trace = 1,
parscale = pscale,
factr = 1e+08)
ctrl <- modifyList(ctrl, list(...))
## optimization
if (all(Copula.Types %in% c(0, 1, 2, 3:6, 13, 14, 16, 23, 24, 26, 33, 34, 36, 43, 44)) && grad == TRUE) {
if (hessian == TRUE || se == TRUE) {
out1 <- optim(par = startpar,
fn = optim_LL,
gr = gradient,
data = data,
posParams = posParams,
posParams2 = posParams2, Copula.Types = Copula.Types,
start = start,
start2 = start2,
RVM = RVM,
method = "L-BFGS-B",
control = ctrl,
lower = lb,
upper = ub,
hessian = TRUE)
} else {
out1 <- optim(par = startpar,
fn = optim_LL,
gr = gradient,
data = data,
posParams = posParams,
posParams2 = posParams2,
Copula.Types = Copula.Types,
start = start,
start2 = start2,
RVM = RVM,
method = "L-BFGS-B",
control = ctrl,
lower = lb,
upper = ub)
}
} else {
if (hessian == TRUE || se == TRUE) {
out1 <- optim(par = startpar,
fn = optim_LL,
data = data,
posParams = posParams,
posParams2 = posParams2,
Copula.Types = Copula.Types,
start = start,
start2 = start2,
RVM = RVM,
calcupdate = NA,
method = "L-BFGS-B",
control = ctrl,
lower = lb,
upper = ub,
hessian = TRUE)
} else {
out1 <- optim(par = startpar,
fn = optim_LL,
data = data,
posParams = posParams,
posParams2 = posParams2,
Copula.Types = Copula.Types,
start = start,
start2 = start2,
RVM = RVM,
calcupdate = NA,
method = "L-BFGS-B",
control = ctrl,
lower = lb,
upper = ub)
}
}
# list for final output
out <- list()
# store results in out
out$value <- out1$value
out$convergence <- out1$convergence
out$message <- out1$message
out$counts <- out1$counts
if (hessian == TRUE)
out$hessian <- out1$hessian
if (se == TRUE)
out1$se <- sqrt((diag(solve(-out1$hessian))))
# create parameter matrices
kk <- 1
for (ll in seq_len(nParams)) {
out1$par[ll] <- out1$par[ll]
if (Copula.Types[ll] %in% c(2, 7:10, 17:20, 27:30, 37:40,
104, 114, 124, 134, 204, 214, 224, 234)) {
out1$par[nParams + kk] <- out1$par[nParams + kk]
kk <- kk + 1
}
}
newpar <- newpar2 <- matrix(0, d, d)
newpar[posParams] <- out1$par[1:nParams]
newpar2[posParams2] <- out1$par[nParams + seq_len(nParams2)]
newpar[(oldRVM$family == 5) & (newpar == 0)] <- 1e-10
# create RVineMatrix object
out$RVM <- RVineMatrix(Matrix = oldRVM$Matrix,
family = oldRVM$family,
par = newpar,
par2 = newpar2,
names = oldRVM$names)
# add standad errors
if (se == TRUE) {
out$RVM$se <- matrix(0, d, d)
out$RVM$se2 <- matrix(0, d, d)
out$RVM$se[posParams] <- out1$se[1:nParams]
out$RVM$se2[posParams2] <- out1$se[nParams + seq_len(nParams2)]
}
# add summary statistics
like <- RVineLogLik(data[, cnms], out$RVM)
out$RVM$logLik <- like$loglik
out$RVM$pair.logLik <- like$V$value
npar <- sum(out$RVM$family %in% allfams[onepar], na.rm = TRUE) +
2 * sum(out$RVM$family %in% allfams[twopar], na.rm = TRUE)
npar_pair <- (out$RVM$family %in% allfams[onepar]) +
2 * (out$RVM$family %in% allfams[twopar])
out$RVM$AIC <- -2 * like$loglik + 2 * npar
out$RVM$pair.AIC <- -2 * like$V$value + 2 * npar_pair
out$RVM$BIC <- -2 * like$loglik + log(N) * npar
out$RVM$pair.BIC <- -2 * like$V$value + log(N) * npar_pair
out$RVM$emptau <- oldRVM$emptau
out$RVM$p.value.indeptest <- oldRVM$p.value.indeptest
## return results
out
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.