Additional functionalities

The main scope of the package is the estimation of the parameters of tail dependence in the framework of a particular parametric model, namely using a family of Huesler-Reiss distributions. The package provides some few other functionalities which may be useful for post-estimation analysis or for simulations. Here we explain these functionalities.

library(gremes)
#> Loading required package: igraph
#> 
#> Attaching package: 'igraph'
#> The following objects are masked from 'package:stats':
#> 
#>     decompose, spectrum
#> The following object is masked from 'package:base':
#> 
#>     union
#> 
#> Attaching package: 'gremes'
#> The following object is masked from 'package:stats':
#> 
#>     sigma
#> The following object is masked from 'package:base':
#> 
#>     subset

Generate random sample from a Markov tree parameterized cliquewise by Huesler-Reiss distributions

Estimation in this package is based on different limit results of a random vector \(X=(X_v, v\in V)\) which satisfies the global Markov property with respect to a tree \(T=(V,E)\), and every two adjacent nodes have Huesler-Reiss distribution with some parameter \(\theta_e, e\in E\). The univariate margins of \(X\) are unit Frechet, \(F(x)=\exp(-1/x)\). For simulation purposes it is useful to be able to generate from the model for \(X\). The method that does this is rHRM.HRMtree.

seg<- make_tree(8,3, mode = "undirected") # create the tree
seg<- set.vertex.attribute(seg, "name", V(seg), letters[1:8]) # name the nodes
hrm<- HRMtree(seg) # initialize the object of of class HRMtree 
#> From HRMnetwork: Edges have been assigned names
hrm<- setParams(hrm, seq(0.1, 0.7, 0.1)) # set its parameters
#> From setParams.HRMtree: Names have been attributed to the vector 'value' in the order corresponding to the order of the edges: The fist element has the name of the first edge, the second element the name of the second edge, etc.
#> From setParams.HRMtree: The parameters have been attached to the edges according to their names
X<- rHRM(hrm, 1000) # generate a random sample 
round(head(X), 4)
#>            a       b       c       d       e       f       g      h
#> [1,]  1.3932  1.4150  1.0126  1.8260  1.4660  2.5956  1.9160 1.2930
#> [2,]  0.4412  0.4416  0.4772  0.4327  0.4597  0.4188  0.4232 1.5094
#> [3,] 19.9333 19.3102 14.6455 16.2861 26.0415 11.8430 12.9218 6.7473
#> [4,]  0.6681  0.5925  0.5896  0.6870  1.2189  0.7027  0.6339 0.4437
#> [5,]  2.8759  2.9952  2.3163  4.6884  2.8969  1.4057  2.3249 2.3464
#> [6,]  1.4125  1.5517  1.3106  1.1067  1.0615  1.3367  1.7118 2.9851
XX<- rHRM(hrm, 1000, noise = TRUE) # generate a random sample with independent normal noise
round(head(XX), 4)
#>           a      b      c      d       e      f      g      h
#> [1,] 1.4916 2.0960 0.9048 2.4574  1.5095 2.0003 2.5978 0.3779
#> [2,] 2.2289 3.2111 1.0032 1.0817  1.5838 1.8361 2.1998 1.3477
#> [3,] 2.5064 1.2841 0.7836 1.2402  3.1074 1.9519 1.9862 1.0273
#> [4,] 1.6602 2.3554 0.7867 1.8187  1.7332 1.0312 1.0765 1.2696
#> [5,] 2.1833 2.6961 2.2378 2.0931  1.0725 1.7687 2.9677 4.1910
#> [6,] 7.4321 5.8681 8.4692 5.9956 12.9982 6.4727 8.8752 2.3917
plot(seg)

For any \(u\in V\) the joint density function of \(X\) is \[\begin{equation} \label{eq:jdf} f(x)= f_{u}(x_u)\prod_{(v,j)\in E_u}f_{j\mid v}(x_j\mid x_v), \end{equation}\] with \(E_u \subseteq E\) the set of edges directed away from \(u\), i.e., \((v,j) \in E_u\) if and only if \(v = u\) or \(v\) separates \(u\) and \(j\). The joint density \(f\) is determined by \(d-1\) bivariate densities \(f_{vj}\). The univariate margins \(f_u\) are unit Frechet densities, \(f_j(x_j)=\exp(-1/x_j)/x_j^2\) for \(x_j \in (0, \infty)\), and the bivariate margins for each pair of variables on adjacent vertices \(j,v\) are Huesler–Reiss distributions with parameter \(\theta_{jv}\).

To generate an observation from the left hand-side of the equation above we use the right hand-side of that equation, proceeding iteratively, walking along paths starting from \(u\) using the conditional densities. An observation of \(X_j\) given \(X_v = x_v\) is generated via the inverse function of the cdf \(x_j \mapsto F_{j|v}(x_j \mid x_v)\), the conditional cdf of \(x_j\) given \(x_v\). To do so, the equation \(F_{j|v}(x_j \mid x_v)-p=0\) is solved numerically as a function in \(x_j\) for fixed \(p\in(0,1)\). The choice of the Huesler–Reiss bivariate distribution gives the following expression for \(F_{j\mid v}(x_j\mid x_v)\): \[\begin{equation*} \begin{split} % F_{l|k}(x_l|x_k)&= &\Phi\left( \frac{\theta_{jv}}{2}+\frac{1}{\theta_{jv}}\ln\frac{x_j}{x_v} \right) \cdot \exp\left[ -\frac{1}{x_v}\left\{\Phi\left(\frac{\theta_{jv}}{2}+\frac{1}{\theta_{jv}}\ln\frac{x_j}{x_v}\right)-1\right\} - \frac{1}{x_j}\Phi\left(\frac{\theta_{jv}}{2}+\frac{1}{\theta_{jv}}\ln\frac{x_v}{x_j}\right) \right]. \end{split} \end{equation*}\]

Diagnostic for the random sample generated as above

As a diagnostic tool we provide comparison between the bivariate true copula and the bivariate empirical copula. We fix a node, say \(a\) and if \(b\) and \(c\) are adjacent to \(a\) we compute the true bivariate copulas and the empirical copulas of pairs \((X_a, X_b)\) and \((X_a, X_c)\) and compare them. We also provide plot of the true and the empirical cumulative distribution function of the selected variable, e.g., of \(a\) in the preceding example.

diagnost(hrm, X, "b", y = c(0.3,0.5))

#> [1] "variable  b adj. variable  a ; true copula  0.299999999921372 ; empirical copula  0.3"
#> [1] "variable  b adj. variable  e ; true copula  0.29589026650931 ; empirical copula  0.296"
#> [1] "variable  b adj. variable  f ; true copula  0.291002999050035 ; empirical copula  0.294"
#> [1] "variable  b adj. variable  g ; true copula  0.28491676575979 ; empirical copula  0.284"

When we add some noise to the model the univariate empirical curve is a bit below the true one, but in the tail the difference diminishes.

Similarly for the tail of the bivariate copulas: for higher coordinates, (0.8, 0.9), the bivariate true and empirical copulas become closer to each other than for smaller values of the coordinates (0.3, 0.5):

diagnost(hrm, XX, "b", y = c(0.3, 0.5))

#> [1] "variable  b adj. variable  a ; true copula  0.299999999921372 ; empirical copula  0.261"
#> [1] "variable  b adj. variable  e ; true copula  0.29589026650931 ; empirical copula  0.264"
#> [1] "variable  b adj. variable  f ; true copula  0.291002999050035 ; empirical copula  0.264"
#> [1] "variable  b adj. variable  g ; true copula  0.28491676575979 ; empirical copula  0.241"
diagnost(hrm, XX, "b", y = c(0.8, 0.9))

#> [1] "variable  b adj. variable  a ; true copula  0.8 ; empirical copula  0.8"
#> [1] "variable  b adj. variable  e ; true copula  0.799431996325495 ; empirical copula  0.8"
#> [1] "variable  b adj. variable  f ; true copula  0.798248020703279 ; empirical copula  0.797"
#> [1] "variable  b adj. variable  g ; true copula  0.796402838196817 ; empirical copula  0.795"

Generate a random sample from max-stable Huesler-Reiss distribution with structured matrix

When \(X\) is Markov with respect to a block graph and parameterized cliquewise by a family with Huesler-Reiss distributions we do not dispose of a method to generate from the exact distribution of such a model as we do when the graph is a tree. We can generate from the max-stable attractor of this model. To do this we use the package mev Belzile et al. (2020).

bg<- graph(c(1,2,2,3,1,3,
             3,4,3,5,4,5,
             3,7,3,6,6,7), directed = FALSE) # create the graph
bg<- set.vertex.attribute(bg, "name", V(bg), letters[1:7]) # name the nodes
hrbg<- HRMBG(bg) # initialize an object with zero dependence parameters
#> From HRMnetwork: Edges have been assigned names
hrbg<- setParams(hrbg, seq(0.1, 0.9, 0.1)) # set the parameters
#> From setParams.HRMtree: Names have been attributed to the vector 'value' in the order corresponding to the order of the edges: The fist element has the name of the first edge, the second element the name of the second edge, etc.
#> From setParams.HRMtree: The parameters have been attached to the edges according to their names
lam<- HRLambda(hrbg) # compute the structured matrix Lambda, see Vignette "Huesler-Reiss distributions"
XB<- rHRM(hrbg, lam, 1000, noise = TRUE) 
round(head(XB), 4)
#>           a      b       c      d       e      f       g
#> [1,] 4.0174 3.4109  2.1142 1.4723  3.5631 0.6850  1.3459
#> [2,] 7.0732 5.6515  7.9537 1.3179  7.7321 7.2283 14.6200
#> [3,] 4.8353 5.9588 11.2259 4.2713 16.4408 7.4186  2.6119
#> [4,] 1.1597 4.5375  1.7774 1.1395  4.3121 2.2265  2.9631
#> [5,] 0.6790 0.9906  1.8080 2.9548  4.1811 2.6307  1.5099
#> [6,] 2.4132 1.7651  1.7561 1.3590  1.4822 1.7507  1.4831
plot(bg)

This code samples from a max-stable Huesler-Reiss copula with unit Frechet univariate margins and with parameter matrix \[ \lambda_{ij}^2=\sum_{e\in p(i,j)}\delta_e^2, \] where \(p(i,j)\) is the unique shortest path between two nodes \(i,j\). See Vignette “Huesler-Reiss distributions” also.

Of course the matrix Lambda can be one corresponding to a tree, given that tree is a special case of a block graph.

The stable tail dependence function (stdf)

The stdf is a key quantity in the topics of multivariate extremal dependence. It is defined as the following limit \[ l(x_v, v\in V)=\lim_{t\rightarrow\infty}t P\Big(\bigcup_{v\in V} \{X_v>t/x_v\}\Big), x\in (0,\infty)^{|V|} \] where \(X\) has unit Pareto univariate margins. The relation between the stdf and a max-stable (extreme value) distribution is given by \[ G(x_v, v\in V)=\exp\big\{-l(1/x_v, v\in V)\big\}, \] if \(G\) is an extreme value distribution with unit Frechet margins. More about the stdf can be read in Drees and Huang (1998), Beirlant et al. (2004), Haan and Ferreira (2007). The distribution \(G\) can have a parametric model such as the Huesler-Reiss model used throughout this package (see Vignette “Huesler-Reiss distributions” for their particular parameterizations).

The package provides tools for computing parametric and non-parametric stdf for models on trees and block graphs.

The parametric estimate of \(l\) when \(l\) is parameterized using the Huesler-Reiss distribution is \(l(x_v, v\in V; \hat{\theta}_{n,k})\) in case of tree models and \(l(x_v, v\in V; \hat{\delta}_{n,k})\) in case of models on block graphs. The vector \(\hat{\theta}_{n,k}=(\hat{\theta}_{e; n,k}, e\in E)\) is any estimate of the edge weights on a tree, and \(\hat{\delta}_{n,k}=(\hat{\delta}_{e; n,k}, e\in E)\) is any estimate of the edge weights on a block graph.

The non-parametric estimate of the stdf is given in Drees and Huang (1998) \[\begin{equation} \hat{l}_{n,k}(x_v, v\in V) = \frac{1}{k}\sum_{i=1}^n 1\left( \bigcup_{v \in V}\Big\{ n\hat{F}_{v,n}(X_{v,i}) >n+1/2-kx_v\Big\} \right), \end{equation}\] where \(\hat{F}_{v,n}(x)=(1/n)\sum_{i=1}^n1(X_{v,i}\leq x)\), i.e., the non-parametric estimate of the cumulative distribution function.

The methods which compute stdf are

Here we have some examples for these.

# non-parametric estimates on an object containing a tree 
x<- runif(8)
names(x)<- letters[1:8]
tobj<- Tree(seg, XX) # an object of containing block graph and the data associated to it
#> From validate.Network: Edges have been assigned names
#> From validate.Network: No latent variables
#> From validate.Network: Edges have been assigned names
#> From validate.Network: No latent variables
stdf(tobj, x, 0.2)
#> [1] 1.135
x<- x[3:8] # with latent variables on nodes "a" and "b"
names(x)<- letters[3:8] 
XU<- X[,3:8] 
tobjU<- Tree(seg, XU) # an object containing a tree an the data associated to it
#> From validate.Network: Edges have been assigned names
#> From validate.Network: There are nodes with latent variables
#> From validate.Network: Edges have been assigned names
#> From validate.Network: There are nodes with latent variables
stdf(tobjU, x, 0.25) 
#> [1] 1.112
x<- matrix(runif(40), 5,8)
colnames(x)<- letters[1:8]
stdf(tobj, x, 0.25)
#> [1] 1.072 1.060 1.296 1.132 1.168

# non-parametric estimates on an object containing a block graph
x<- runif(7)
names(x)<- letters[1:7]
bgobj<- BlockGraph(bg, XB) # an object containing a tree an the data associated to it
#> From validate.Network: Edges have been assigned names
#> From validate.Network: No latent variables
#> From validate.Network: Edges have been assigned names
#> From validate.Network: No latent variables
stdf(bgobj, x, 0.15)
#> [1] 1.566667

# parametric estimates on model with respect to a tree
x<- c(0, 0.1, 0, 2.5, 0, 1.3, 2.3, 1.5)
names(x)<- letters[1:8]
stdf(hrm, x )
#> From setRootDepSet.RootDepSet: The order of the subset must correspond to the
#>             order of its corresponding root
#>        b 
#> 3.253938 
#> attr(,"error")
#> [1] 2.349214e-43
#> attr(,"msg")
#> [1] "Normal Completion"

# parametric estimates on a model with respect to a block graph
x<- runif(7)
names(x)<- letters[1:7]
stdf(hrbg, x)
#> From setRootDepSet.RootDepSet: The order of the subset must correspond to the
#>             order of its corresponding root
#> From setRootDepSet.RootDepSet: The order of the subset must correspond to the
#>             order of its corresponding root
#> From setRootDepSet.RootDepSet: The order of the subset must correspond to the
#>             order of its corresponding root
#> From setRootDepSet.RootDepSet: The order of the subset must correspond to the
#>             order of its corresponding root
#> From setRootDepSet.RootDepSet: The order of the subset must correspond to the
#>             order of its corresponding root
#> From setRootDepSet.RootDepSet: The order of the subset must correspond to the
#>             order of its corresponding root
#> From setRootDepSet.RootDepSet: The order of the subset must correspond to the
#>             order of its corresponding root
#>        a 
#> 1.688213 
#> attr(,"error")
#> [1] 3.305471e-06
#> attr(,"msg")
#> [1] "Normal Completion"

Extremal coefficients

We look at extremal coefficients as some tools for post-estimation analysis. A typical analysis includes a comparison between estimates of extremal coefficients using the estimates of the edge weights and non-parametric extremal coefficients.

The extremal coefficient of variables in the set \(J\subseteq V\) is given by \[ l_J=l(1_J), \] where \(l\) is the stable tail dependence function and \(1_J=(1_{i\in J}, i\in V)\), i.e., a vector of length \(|V|\) whose elements are zero or one. An element \(i\) will be one if it belongs to \(J\) and zero otherwise. The range of a \(J\)-variate extremal coefficient is between 1 and \(|J|\) with dependence decreasing for large value of the coefficient.

There are three methods for computing extremal coefficients:

Non-parametric extremal coefficients - trees and block graphs

Given two objects of classes BlockGraph or Tree respectively we can compute the bivariate extremal coefficients as follows (k-ratio = 0.2).

extrCoeff(bgobj, 0.2) 
#>   a     b     c     d     e    f     g
#> a 0 1.205 1.360 1.495 1.505 1.61 1.545
#> b 0 0.000 1.295 1.520 1.500 1.59 1.560
#> c 0 0.000 0.000 1.460 1.430 1.56 1.515
#> d 0 0.000 0.000 0.000 1.495 1.63 1.605
#> e 0 0.000 0.000 0.000 0.000 1.60 1.610
#> f 0 0.000 0.000 0.000 0.000 0.00 1.560
#> g 0 0.000 0.000 0.000 0.000 0.00 0.000
extrCoeff(tobj, 0.2)  
#>   a     b     c     d     e     f     g     h
#> a 0 1.085 1.105 1.165 1.175 1.210 1.225 1.265
#> b 0 0.000 1.125 1.150 1.180 1.210 1.250 1.260
#> c 0 0.000 0.000 1.180 1.215 1.225 1.245 1.220
#> d 0 0.000 0.000 0.000 1.210 1.225 1.300 1.320
#> e 0 0.000 0.000 0.000 0.000 1.245 1.260 1.315
#> f 0 0.000 0.000 0.000 0.000 0.000 1.310 1.325
#> g 0 0.000 0.000 0.000 0.000 0.000 0.000 1.320
#> h 0 0.000 0.000 0.000 0.000 0.000 0.000 0.000

For an arbitrary dimension of the extremal coefficient, we need to create the vector of coordinates and name it according to the names of the nodes. If we want a four variate extremal coefficient of the variables \((X_a, X_c, X_d, X_e)\) we need to do

y<- c(1, 0, 1, 1, 1, 0, 0) 
names(y)<- letters[1:7]
extrCoeff(bgobj, 0.25, y) 
#> [1] 1.908
extrCoeff(tobj, 0.25, y)
#> [1] 1.376

If there are latent variables the bivariate extremal coefficients are computed between only observable variables.

# on the tree 
extrCoeff(tobjU, 0.2)
#>   c    d     e    f    g     h
#> c 0 1.13 1.200 1.21 1.23 1.225
#> d 0 0.00 1.195 1.18 1.23 1.280
#> e 0 0.00 0.000 1.24 1.22 1.315
#> f 0 0.00 0.000 0.00 1.26 1.330
#> g 0 0.00 0.000 0.00 0.00 1.340
#> h 0 0.00 0.000 0.00 0.00 0.000

# on the block graph
XBU<- XB[, -3] 
bgobjU<- BlockGraph(bg, XBU)
#> From validate.Network: Edges have been assigned names
#> From validate.Network: There are nodes with latent variables
#> From validate.Network: Edges have been assigned names
#> From validate.Network: There are nodes with latent variables
extrCoeff(bgobjU, 0.3)
#>   a        b        d        e        f        g
#> a 0 1.193333 1.446667 1.456667 1.530000 1.510000
#> b 0 0.000000 1.433333 1.426667 1.526667 1.490000
#> d 0 0.000000 0.000000 1.446667 1.576667 1.553333
#> e 0 0.000000 0.000000 0.000000 1.543333 1.533333
#> f 0 0.000000 0.000000 0.000000 0.000000 1.516667
#> g 0 0.000000 0.000000 0.000000 0.000000 0.000000

Note that for the parameters to be all identifiable on the tree, only variables on nodes \(a,b\) can be latent. On the block graph only variable \(c\) can be latent. This respects the requirement that all latent variables should be connected to at least three cliques Asenova, Mazo, and Segers (2021).

If we want an extremal coefficient of other dimension we need to pass a vector with non-zero coordinates only for observed variables.

v<- c(0,0,1,1,0,1,0,1)
names(v)<- letters[1:8]
extrCoeff(tobjU, 0.2, v)
#> [1] 1.455

v<- c(1, 1, 0, 1, 0, 0, 1)
names(v)<- letters[1:7]
extrCoeff(bgobjU, 0.15, v)
#> [1] 2.2

Parametric extremal coefficients

The method extrCoeff.HRMtree is used for models on trees.

extrCoeff(hrm) # bivariate
#>   a        b        c        d        e        f        g        h
#> a 0 1.039878 1.079656 1.119235 1.163330 1.201239 1.238977 1.284146
#> b 0 0.000000 1.089021 1.125633 1.158519 1.197413 1.235823 1.286697
#> c 0 0.000000 0.000000 1.143065 1.181231 1.215809 1.251150 1.273661
#> d 0 0.000000 0.000000 0.000000 1.201239 1.232620 1.265478 1.306198
#> e 0 0.000000 0.000000 0.000000 0.000000 1.251150 1.281568 1.324294
#> f 0 0.000000 0.000000 0.000000 0.000000 0.000000 1.303842 1.343254
#> g 0 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.364744
#> h 0 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
v<- c(0,0,1,1,0,1,0,1)
names(v)<- letters[1:8]
extrCoeff(hrm, v) # for a particular set of variables 
#> From setRootDepSet.RootDepSet: The order of the subset must correspond to the
#>             order of its corresponding root
#>        c 
#> 1.518242 
#> attr(,"error")
#> [1] 7.831616e-05
#> attr(,"msg")
#> [1] "Normal Completion"

The method extrCoeff.HRMBG is used for models on block graphs.

extrCoeff(hrbg) # bivariate
#>   a       b        c        d        e        f        g
#> a 0 1.24817 1.416118 1.597216 1.628907 1.705734 1.682689
#> b 0 0.00000 1.345279 1.561422 1.597216 1.682689 1.657218
#> c 0 0.00000 0.000000 1.472911 1.520500 1.628907 1.597216
#> d 0 0.00000 0.000000 0.000000 1.561422 1.726678 1.705734
#> e 0 0.00000 0.000000 0.000000 0.000000 1.745787 1.726678
#> f 0 0.00000 0.000000 0.000000 0.000000 0.000000 1.657218
#> g 0 0.00000 0.000000 0.000000 0.000000 0.000000 0.000000
v<- c(0,0,1,1,0,1,0)
names(v)<- letters[1:7]
extrCoeff(hrbg, v) # for a particular set of variables 
#> From setRootDepSet.RootDepSet: The order of the subset must correspond to the
#>             order of its corresponding root
#> From setRootDepSet.RootDepSet: The order of the subset must correspond to the
#>             order of its corresponding root
#> From setRootDepSet.RootDepSet: The order of the subset must correspond to the
#>             order of its corresponding root
#>        c 
#> 2.062956 
#> attr(,"error")
#> [1] 1e-15
#> attr(,"msg")
#> [1] "Normal Completion"

Tail dependence coefficients for models on trees

The tail dependence coefficient (tdc) of a subset \(W\subseteq V\) and a node \(u\notin W\) is defined as \[ \lim_{t\rightarrow\infty}\chi_{W|u}(t) := \lim_{t\rightarrow\infty}\Pr(X_W > t \mid X_u > t) \qquad u \not\in W, t>1. \] The event \(\{X_W > t\}\) is to be read as the intersection \(\bigcap_{w \in W}\{X_w > t\}\). When \(X\) has unit Pareto univariate margins \(\Pr(X_u > t) = 1/t\) for \(t > 1\), so that \[ \Pr(X_W > t \mid X_u > t)= t \, \Pr(X_{W \cup u} > t) \] If we set \(\overline{W}:=W\cup u\) we have \[\begin{equation} \label{eq:chiW} \chi_{\overline{W}} = \lim_{t \to \infty} \chi_{\overline{W}}(t) \end{equation}\] as a bivariate or multivariate tail dependence coefficient (tdc). In terms of the stable tail dependence function \(l\) of \(X_U\), the inclusion–exclusion formula yields \[ \chi_{\overline{W}} = \sum_{i=1}^{|\overline{W}|} (-1)^{i-1} \sum_{J \subseteq \overline{W}, |J|=i} l(1_J), \] where \(1_J = (1_{j \in J}, j \in U)\) is a vector of zeroes and ones. Note that \(l(1_J)\) is the extremal coefficient as considered above.

The parametric estimate of \(\chi_{\overline{W}}\) involves the parametric expressions of the stdf evaluated at a particular estimate of the parameter vector \(\theta=(\theta_e, e\in E)\).

Non-parametrically we estimate \(\chi_{\overline{W}}(t)\) at \(t = n/k\) via \[ \hat{\chi}_{\overline{W}} = \frac{n}{k} \frac{1}{n} \sum_{i=1}^n 1 \left\{ n\hat{F}_{v,n}(X_{v,i})>n-k, \, \forall v \in \overline{W} \right\}. \]

For the set \(\overline{W}=(a,c,d,e)\) we should do

v<- c(1,0,1,1,1,0,0,0)
names(v)<- letters[1:8]
suppressMessages(taildepCoeff(hrm, v)) # parametric tdc
#>         a 
#> 0.7420265 
#> attr(,"error")
#> [1] 0
#> attr(,"msg")
#> [1] "univariate: using pnorm"
taildepCoeff(tobj, 0.2, v) # non-parametric tdc
#> [1] 0.69

Confidence intervals for pairwise ECE - models on trees

Confidence intervals can be computed for the edge weights for models on trees if the pairwise extremal coefficients estimator is used. This is thanks to the distribution available of this estimator in Einmahl, Kiriliouk, and Segers (2017). The method that computes the confidence intervals is confInt.EKS.

Let \(\hat{\theta}_{n,k} = \hat{\theta}_{n,k}^{\mathrm{ECE}}\) denote the pairwise (bivariate) extremal coefficient estimator and let \(\theta_0\) denote the true vector of parameters. By Theorem 2 in Einmahl, Kiriliouk, and Segers (2017) with \(\Omega\) equal to the identity matrix, the ECE is asymptotically normal, \[ \sqrt{k}(\hat{\theta}_{n,k}-\theta_0)\sim \mathcal{N}_{|E|}\bigl(0, M(\theta_0)\bigr), \qquad n \to \infty, \] The asymptotic covariance matrix takes the form \[ M(\theta_0) = (\dot{L}^\top\dot{L})^{-1} \dot{L}^\top\Sigma_L\dot{L} (\dot{L}^\top\dot{L})^{-1}\, . \] The matrices \(\dot{L}\) and \(\Sigma_L\) depend on \(\theta_0\) and are based on partial derivatives of the stable tail dependence function.

For every \(k\) and every \(e \in E\), an asymptotic 95% confidence interval for the edge parameter \(\theta_{0,e}\) is given by \[ \theta_{0,e} \in \left[\hat{\theta}_{k,n;e}\pm 1.96\sqrt{\{M(\hat{\theta}_{k,n})\}_{ee}/k}\right] . \]

For more details on this interval we refer to section A.5 in Asenova, Mazo, and Segers (2021).

In the example below we suppose that the estimates obtained from the pairwise extremal coefficient estimator are given by the sequence \((0.1, 0.2, \ldots, 0.8)\) and the matrix of evaluation points based on pairs is the one based on all pairs. We suppose also that the estimates are based on \(k=150\).

# create the matrix of evaluation points 
tup<- Tuples()
x<- rep(1, 8)
names(x)<- letters[1:8]
pair<- evalPoints(tup, tobj, x)
# create an object of class EKS with the supposed estimates of the parameters
eks<- EKS(seg)
#> From HRMnetwork: Edges have been assigned names
eks<- setParams(eks, seq(0.1, 0.8, 0.1))
#> From setParams.HRMtree: Names have been attributed to the vector 'value' in the order corresponding to the order of the edges: The fist element has the name of the first edge, the second element the name of the second edge, etc.
#> From setParams.HRMtree: The parameters have been attached to the edges according to their names
suppressMessages(confInt(eks, pair, 150))
#>         [,1]      [,2]
#> e1 0.0432927 0.1567073
#> e2 0.1248752 0.2751248
#> e3 0.2090131 0.3909869
#> e4 0.2954711 0.5045289
#> e5 0.3832070 0.6167930
#> e6 0.4715271 0.7284729
#> e7 0.5579068 0.8420932

When there are latent variables we should use a matrix of evaluation points where all pairs are between observed variables only. Hence we should use this matrix that should have been used in estimation for computing the confidence intervals too.

x<- rep(1, 6)
names(x)<- letters[3:8]
pairU<- evalPoints(tup, tobjU, x )
suppressMessages(confInt(eks, pairU, 150))
#>           [,1]      [,2]
#> e1 -0.18534974 0.3853497
#> e2  0.08125466 0.3187453
#> e3  0.19472686 0.4052731
#> e4  0.28504080 0.5149592
#> e5  0.37882562 0.6211744
#> e6  0.46954141 0.7304586
#> e7  0.55860906 0.8413909

References

Asenova, Stefka, Gildas Mazo, and Johan Segers. 2021. “Inference on Extremal Dependence in the Domain of Attraction of a Structured Hüsler–Reiss Distribution Motivated by a Markov Tree with Latent Variables.” Extremes. https://doi.org/10.1007/s10687-021-00407-5.

Beirlant, J., Y. Goegebeur, J. Segers, J. L. Teugels, D. De Waal, and C. Ferro. 2004. Statistics of Extremes: Theory and Applications. Wiley Series in Probability and Statistics. Wiley. https://books.google.bg/books?id=GtIYLAlTcKEC.

Belzile, Leo, Jennifer L. Wadsworth, Paul J. Northrop, Scott D. Grimshaw, and Raphael Huser. 2020. Mev: Multivariate Extreme Value Distributions. https://CRAN.R-project.org/package=mev.

Drees, Holger, and Xin Huang. 1998. “Best Attainable Rates of Convergence for Estimators of the Stable Tail Dependence Function.” Journal of Multivariate Analysis 64 (1): 25–46. https://doi.org/https://doi.org/10.1006/jmva.1997.1708.

Einmahl, J., A. Kiriliouk, and J. Segers. 2017. “A Continuous Updating Weighted Least Squares Estimator of Tail Dependence in High Dimensions.” Extremes.

Haan, L. de, and A. Ferreira. 2007. Extreme Value Theory: An Introduction. Springer Series in Operations Research and Financial Engineering. Springer New York. https://books.google.bg/books?id=t6tfXnykazEC.