Notebook: fitting vine copula models

Model for simulating data to use in fitting.

The model we are using to simultate the data is based on the simple Bayesian network shown in Figure 1. We’ve decomposed this in such a way that \(C_{12}\), \(C_{13}\), \(C_{34}\), and \(C_{23;1}\)are all part of the decomposition.
These three marginal distributions would be specified in the BN, and the conditional distribution is independent by construction.

Figure 1: The Bayesian network we intend to approximate as a vine copula model.

The model also requires \(C_{24;13}\), which we take just to be equal to \(C_{24}\). The last joint conditional needed is \(C_{14;3}\). For this we assume that \(\rho_{14;3} = \rho_{12;3} \rho_{24;3} = \rho_{12}\rho_{24}\), with the latter two specified as part of the BN.

The joint distributions attached to the edges of the BN graph are (with \(G(\rho)\) the bivariate Gaussian copula with correlation \(\rho\) and \(t(\rho, \nu)\) the bivariate t-copula with correlation \(\rho\) and degrees of freedom \(\nu\)):

Edge	Dist
1 -> 2	G(0.5)
1 -> 3	t(0.5, 4)
2 -> 4	G(0.25)
3 -> 4	t(-0.25, 2)

So, notionally, 1 is activating both 2 and 3, with 2 activating 4 and 3 inhibiting 4. Since both of the joint distributions along \(1 \rightarrow 2 \rightarrow 4\) are Gaussian, we take $C_{14;3} to be Gaussian as well.

## Data for fitting.
vmat <- matrix(c(3,1,2,4, 1,2,3,0, 1,2,0,0, 1,0,0,0), ncol=4)
vstruct <- as_rvine_structure(vmat)
ci <- bicop_dist()
## Pair copula list
pclist <- list(
  list(bicop_dist('t', parameters=c(-0.25, 2)), bicop_dist('t', parameters=c(0.5, 4)), bicop_dist('gaussian', parameters=0.5)),
  list(bicop_dist('gaussian', parameters=0.125), ci),
  list(bicop_dist('gaussian', parameters=0.25))
)
vcop_true <- vinecop_dist(pclist, vstruct)

plot(vcop_true)

contour(vcop_true)

## generate samples and transform to std. normal margins.
z <- qnorm(rvinecop(500, vcop_true, qrng=TRUE, cores=2))

colnames(z) <- paste0('z', seq(1,4))
zf <- as.data.frame(z)
GGally::ggpairs(zf, progress=FALSE)

We see that \(\rho_{14} \approx 0\), which makes sense because \(Z_1\) is activating \(Z_2\) and \(Z_3\) about equally, and \(Z_2\) is activating \(Z_4\) with about the same strength that \(Z_3\) is inhibiting it. We also have $(Z_2Z_4) = $ 0.21, which is a bit lower than our intended value of 0.25, suggesting that setting \(C_{24;3} = C_{24}\) is not giving us exactly what we wanted (but maybe not; the partial correlations don’t look too bad). Finally, we can estimate the partial correlation \(\rho_{14;23}\), which we intend to be zero.

partialcor(~z2+z3, zf)

              cor          z      pval     lowerCI    upperCI
z1~z4 0.003236339 0.07193148 0.9426564 -0.08474286 0.09116547

Commentary

These results suggest that we did a pretty good job of capturing the essential characteristics of the Bayesian network in a vine copula model. In particular, \(Z_1 \perp \!\!\! \perp Z_4 |Z_2,Z_3\), and \(Z_2 \perp \!\!\! \perp Z_3 | Z_1\), which are the sine qua non of the network in Figure 1. Thus, taking \(C_{14;3}\) to be the composition of \(C_{12}\) and \(C_{24}\) seems to have been the right idea. This distribution was easy to construct from our two Gaussian copula distributions, but in principle you could do it for any combination of distributions by constructing conditional distributions and using the fact that \(p(4|1) = p(4|2) p(2|1)\). Albeit, this involves some rather tedious partial derivatives on the front end and integration on the back.

Although we were pretty successful at producing a good approximation to a simple Bayesian network using vine copulas, that was not actually our primary purpose for this exercise. What we really want to do is to try to fit a vine distribution to the data sampled from this distribution. We turn our attention to this task next.

Fitting a vine copula model

The vinecop function provides for automatic fitting and model selection for vine copula models. We have a lot of different options for what to estimate and what to hold fixed, so we will try several of them to see what kind of results we get.

We gave our observations standard normal margins, but we wouldn’t know that in practice, so we will estimate the pseudo-observations from the empirical marginal DFs.

u <- pseudo_obs(z)
colnames(u) <- paste0('u', seq(1,4))

Fitting vine copulas to the sample dataset

Structure fixed

Pair copulas from elliptical family

We will start with the simplest thing we can do, which is to give the algorithm the structure that we used and have it fit elliptical (or independence) pair copulas. I tried several variations on this method, including changing the selection criterion, using maximum likelihood instead of inverse-tau for the parameter estimation, and changing the independence prior psi0, but these made only trivial differences in the model produced.

fit1 <- vinecop(u, family_set=c('indep','elliptical'), structure=vstruct, selcrit='mbicv', 
                par_method='itau', psi0=0.5, show_trace=TRUE, cores=4)

** Tree: 0
4,3 <-> Student, parameters = -0.253297
  2.09253
3,1 <-> Student, parameters = 0.501381
 4.04004
2,1 <-> Gaussian, parameters = 0.494742
** Tree: 1
4,1 | 3 <-> Gaussian, parameters = 0.129754
3,2 | 1 <-> Independence
** Tree: 2
4,2 | 1,3 <-> Gaussian, parameters = 0.250494

fit1

4-dimensional vine copula fit ('vinecop')
nobs = 500   logLik = 223.46   npars = 7   AIC = -432.91   BIC = -403.41   

As we can see, this fit was very close to the model that we started with.

Pair copulas estimated from a broader set

Here we broaden the available set of copulas to include all families for which a Kendall’s tau inversion formula is available. This is the elliptical families, plus the single-parameter Archimedean families. Since one of our relationships has a negative \(\tau\) value, an Archimedean copula would have to be rotated to fit that data. It’s not clear from the documentation if the fitting can do this automatically.

fit2 <- vinecop(u, family_set=c('itau'), structure=vstruct, selcrit='mbicv', 
                par_method='itau', psi0=0.5, show_trace=TRUE, cores=4)

** Tree: 0
4,3 <-> Student, parameters = -0.253297
  2.09253
3,1 <-> Student, parameters = 0.501381
 4.04004
2,1 <-> Gaussian, parameters = 0.494742
** Tree: 1
4,1 | 3 <-> Frank, parameters = 0.749713
3,2 | 1 <-> Independence
** Tree: 2
4,2 | 1,3 <-> Gaussian, parameters = 0.252249

fit2

4-dimensional vine copula fit ('vinecop')
nobs = 500   logLik = 223.82   npars = 7   AIC = -433.63   BIC = -404.13   

For \(C_{14;3}\) this model went with a Frank copula, instead of the Gaussian we used in the base model. The Frank copula, like the Gaussian, is a single-parameter copula with no tail dependence. We can compare the bivariate versions graphically:

contour(bicop_dist('gaussian', parameters=0.13), main='Gaussian: rho=0.13')

contour(bicop_dist('frank', parameters=0.75), main='Frank: theta=0.75')

As we can see, the differences between the two are barely perceptible.

Pair copulas from Archimedean copulas only

The main purpose of this experiment is to see if the fitting algorithms will automatically produce a rotated copula when appropriate.

fit3 <- vinecop(u, family_set=c('indep','archimedean'), structure=vstruct, selcrit='mbicv', 
                par_method='itau', psi0=0.5, show_trace=TRUE, cores=4)

** Tree: 0
4,3 <-> Gumbel 270°, parameters = 1.19479
3,1 <-> Gumbel 180°, parameters = 1.50229
2,1 <-> Frank, parameters = 3.25965
** Tree: 1
4,1 | 3 <-> Frank, parameters = 0.730207
3,2 | 1 <-> Independence
** Tree: 2
4,2 | 1,3 <-> Gumbel, parameters = 1.19309

fit3

4-dimensional vine copula fit ('vinecop')
nobs = 500   logLik = 183.85   npars = 5   AIC = -357.7   BIC = -336.63   

contour(fit3)

It appears that they will. The fit is not as good (as evidenced by the higher AIC and BIC), mostly because the Gumbel copula has only upper tail dependence (or, when rotated 180 degrees, lower); whereas, the t-copula has both upper and lower tail dependence. Therefore, we see a large discrepancy in the contour plots for \(C_{34}\) and \(C_{13}\), where there is not enough tail dependence, as compared to our input model. Curiously, \(C_{12}\) appears to have produced too much tail dependence.

Nonparametric copulas

Restrict the modeling to nonparametric copulas. Trying to fit a mixed set of parametric and nonparametric pair copulas causes only the parametric copulas to be chosen. It’s not clear why this is, since the nonparametric copulas produce a better fit to the data, and the nonparametric fits don’t seem to be getting penalized in AIC or BIC for the additional effective parameters.

fit4 <- vinecop(u, family_set=c('tll'), structure=vstruct, selcrit='mbicv', 
                psi0=0.5, show_trace=TRUE, cores=4)

** Tree: 0
4,3 <-> TLL, parameters = [30x30 grid]
3,1 <-> TLL, parameters = [30x30 grid]
2,1 <-> TLL, parameters = [30x30 grid]
** Tree: 1
4,1 | 3 <-> TLL, parameters = [30x30 grid]
3,2 | 1 <-> TLL, parameters = [30x30 grid]
** Tree: 2
4,2 | 1,3 <-> TLL, parameters = [30x30 grid]

fit4

4-dimensional vine copula fit ('vinecop')
nobs = 500   logLik = 271.73   npars = 0   AIC = -543.46   BIC = -543.46   

contour(fit4)

These look like pretty plausible approximations to our input model. The AIC and BIC are quite a bit lower than the best parametric fit, but the effective number of parameters is given as 0, which makes me suspect that the AIC and BIC are not being adjusted for the number of parameters (note that AIC = -2*logLik, which is another red flag). On the other hand, the fit doesn’t choose the nonparametric pair copulas when given a choice, so perhaps the adjustment is being applied internally and not reported in the summary output.

Automatic structure selection

The vinecop function can also attempt to select a graph structure automatically.

Elliptical pair copulas

fit5 <- vinecop(u, family_set=c('indep','elliptical'), selcrit='mbicv', 
                par_method='itau', psi0=0.5, show_trace=TRUE, cores=4)

** Tree: 0
3,1 <-> Student, parameters = 0.501381
 4.04004
2,1 <-> Gaussian, parameters = 0.494742
4,3 <-> Student, parameters = -0.253297
  2.09253
** Tree: 1
2,3 | 1 <-> Independence
4,1 | 3 <-> Gaussian, parameters = 0.129754
** Tree: 2
4,2 | 1,3 <-> Gaussian, parameters = 0.250494

fit5

4-dimensional vine copula fit ('vinecop')
nobs = 500   logLik = 223.46   npars = 7   AIC = -432.91   BIC = -403.41   

plot(fit5, tree='ALL')

contour(fit5)

The automatic selection draws the tree a bit differently, but it’s isomorphic to our original.

General pair copulas

fit6 <- vinecop(u, family_set=c('itau'), selcrit='mbicv', 
                par_method='itau', psi0=0.5, show_trace=TRUE, cores=4)

** Tree: 0
2,1 <-> Gaussian, parameters = 0.494742
3,1 <-> Student, parameters = 0.501381
 4.04004
4,3 <-> Student, parameters = -0.253297
  2.09253
** Tree: 1
3,2 | 1 <-> Independence
4,1 | 3 <-> Frank, parameters = 0.749713
** Tree: 2
4,2 | 1,3 <-> Gaussian, parameters = 0.252249

fit6

4-dimensional vine copula fit ('vinecop')
nobs = 500   logLik = 223.82   npars = 7   AIC = -433.63   BIC = -404.13   

plot(fit6, tree='ALL')

contour(fit6)

Again, essentially the same result as when we specified the structure.

Fitting a vine copula to a grouped-t model

The copula package has a mechanism for creating a hierarchically grouped meta-t model. That sort of model might provide a more interesting test for estimating structure.

This example is taken from pp. 95-97 of Elements of Copula Modeling with R.

library('copula')
set.seed(867-5309)

sim_grp_t <- function(N) {
  ## N is sample size
  dims <- seq(1,4)
  dimtot <- sum(dims)
  nu <- rep(4^seq(2,-1), dims)  # degrees of freedom for each grouping
  Z <- matrix(rnorm(N*dimtot), ncol=N)
  sig <- matrix(0.5, nrow=dimtot, ncol=dimtot)  # correlation matrix
  diag(sig) <- 1
  A <- t(chol(sig))
  Y <- t(A %*% Z)             # n x d matrix containing vectors distributed as N(0, sig)
  U <- runif(N)
  W <- sapply(nu, function(nx) {1/qgamma(U, shape=nx/2, rate=nx/2)})
  X <- sqrt(W) * Y            # These are the observations with t-distributed margins
  U <- sapply(seq(1,dimtot), function(j) {pt(X[,j], df=nu[j])})   # transform to uniform margins
  
  ## Return normal margins
  zobs <- qnorm(U)
  colnames(zobs) <- paste0('z',seq(1,dimtot))
  zobs
}

zgrpt <- sim_grp_t(1000)
ugrpt <- pseudo_obs(zgrpt)
ggpairs(as.data.frame(zgrpt),
        upper = list(continuous = wrap('cor', size=2.5)),
        lower = list(continuous = wrap('points', alpha=0.3, size=0.25)),
        progress=FALSE) + theme_bw()

Fit using all copula families with inverse-tau relations.

fit_grpt1 <- vinecop(ugrpt, family_set=c('itau'), selcrit='mbicv', 
                par_method='itau', psi0=0.5, show_trace=TRUE, cores=4)

** Tree: 0
4,1 <-> Student, parameters = 0.514564
 5.29693
4,2 <-> Student, parameters = 0.507129
  2.6696
4,3 <-> Student, parameters = 0.533786
 2.86814
6,4 <-> Student, parameters = 0.561683
 2.00001
7,5 <-> Student, parameters = 0.516288
 2.00001
7,6 <-> Student, parameters = 0.513101
 2.00001
9,7 <-> Student, parameters = 0.53187
2.00001
10,7 <-> Student, parameters = 0.525458
 2.00001
8,7 <-> Student, parameters = 0.525324
 2.00001
** Tree: 1
7,4 | 6 <-> Student, parameters = 0.290302
 2.00001
2,1 | 4 <-> Gumbel 180°, parameters = 1.28994
3,1 | 4 <-> Gaussian, parameters = 0.338643
6,2 | 4 <-> Student, parameters = 0.275081
 2.68418
8,5 | 7 <-> Student, parameters = 0.365936
 2.00001
8,6 | 7 <-> Student, parameters = 0.332642
 2.00001
8,9 | 7 <-> Student, parameters = 0.327435
 2.00001
8,10 | 7 <-> Student, parameters = 0.347717
 2.00001
** Tree: 2
3,2 | 1,4 <-> Student, parameters = 0.239712
 6.83347
2,7 | 4,6 <-> Joe 180°, parameters = 1.21154
8,4 | 6,7 <-> Student, parameters = 0.245893
 2.56118
9,6 | 7,8 <-> Student, parameters = 0.225882
  2.9277
6,1 | 2,4 <-> Gumbel 180°, parameters = 1.12553
10,6 | 7,8 <-> Student, parameters = 0.241739
 2.97688
6,5 | 7,8 <-> Gumbel 180°, parameters = 1.17609
** Tree: 3
9,4 | 6,7,8 <-> Student, parameters = 0.194241
 5.46649
6,3 | 1,2,4 <-> Student, parameters = 0.141228
 3.70006
8,2 | 4,6,7 <-> Gumbel, parameters = 1.12471
10,9 | 6,7,8 <-> Student, parameters = 0.175735
 2.00001
5,4 | 6,7,8 <-> Student, parameters = 0.230445
 4.92458
1,7 | 2,4,6 <-> Gumbel, parameters = 1.10015
** Tree: 4
5,9 | 4,6,7,8 <-> Student, parameters = 0.139466
 5.80893
10,4 | 6,7,8,9 <-> Student, parameters = 0.141851
  8.5237
5,2 | 4,6,7,8 <-> Gumbel 180°, parameters = 1.15983
7,3 | 1,2,4,6 <-> Student, parameters = 0.0880204
  7.42971
1,8 | 2,4,6,7 <-> Gaussian, parameters = 0.0856895
** Tree: 5
1,5 | 2,4,6,7,8 <-> Student, parameters = 0.10718
7.70394
2,9 | 4,5,6,7,8 <-> Independence
10,5 | 4,6,7,8,9 <-> Gumbel, parameters = 1.05912
8,3 | 1,2,4,6,7 <-> Gumbel, parameters = 1.05708
** Tree: 6
10,2 | 4,5,6,7,8,9 <-> Frank, parameters = 1.01307
9,1 | 2,4,5,6,7,8 <-> Independence
3,5 | 1,2,4,6,7,8 <-> Frank, parameters = 0.663804
** Tree: 7
1,10 | 2,4,5,6,7,8,9 <-> Independence
3,9 | 1,2,4,5,6,7,8 <-> Clayton, parameters = 0.150755
** Tree: 8
3,10 | 1,2,4,5,6,7,8,9 <-> Independence

fit_grpt1

10-dimensional vine copula fit ('vinecop')
nobs = 1000   logLik = 4524.08   npars = 68   AIC = -8912.15   BIC = -8578.43   

plot(fit_grpt1, tree=1:4)

contour(fit_grpt1)

The structure seems reasonably plausible, albeit with some odd choices. For example, 7-10 are connected to each other, but 5 hangs off of 7 for some reason, and isn’t connected to 4 and 6. Similarly, 2 and 3 are separated by 4.

The distributions in Tree 1 seem about right, but rvinecopulib has a lower limit of 2 df for its t distributions, so you can’t get the really heavy tail dependence that we see in z4-10, which were based on t-copulas with df = 1 (for 4-6) and df = 0.25 (for 7-10).

We can see how well we did by simulating from the vine distribution and making the pairs plot

usim1 <- rvinecop(Nsim, fit_grpt1)

zsim1 <- qnorm(usim1)
colnames(zsim1) <- paste0('z', seq(1,10))
ggpairs(as.data.frame(zsim1),
        upper = list(continuous = wrap('cor', size=2.5)),
        lower = list(continuous = wrap('points', alpha=0.3, size=0.25)),
        progress=FALSE) + theme_bw()

As expected, the distributions look pretty good for \(Z_1\)-\(Z_3\), but for \(Z_6\) - \(Z_{10}\) the tail dependence (i.e., the arms of the crosses) just isn’t there; however, this is probably due to the limitations on the allowable \(\nu\) values in rvinecopulib’s bivariate t.

Comparison to fit with prescribed structure

Same process as before, but we will force the tree to be hierarchical, with 1 at the root, with 2,4, and 7 (the group leaders) below that, and the members of each group below their group leaders.

mgrpt <- matrix(c(
    rep(7,3), 4,2,4, rep(1,4),
    rep(1,6), 2,2,2, 0,
    rep(2,4), 4,2,4,4,0,0,
    rep(4,3), rep(7,4), rep(0,3),
    rep(5,6), rep(0,4),
    rep(3,5), rep(0,5),
    rep(6,4), rep(0,6),
    rep(8,3), rep(0,7),
    rep(9,2), rep(0,8),
    10, rep(0,9)),
  nrow=10, ncol=10, byrow=TRUE)

vgrpt <- as_rvine_structure(mgrpt)
plot(vgrpt, tree = seq(1,4))

When fitting the pair copulas, we will also open up the copula families to allow some of the two-parameter Archimedean copulas.

fit_grpt2 <- vinecop(ugrpt, family_set=c('parametric'), structure=vgrpt, selcrit='mbicv', 
                par_method='mle', psi0=0.5, show_trace=TRUE, cores=4)

** Tree: 0
10,7 <-> Student, parameters = 0.811431
       2
9,7 <-> Student, parameters = 0.802158
       2
8,7 <-> Student, parameters = 0.80792
      2
6,4 <-> Student, parameters = 0.622843
       2
3,2 <-> Student, parameters = 0.504404
 3.92374
5,4 <-> Student, parameters = 0.563438
       2
7,1 <-> Student, parameters = 0.451049
 4.07094
4,1 <-> Student, parameters = 0.510615
 5.26557
2,1 <-> Student, parameters = 0.50138
5.77521
** Tree: 1
10,1 | 7 <-> Frank, parameters = 1.44937
9,1 | 7 <-> Gaussian, parameters = 0.207893
8,1 | 7 <-> Gaussian, parameters = 0.238788
6,1 | 4 <-> Gumbel 180°, parameters = 1.18613
3,1 | 2 <-> Gaussian, parameters = 0.348033
5,1 | 4 <-> Gumbel 180°, parameters = 1.22194
7,2 | 1 <-> Student, parameters = 0.293591
 2.36656
4,2 | 1 <-> Student, parameters = 0.336398
 3.16134
** Tree: 2
10,2 | 1,7 <-> Gumbel 180°, parameters = 1.16162
9,2 | 1,7 <-> Clayton, parameters = 0.200266
8,2 | 1,7 <-> Student, parameters = 0.206025
 10.1145
6,2 | 1,4 <-> Student, parameters = 0.206711
 3.41924
3,4 | 1,2 <-> Student, parameters = 0.297607
 4.83687
5,2 | 1,4 <-> Student, parameters = 0.246325
 4.62978
7,4 | 1,2 <-> Student, parameters = 0.389389
       2
** Tree: 3
10,4 | 1,2,7 <-> Student, parameters = 0.182316
 4.52844
9,4 | 1,2,7 <-> Student, parameters = 0.219031
 5.44472
8,4 | 1,2,7 <-> Student, parameters = 0.221756
 6.91201
6,7 | 1,2,4 <-> Student, parameters = 0.310783
 2.07637
3,7 | 1,2,4 <-> Student, parameters = 0.148369
 4.52282
5,7 | 1,2,4 <-> Student, parameters = 0.263032
 2.35066
** Tree: 4
10,5 | 1,2,4,7 <-> Student, parameters = 0.0758867
  6.91558
9,5 | 1,2,4,7 <-> Gumbel 180°, parameters = 1.09769
8,5 | 1,2,4,7 <-> Gumbel 180°, parameters = 1.13229
6,5 | 1,2,4,7 <-> Student, parameters = 0.128224
 5.28259
3,5 | 1,2,4,7 <-> Student, parameters = 0.133015
 8.83002
** Tree: 5
10,3 | 1,2,4,5,7 <-> Joe, parameters = 1.06305
9,3 | 1,2,4,5,7 <-> Gaussian, parameters = 0.11686
8,3 | 1,2,4,5,7 <-> Clayton 180°, parameters = 0.101146
6,3 | 1,2,4,5,7 <-> Gumbel, parameters = 1.07941
** Tree: 6
10,6 | 1,2,3,4,5,7 <-> Gumbel 180°, parameters = 1.07402
9,6 | 1,2,3,4,5,7 <-> Gumbel 180°, parameters = 1.08135
8,6 | 1,2,3,4,5,7 <-> Gumbel 180°, parameters = 1.08454
** Tree: 7
10,8 | 1,2,3,4,5,6,7 <-> Student, parameters = 0.271489
 3.51745
9,8 | 1,2,3,4,5,6,7 <-> Student, parameters = 0.287039
 2.80448
** Tree: 8
10,9 | 1,2,3,4,5,6,7,8 <-> Gumbel, parameters = 1.09925

fit_grpt2

10-dimensional vine copula fit ('vinecop')
nobs = 1000   logLik = 4401.97   npars = 72   AIC = -8659.94   BIC = -8306.58   

contour(fit_grpt2)

There is a big improvement in the AIC and BIC here. None of the two-parameter Archimedean copulas got used, so it’s some combination of the structure and the mle solver method. We’ll once again restrict to inverse-tau capable families so we can do a rigorous comparison between the two structures.

fit_grpt2a <- vinecop(ugrpt, family_set=c('itau'), structure=vgrpt, selcrit='mbicv', 
                par_method='itau', psi0=0.5, show_trace=TRUE, cores=4)

** Tree: 0
10,7 <-> Student, parameters = 0.525458
 2.00001
9,7 <-> Student, parameters = 0.53187
2.00001
8,7 <-> Student, parameters = 0.525324
 2.00001
6,4 <-> Student, parameters = 0.561683
 2.00001
3,2 <-> Student, parameters = 0.504807
 3.92658
5,4 <-> Student, parameters = 0.495106
 2.00001
7,1 <-> Student, parameters = 0.461495
 4.11173
4,1 <-> Student, parameters = 0.514564
 5.29693
2,1 <-> Student, parameters = 0.504932
 5.81607
** Tree: 1
10,1 | 7 <-> Frank, parameters = 1.5675
9,1 | 7 <-> Gaussian, parameters = 0.257549
8,1 | 7 <-> Gaussian, parameters = 0.29035
6,1 | 4 <-> Gumbel 180°, parameters = 1.18802
3,1 | 2 <-> Gaussian, parameters = 0.34456
5,1 | 4 <-> Gumbel 180°, parameters = 1.21019
7,2 | 1 <-> Student, parameters = 0.289525
 2.38584
4,2 | 1 <-> Student, parameters = 0.331277
 3.16882
** Tree: 2
10,2 | 1,7 <-> Gumbel 180°, parameters = 1.22113
9,2 | 1,7 <-> Gumbel 180°, parameters = 1.14679
8,2 | 1,7 <-> Student, parameters = 0.243179
 5.15208
6,2 | 1,4 <-> Student, parameters = 0.20873
3.18884
3,4 | 1,2 <-> Student, parameters = 0.284869
 4.79308
5,2 | 1,4 <-> Student, parameters = 0.241153
 4.30703
7,4 | 1,2 <-> Student, parameters = 0.282884
 2.00001
** Tree: 3
10,4 | 1,2,7 <-> Student, parameters = 0.173258
 2.31843
9,4 | 1,2,7 <-> Student, parameters = 0.21406
 2.5212
8,4 | 1,2,7 <-> Student, parameters = 0.221205
 2.73718
6,7 | 1,2,4 <-> Student, parameters = 0.247435
 2.00001
3,7 | 1,2,4 <-> Student, parameters = 0.133385
 4.29581
5,7 | 1,2,4 <-> Student, parameters = 0.221058
 2.14312
** Tree: 4
10,5 | 1,2,4,7 <-> Student, parameters = 0.0763741
  3.96541
9,5 | 1,2,4,7 <-> Student, parameters = 0.142361
 4.61162
8,5 | 1,2,4,7 <-> Joe 180°, parameters = 1.23567
6,5 | 1,2,4,7 <-> Student, parameters = 0.123024
 4.46591
3,5 | 1,2,4,7 <-> Gumbel 180°, parameters = 1.09828
** Tree: 5
10,3 | 1,2,4,5,7 <-> Clayton 180°, parameters = 0.0985896
9,3 | 1,2,4,5,7 <-> Clayton, parameters = 0.169023
8,3 | 1,2,4,5,7 <-> Joe, parameters = 1.08779
6,3 | 1,2,4,5,7 <-> Gumbel, parameters = 1.07391
** Tree: 6
10,6 | 1,2,3,4,5,7 <-> Student, parameters = 0.117117
 6.16926
9,6 | 1,2,3,4,5,7 <-> Joe 180°, parameters = 1.13494
8,6 | 1,2,3,4,5,7 <-> Student, parameters = 0.134332
 5.38902
** Tree: 7
10,8 | 1,2,3,4,5,6,7 <-> Student, parameters = 0.208102
 2.00001
9,8 | 1,2,3,4,5,6,7 <-> Student, parameters = 0.170346
 2.00001
** Tree: 8
10,9 | 1,2,3,4,5,6,7,8 <-> Student, parameters = 0.128771
 3.61083

fit_grpt2a

10-dimensional vine copula fit ('vinecop')
nobs = 1000   logLik = 4327.48   npars = 75   AIC = -8504.97   BIC = -8136.89   

contour(fit_grpt2a)

This one has AIC and BIC that are even better than the last, but the log-likelihood is actually a little lower. (How is this even possible, given that this model has more parameters than the previous one?)

Given that the inverse-tau version of the model is a little suspect, we will plot some samples from the MLE version. One thing of particular interest is that in tree 0 the MLE model increased the correlation parameter in the t-distributions, presumably to make up for not being able to get the degrees of freedom low enough to accurately represent the dats.

usim2 <- rvinecop(Nsim, fit_grpt2)

zsim2 <- qnorm(usim2)
colnames(zsim2) <- paste0('z', seq(1,10))
ggpairs(as.data.frame(zsim1),
        upper = list(continuous = wrap('cor', size=2.5)),
        lower = list(continuous = wrap('points', alpha=0.3, size=0.25)),
        progress=FALSE) + theme_bw()