In lizhengxiao/rMGLReg: MGL copula regressions

Case I - modelling Danish fire insurance data with covariates

This example implements the MGL and MGL-EV copula regression models of Zhengxiao Li et al. (2021, “MGL copulas”).

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

read the danish data set

  library(rMGLReg)
  library(fitdistrplus)
  library(splines)
  library(snpar)
  data("danishmulti")
  dt <- data.table::data.table(danishmulti)
  dtnew <- dt[Building>0&Contents>0]
  y1 <- dtnew$Building
  y2 <- dtnew$Contents
  y <- cbind(y1, y2)

  # empirical cdf
  u1 <- snpar::kde(y[,1], kernel = "gaus", 
             xgrid = y[,1],
             h = 0.2)$Fhat
  u2 <- snpar::kde(y[,2], kernel = "gaus", 
             xgrid = y[,2],
             h = 0.2)$Fhat
  U <- cbind(u1, u2) # two-dim

library(ggplot2)
library(latex2exp)
Usample <- U
XY <- y
newtheme <-   theme_bw() + theme(axis.line = element_line(colour = "black"),
                                 axis.text.x = element_text(margin = margin(t = 0.25, unit = "cm")),
                                 axis.text.y = element_text(margin = margin(r = 0.25, unit = "cm"),
                                                            size = 10, 
                                                            vjust = 0.5, 
                                                            hjust = 0.5),
                                 axis.ticks.length = unit(-0.1, "cm"),
                                 plot.title = element_text(hjust = 0.5),
                                 legend.direction = 'vertical',
                                 legend.position = c(.95, .99),
                                 legend.justification = c("right", "top"),
                                 legend.box.just = "right",
                                 legend.text = element_text(size = 10)) 

dtplot <- data.frame(U1 = Usample[,1], U2 = Usample[,2],
                     logY1 = log(XY[,1]), logY2 = log(XY[,2]))

p1 <- ggplot(data = dtplot, mapping = aes(y = logY1, 
                                          x = logY2)) + 
  newtheme +
  geom_point() + 
  labs(title = "", 
       x = TeX("$log Y_1$"), 
       y = TeX("$log Y_2$")) +  
  scale_x_continuous(limits = c(-7.5, 5),
                     breaks = seq(-7.5, 5, by = 2.5)) +
  scale_y_continuous(limits = c(-5, 5),
                     breaks = seq(-5, 5, by = 2.5))

p1

p2 <- ggplot(data = dtplot, mapping = aes(y = U2, 
                                           x = U1)) + 
  newtheme +
  geom_point() + 
  labs(title = "", 
     x = TeX("$u_1$"), 
     y = TeX("$u_2$")) +  
  scale_x_continuous(limits = c(0, 1),
                     breaks = seq(0, 1, by = 0.2)) +
  scale_y_continuous(limits = c(0, 1),
                     breaks = seq(0, 1, by = 0.2))

p2
library(patchwork)
p0 <- p1 + p2 + plot_layout(ncol = 2)
p0

survival MGL copual and survival MGL-EV regression

  dtnew[, year := as.numeric(substr(Date, start = 1, stop = 4))]
  dtnew[, x1 := (year - mean(year))/(sd(year))]
  dtnew[, day := difftime(as.Date(Date), as.Date("1980-01-03"), units="days")]

  X <- splines::ns(dtnew$year, knots = quantile(dtnew$year, c(0.5)), intercept = T)

  # survival MGL copual regression
  m.MGL180 <- MGL.reg(U = U, copula = "MGL180",
                                 X = X, method = "Nelder-Mead",
                                 initpar = c(-0.32, 0.001, 0.001)
  )
  # survival MGL-EV copual regression
  m.MGLEV180 <- MGL.reg(U = U, copula = "MGL-EV180",
                         method = "Nelder-Mead",
                         X = X,
                         initpar = c(-0.32, 0.001, 0.001)
  )
  # gumbel regression
  m.gumbel <- MGL.reg(U = U, copula = "Gumbel",
                         method = "Nelder-Mead",
                         X = X,
                         initpar = c(-0.32, 0.001, 0.001)
  )
  # MGB2 regression
  m.MGB2 <- MGL.reg(U = U, copula = "MGB2",
                         method = "Nelder-Mead",
                         X = X,
                         initpar = c(0.1, -0.5, -0.5, -0.5, -0.5)
  )

Estimates and goodness fit

  estimates.copula <- rbind(c(m.MGL180$estimates, NA, NA), 
                            c(m.MGLEV180$estimates, NA, NA), 
                            c(m.gumbel$estimates, NA, NA),
                            c(m.MGB2$estimates)
                            )
  sd.copula <- rbind(c(m.MGL180$se, NA, NA), 
                     c(m.MGLEV180$se, NA, NA), 
                     c(m.gumbel$se, NA, NA),
                     c(m.MGB2$se)
                     )
  ll.copula <- rbind((m.MGL180)$loglike, (m.MGLEV180)$loglike, (m.gumbel)$loglike, (m.MGB2)$loglike)
  AIC.copula <- rbind((m.MGL180)$AIC, (m.MGLEV180)$AIC, (m.gumbel)$AIC, (m.MGB2)$AIC)
  BIC.copula <- rbind((m.MGL180)$BIC, (m.MGLEV180)$BIC, (m.gumbel)$BIC, (m.MGB2)$BIC)
  table1 <- cbind(estimates.copula, sd.copula, ll.copula, AIC.copula, BIC.copula)

  row.names(table1) <- c('MGL180', 'MGLEV180', 'Gumbel', 'MGB2') 
  colnames(table1) <- c("estimates", "estimates","estimates", "estimates","estimates", 
                        "se", "se","se","se","se",
                        'loglike', 
                        'AIC', 'BIC')
  knitr::kable(t(table1), digits = 3)

The predicted value of copula parameter with different value of the covariate \textit{Year} for the Danish fire insurance data set.

# Survival MGL regression model
par.est <- m.MGL180$estimates
agepred <- seq(from = 1980, to = 1990)
Xcopula <-  splines::ns(agepred, knots = quantile(dtnew$year, c(0.5)), intercept = T)
delta.est <- exp(Xcopula%*%par.est)
cov.est <- -solve(m.MGL180$hessian)
beta.sim <- matrix(0, nrow = 100, ncol = 3)
delta.mat <- matrix(0, nrow = length(agepred), ncol = 100)
beta.sim <- mvrnorm(100, mu = par.est, Sigma = cov.est)
for(i in 1:100){
  delta.mat[,i] <- exp(Xcopula%*%beta.sim[i,])
}
matplot(delta.mat, col = 'gray', type = 'l', ylim = c(0, 2), main = 'Survival MGL')
lines(delta.est, col = 'red', lwd = 2)


# Surivial MGL-EV regression
par.est <- m.MGLEV180$estimates
delta.est <- exp(Xcopula%*%par.est)
cov.est <- -solve(m.MGLEV180$hessian)
beta.sim <- matrix(0, nrow = 100, ncol = 3)
delta.mat <- matrix(0, nrow = length(agepred), ncol = 100)
beta.sim <- mvrnorm(100, mu = par.est, Sigma = cov.est)
for(i in 1:100){
  delta.mat[,i] <- exp(Xcopula%*%beta.sim[i,])
}
matplot(delta.mat, col = 'gray', type = 'l', ylim = c(0, 2), main = 'Survival MGL-EV')
lines(delta.est, col = 'red', lwd = 2)

# Gumbel regression
par.est <- m.gumbel$estimates
delta.est <- exp(Xcopula%*%par.est) + 1

cov.est <- -solve(m.gumbel$hessian)
beta.sim <- matrix(0, nrow = 100, ncol = 3)
delta.mat <- matrix(0, nrow = length(agepred), ncol = 100)
beta.sim <- mvrnorm(100, mu = par.est, Sigma = cov.est)
for(i in 1:100){
  delta.mat[,i] <- exp(Xcopula%*%beta.sim[i,]) + 1
}
matplot(delta.mat, col = 'gray', type = 'l', ylim = c(0, 2), main = 'Gumbel')
lines(delta.est, col = 'red', lwd = 2)

# MGB2 regression
par.est <- m.MGB2$estimates
delta.est <- exp(Xcopula%*%par.est[1:ncol(X)])

cov.est <- -solve(m.MGB2$hessian)
beta.sim <- matrix(0, nrow = 100, ncol = 5)
delta.mat <- matrix(0, nrow = length(agepred), ncol = 100)
pars.sim <- mvrnorm(100, mu = par.est, Sigma = cov.est)
beta.sim <- pars.sim[,1:ncol(X)]
for(i in 1:100){
  delta.mat[,i] <- exp(Xcopula%*%beta.sim[i,])
}
matplot(delta.mat, col = 'gray', type = 'l', ylim = c(0, 4), main = 'MGB2')
lines(delta.est, col = 'red', lwd = 2)

# Survival MGL regression model
library(reshape2)
par.est <- m.MGL180$estimates
agepred <- seq(from = 1980, to = 1990)
Xcopula <-  splines::ns(agepred, knots = quantile(dtnew$year, c(0.5)), intercept = T)
delta.est <- exp(Xcopula%*%par.est)
cov.est <- -solve(m.MGL180$hessian)
beta.sim <- matrix(0, nrow = 100, ncol = 3)
delta.mat <- matrix(0, nrow = length(agepred), ncol = 100)
beta.sim <- mvrnorm(100, mu = par.est, Sigma = cov.est)
for(i in 1:100){
  delta.mat[,i] <- exp(Xcopula%*%beta.sim[i,])
}
data <- as.data.frame(delta.mat)
# data <- data.table(data)
#id variable for position in matrix
data$id <- agepred
#reshape to long format
plot_data <- melt(data, id.var = "id")
p1  <- ggplot() +
  theme_bw() + xlab('Year') + ylab(expression(paste('Predicted', "  ",delta))) +
  geom_line(data = plot_data, mapping = aes(y = value,
                                            x = id,
                                            group = variable),
            size = 0.3, linetype = 1, col = 'gray') +
  geom_path(aes_(x = agepred, y = delta.est),
            size = 1, col = 'red'
  ) +
  scale_x_continuous(limits = c(1980, 1990),
                     breaks = seq(1980, 1990, by = 2)) +
  scale_y_continuous(limits = c(0, 2),
                     breaks = seq(0, 2, by = 0.5)) +
  # ggtitle("Survival MGL copula regression") +
  facet_grid(. ~ "Survival MGL copula regression") +
  theme(plot.title = element_text(hjust = 0.5))

p1



# Surivial MGL-EV regression
par.est <- m.MGLEV180$estimates
agepred <- seq(from = 1980, to = 1990)
Xcopula <-  ns(agepred, knots = quantile(dtnew$year, c(0.5)), intercept = T)
delta.est <- exp(Xcopula%*%par.est)
cov.est <- -solve(m.MGLEV180$hessian)
beta.sim <- matrix(0, nrow = 100, ncol = 3)
delta.mat <- matrix(0, nrow = length(agepred), ncol = 100)
beta.sim <- mvrnorm(100, mu = par.est, Sigma = cov.est)
for(i in 1:100){
  delta.mat[,i] <- exp(Xcopula%*%beta.sim[i,])
}

data <- as.data.frame(delta.mat)
#id variable for position in matrix
data$id <- agepred
#reshape to long format
plot_data <- melt(data, id.var = "id")
p2  <- ggplot() +
  theme_bw() + xlab('Year') + ylab(expression(paste('Predicted', "  ",delta))) +
  geom_line(data = plot_data, mapping = aes(y = value,
                                            x = id,
                                            group = variable),
            size = 0.3, linetype = 1, col = 'gray') +
  geom_path(aes_(x = agepred, y = delta.est),
            size = 1, col = 'red'
  ) +
  scale_x_continuous(limits = c(1980, 1990),
                     breaks = seq(1980, 1990, by = 2)) +
  scale_y_continuous(limits = c(0, 2),
                     breaks = seq(0, 2, by = 0.5)) +
  # ggtitle("Survival MGL-EV copula regression") +
  facet_grid(. ~ "Survival MGL-EV copula regression") +
  theme(plot.title = element_text(hjust = 0.5))

p2


# Gumbel regression
par.est <- m.gumbel$estimates
delta.est <- exp(Xcopula%*%par.est) + 1

cov.est <- -solve(m.gumbel$hessian)
beta.sim <- matrix(0, nrow = 100, ncol = 3)
delta.mat <- matrix(0, nrow = length(agepred), ncol = 100)
beta.sim <- mvrnorm(100, mu = par.est, Sigma = cov.est)
for(i in 1:100){
  delta.mat[,i] <- exp(Xcopula%*%beta.sim[i,]) + 1
}
data <- as.data.frame(delta.mat)
#id variable for position in matrix
data$id <- agepred
#reshape to long format
plot_data <- melt(data, id.var = "id")
p3  <- ggplot() +
  theme_bw() + xlab('Year') + ylab(expression(paste('Predicted', "  ",delta))) +
  geom_line(data = plot_data, mapping = aes(y = value,
                                            x = id,
                                            group = variable),
            size = 0.3, linetype = 1, col = 'gray') +
  geom_path(aes_(x = agepred, y = delta.est),
            size = 1, col = 'red'
  ) +
  scale_x_continuous(limits = c(1980, 1990),
                     breaks = seq(1980, 1990, by = 2)) +
  scale_y_continuous(limits = c(0, 2),
                     breaks = seq(0, 2, by = 0.5)) +
  # scale_y_continuous(limits = c(1, 1.6),
  #                    breaks = seq(1, 1.6, by = 0.2)) +

  # ggtitle("Gumbel copula regression") +
  facet_grid(. ~ "Gumbel copula regression") +
  theme(plot.title = element_text(hjust = 0.5))

p3


# MGB2 regression
par.est <- m.MGB2$estimates
delta.est <- exp(Xcopula%*%par.est[1:ncol(X)])
cov.est <- -solve(m.MGB2$hessian)
beta.sim <- matrix(0, nrow = 100, ncol = 5)
delta.mat <- matrix(0, nrow = length(agepred), ncol = 100)
pars.sim <- mvrnorm(100, mu = par.est, Sigma = cov.est)
beta.sim <- pars.sim[,1:ncol(X)]
for(i in 1:100){
  delta.mat[,i] <- exp(Xcopula%*%beta.sim[i,])
}


data <- as.data.frame(delta.mat)
#id variable for position in matrix
data$id <- agepred
#reshape to long format
plot_data <- melt(data, id.var = "id")
p4  <- ggplot() +
  theme_bw() + xlab('Year') + ylab(expression(paste('Predicted', "  ", q))) +
  geom_line(data = plot_data, mapping = aes(y = value,
                                            x = id,
                                            group = variable),
            size = 0.3, linetype = 1, col = 'gray') +
  geom_path(aes_(x = agepred, y = delta.est),
            size = 1, col = 'red'
  ) +
  scale_x_continuous(limits = c(1980, 1990),
                     breaks = seq(1980, 1990, by = 2)) +
  scale_y_continuous(limits = c(0, 2),
                     breaks = seq(0, 2, by = 0.5)) +
  facet_grid(. ~ "MGB2 copula regression") +
  theme(plot.title = element_text(hjust = 0.5))

p4


p <- p1 + p2 + p3 + p4
p

# ggsave(p, file = 'C:/Users/Lee/Desktop/pic_name.pdf', width = 12, height = 8)