Code - Note 5"

In gremes: Estimation of Tail Dependence in Graphical Models

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

Generate the data first

We do not have a real dataset according to a given block graph which is not a tree. Therefore we first generate the observations. How is the data generated is described in Vignette "Additional functionalities".

Create the block graph.

g<- graph(c(1,3,1,2,2,3,
            3,4,4,5,5,3,
            3,7,3,6,6,7), directed=FALSE)
g<- set.vertex.attribute(g, "name", V(g), c("a", "b", "c", "d", "e", "f", "g"))
plot(g)

Create the edge weights to be assigned to the edges of the graph.

# all deltas are squares already
C1<- c(0.2, 0.8, 0.6)     # d_13^2, d_12^2, d_23^2
C2<- c(0.3, 0.5, 0.1)     # d_34^2, d_45^2, d_35^2
C3<- c(0.4, 0.05, 0.25)   # d_37^2, d_36^2, d_67^2

Attach the edge weights to the edges.

hrmbgobj<- HRMBG(g)
hrmbgobj<- setParams(hrmbgobj, c(C1, C2, C3))
hrmbgobj

Create the matrix $\Lambda$, whose entry $\lambda_{ij}$ is the sum of the edge weights on the unique shortest path between node $i$ and node $j$.

hrmlam<- HRLambda(hrmbgobj)
hrmlam

Generate 1000 observations from Huesler-Reiss distribution with parameter matrix $\Lambda$ and some independent random noise.

X<- rHRM(hrmbgobj, hrmlam, 1000, noise = TRUE)

Perform estimation

Create the subsets for local estimation.

rdsobj<- RootDepSet()
rdsobj<- setRootDepSet(rdsobj, subset = list(c("a", "b", "d", "e", "f", "g"),
                                             c("a", "b", "d", "e", "f", "g"),
                                             c("a", "b", "d", "e", "f", "g"),
                                             c("a", "b", "d", "e", "f", "g"),
                                             c("a", "b", "d", "e", "f", "g"),
                                             c("a", "b", "d", "e", "f", "g")), 
                       c("a", "b", "d", "e", "f", "g"))

Estimate the model treating the third variable as latent: create first an object of class HRMBG and then use on it the method estimate.

hrmbg<- HRMBG(g)
hrmbg<- suppressMessages(estimate(hrmbg, X[,-3], rdsobj, k_ratio=0.2))
hrmbg$depParams

We have suppressed the messages. The messages are informative, if they don't stop the process they are not errors.

gremes documentation built on Feb. 16, 2023, 8:06 p.m.