R: Selection and Maximum penalized likelihood estimation of a...

gamBiCopSelect {gamCopula}

R Documentation

Selection and Maximum penalized likelihood estimation of a Generalized Additive model (gam) for the copula parameter or Kendall's tau.

Description

This function selects an appropriate bivariate copula family for given bivariate copula data using one of a range of methods. The corresponding parameter estimates are obtained by maximum penalized likelihood estimation, where each Newton-Raphson iteration is reformulated as a generalized ridge regression solved using the mgcv package.

Usage

gamBiCopSelect(
  udata,
  lin.covs = NULL,
  smooth.covs = NULL,
  familyset = NA,
  rotations = TRUE,
  familycrit = "AIC",
  level = 0.05,
  edf = 1.5,
  tau = TRUE,
  method = "FS",
  tol.rel = 0.001,
  n.iters = 10,
  parallel = FALSE,
  verbose = FALSE,
  select.once = TRUE,
  ...
)

Arguments

`udata`	A matrix or data frame containing the model responses, (u1,u2) in [0,1]x[0,1]
`lin.covs`	A matrix or data frame containing the parametric (i.e., linear) covariates.
`smooth.covs`	A matrix or data frame containing the non-parametric (i.e., smooth) covariates.
`familyset`	(Similar to `BiCopSelect` from the `VineCopula` package) Vector of bivariate copula families to select from. If `familyset = NA` (default), selection among all possible families is performed. Coding of bivariate copula families: `1` Gaussian, `2` Student t, `5` Frank, `301` Double Clayton type I (standard and rotated 90 degrees), `302` Double Clayton type II (standard and rotated 270 degrees), `303` Double Clayton type III (survival and rotated 90 degrees), `304` Double Clayton type IV (survival and rotated 270 degrees), `401` Double Gumbel type I (standard and rotated 90 degrees), `402` Double Gumbel type II (standard and rotated 270 degrees), `403` Double Gumbel type III (survival and rotated 90 degrees), `404` Double Gumbel type IV (survival and rotated 270 degrees).
`rotations`	If `TRUE`, all rotations of the families in familyset are included.
`familycrit`	Character indicating the criterion for bivariate copula selection. Possible choices: `familycrit = 'AIC'` (default) or `'BIC'`, as in `BiCopSelect` from the `VineCopula` package.
`level`	Numerical; significance level of the test for removing individual predictors (default: `level = 0.05`).
`edf`	Numerical; if the estimated EDF for individual predictors is smaller than `edf` but the predictor is still significant, then it is set as linear (default: `edf = 1.5`).
`tau`	`FALSE` for a calibration function specified for the Copula parameter or `TRUE` (default) for a calibration function specified for Kendall's tau.
`method`	`'FS'` for Fisher-scoring (default) and `'NR'` for Newton-Raphson.
`tol.rel`	Relative tolerance for `'FS'`/`'NR'` algorithm.
`n.iters`	Maximal number of iterations for `'FS'`/`'NR'` algorithm.
`parallel`	`TRUE` for a parallel estimation across copula families.
`verbose`	`TRUE` prints informations during the estimation.
`select.once`	if `TRUE` the GAM structure is only selected once, for the family that appears first in `familyset`.
`...`	Additional parameters to be passed to `gam`

Value

gamBiCopFit returns a list consisting of

`res`	S4 `gamBiCop-class` object.
`method`	`'FS'` for Fisher-scoring and `'NR'` for Newton-Raphson.
`tol.rel`	relative tolerance for `'FS'`/`'NR'` algorithm.
`n.iters`	maximal number of iterations for `'FS'`/`'NR'` algorithm.
`trace`	the estimation procedure's trace.
`conv`	`0` if the algorithm converged and `1` otherwise.

Examples

require(copula)
set.seed(0)

## Simulation parameters (sample size, correlation between covariates,
## Student copula with 4 degrees of freedom)
n <- 5e2
rho <- 0.9
fam <- 2
par2 <- 4

## A calibration surface depending on four variables
eta0 <- 1
calib.surf <- list(
  calib.lin <- function(t, Ti = 0, Tf = 1, b = 2) {
    return(-2 + 4 * t)
  },
  calib.quad <- function(t, Ti = 0, Tf = 1, b = 8) {
    Tm <- (Tf - Ti) / 2
    a <- -(b / 3) * (Tf^2 - 3 * Tf * Tm + 3 * Tm^2)
    return(a + b * (t - Tm)^2)
  },
  calib.sin <- function(t, Ti = 0, Tf = 1, b = 1, f = 1) {
    a <- b * (1 - 2 * Tf * pi / (f * Tf * pi +
      cos(2 * f * pi * (Tf - Ti))
      - cos(2 * f * pi * Ti)))
    return((a + b) / 2 + (b - a) * sin(2 * f * pi * (t - Ti)) / 2)
  },
  calib.exp <- function(t, Ti = 0, Tf = 1, b = 2, s = Tf / 8) {
    Tm <- (Tf - Ti) / 2
    a <- (b * s * sqrt(2 * pi) / Tf) * (pnorm(0, Tm, s) - pnorm(Tf, Tm, s))
    return(a + b * exp(-(t - Tm)^2 / (2 * s^2)))
  }
)

## 6-dimensional matrix X of covariates
covariates.distr <- mvdc(normalCopula(rho, dim = 6),
  c("unif"), list(list(min = 0, max = 1)),
  marginsIdentical = TRUE
)
X <- rMvdc(n, covariates.distr)
colnames(X) <- paste("x", 1:6, sep = "")

## U in [0,1]x[0,1] depending on the four first columns of X
U <- condBiCopSim(fam, function(x1, x2, x3, x4) {
  eta0 + sum(mapply(function(f, x)
    f(x), calib.surf, c(x1, x2, x3, x4)))
}, X[, 1:4], par2 = 4, return.par = TRUE)
## Not run: 
## Selection using AIC (about 30sec on single core)
## Use parallel = TRUE to speed-up....
system.time(best <- gamBiCopSelect(U$data, smooth.covs = X))
print(best$res)
EDF(best$res) ## The first function is linear
## Plot only the smooth component
par(mfrow = c(2, 2))
plot(best$res)

## End(Not run)