Sequential maximum penalized likelihood estimation of a GAM-Vine model.

Source: R/gamVineSeqFit.R

This function estimates the parameter(s) of a Generalized Additive model (GAM) Vine model, where GAMs for individual edges are specified either for the copula parameter or Kendall's tau. It solves the maximum penalized likelihood estimation for the copula families supported in this package by reformulating each Newton-Raphson iteration as a generalized ridge regression, which is solved using the mgcv package.

gamVineSeqFit(
  data,
  GVC,
  covariates = NA,
  method = "FS",
  tol.rel = 0.001,
  n.iters = 10,
  verbose = FALSE,
  ...
)

Arguments

data

A matrix or data frame containing the data in [0,1]^d.

GVC

A gamVine object.

covariates

Vector of names for the covariates.

method

'NR' for Newton-Raphson and 'FS' for Fisher-scoring (default).

tol.rel

Relative tolerance for 'FS'/'NR' algorithm.

n.iters

Maximal number of iterations for 'FS'/'NR' algorithm.

verbose

TRUE if informations should be printed during the estimation and FALSE (default) for a silent version.

...

Additional parameters to be passed to gam from mgcv.

Value

gamVineSeqFit returns a gamVine object.

Examples

require(mgcv)
set.seed(0)

##  Simulation parameters
# Sample size
n <- 1e3
# Copula families
familyset <- c(1:2, 301:304, 401:404)
# Define a 4-dimensional R-vine tree structure matrix
d <- 4
Matrix <- c(2, 3, 4, 1, 0, 3, 4, 1, 0, 0, 4, 1, 0, 0, 0, 1)
Matrix <- matrix(Matrix, d, d)
nnames <- paste("X", 1:d, sep = "")

## A function factory
eta0 <- 1
calib.surf <- list(
  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)))
  }
)

##  Create the model
# Define gam-vine model list
count <- 1
model <- vector(mode = "list", length = d * (d - 1) / 2)
sel <- seq(d, d^2 - d, by = d)

# First tree
for (i in 1:(d - 1)) {
  # Select a copula family
  family <- sample(familyset, 1)
  model[[count]]$family <- family

  # Use the canonical link and a randomly generated parameter
  if (is.element(family, c(1, 2))) {
    model[[count]]$par <- tanh(rnorm(1) / 2)
    if (family == 2) {
      model[[count]]$par2 <- 2 + exp(rnorm(1))
    }
  } else {
    if (is.element(family, c(401:404))) {
      rr <- rnorm(1)
      model[[count]]$par <- sign(rr) * (1 + abs(rr))
    } else {
      model[[count]]$par <- rnorm(1)
    }
    model[[count]]$par2 <- 0
  }
  count <- count + 1
}

# A dummy dataset
data <- data.frame(u1 = runif(1e2), u2 = runif(1e2), matrix(runif(1e2 * d), 1e2, d))

# Trees 2 to (d-1)
for (j in 2:(d - 1)) {
  for (i in 1:(d - j)) {
    # Select a copula family
    family <- sample(familyset, 1)

    # Select the conditiong set and create a model formula
    cond <- nnames[sort(Matrix[(d - j + 2):d, i])]
    tmpform <- paste("~", paste(paste("s(", cond, ", k=10, bs='cr')",
      sep = ""
    ), collapse = " + "))
    l <- length(cond)
    temp <- sample(3, l, replace = TRUE)

    # Spline approximation of the true function
    m <- 1e2
    x <- matrix(seq(0, 1, length.out = m), nrow = m, ncol = 1)
    if (l != 1) {
      tmp.fct <- paste("function(x){eta0+",
        paste(sapply(1:l, function(x)
          paste("calib.surf[[", temp[x], "]](x[", x, "])",
            sep = ""
          )), collapse = "+"), "}",
        sep = ""
      )
      tmp.fct <- eval(parse(text = tmp.fct))
      x <- eval(parse(text = paste0("expand.grid(",
        paste0(rep("x", l), collapse = ","), ")",
        collapse = ""
      )))
      y <- apply(x, 1, tmp.fct)
    } else {
      tmp.fct <- function(x) eta0 + calib.surf[[temp]](x)
      colnames(x) <- cond
      y <- tmp.fct(x)
    }

    # Estimate the gam model
    form <- as.formula(paste0("y", tmpform))
    dd <- data.frame(y, x)
    names(dd) <- c("y", cond)
    b <- gam(form, data = dd)
    # plot(x[,1],(y-fitted(b))/y)

    # Create a dummy gamBiCop object
    tmp <- gamBiCopFit(data = data, formula = form, family = 1, n.iters = 1)$res

    # Update the copula family and the model coefficients
    attr(tmp, "model")$coefficients <- coefficients(b)
    attr(tmp, "model")$smooth <- b$smooth
    attr(tmp, "family") <- family
    if (family == 2) {
      attr(tmp, "par2") <- 2 + exp(rnorm(1))
    }
    model[[count]] <- tmp
    count <- count + 1
  }
}

# Create the gamVineCopula object
GVC <- gamVine(Matrix = Matrix, model = model, names = nnames)
print(GVC)
#> GAM-Vine matrix: 
#>      [,1] [,2] [,3] [,4]
#> [1,]    2    0    0    0
#> [2,]    3    3    0    0
#> [3,]    4    4    4    0
#> [4,]    1    1    1    1
#> 
#>  Where 
#> 1  <->  X1 
#> 2  <->  X2 
#> 3  <->  X3 
#> 4  <->  X4 
#> 
#>  Tree 1: 
#> X2,X1: Gumbel type 3 (survival and 90 degrees rotated) 
#> X3,X1: Gaussian 
#> X4,X1: Gumbel type 1 (standard and 90 degrees rotated) 
#> 
#>  Tree 2: 
#> X2,X4|X1 : Gumbel type 1 (standard and 90 degrees rotated) copula with tau(z) = (exp(z)-1)/(exp(z)+1) where 
#> z ~ s(X1, k = 10, bs = "cr")
#> X3,X4|X1 : Gumbel type 2 (standard and 270 degrees rotated) copula with tau(z) = (exp(z)-1)/(exp(z)+1) where 
#> z ~ s(X1, k = 10, bs = "cr")
#> 
#>  Tree 3: 
#> X2,X3|X4,X1 : Gumbel type 4 (survival and 270 degrees rotated) copula with tau(z) = (exp(z)-1)/(exp(z)+1) where 
#> z ~ s(X1, k = 10, bs = "cr") + s(X4, k = 10, bs = "cr")
if (FALSE) { # \dontrun{
## Simulate and fit the model
sim <- gamVineSimulate(n, GVC)
fitGVC <- gamVineSeqFit(sim, GVC, verbose = TRUE)
fitGVC2 <- gamVineCopSelect(sim, Matrix, verbose = TRUE)
(gamVinePDF(GVC, sim[1:10, ]))

## Plot the results
dev.off()
par(mfrow = c(3, 4))
plot(GVC, ylim = c(-2.5, 2.5))

plot(fitGVC, ylim = c(-2.5, 2.5))

plot(fitGVC2, ylim = c(-2.5, 2.5))
} # }