Estimating an Extreme Bayesian Network via Scalings

Introduction

Human society is continuously faced with challenges arising from factors of both uncontrollable and/or synthetic nature. The former is manifested through events such as natural disasters, in particular climate extremes like heavy rainfall or storms, unusually high/low temperatures, or river flooding. Similarly, synthetic factors correspond to those catastrophes influenced by human intervention, for instance industry fire, terrorist attacks, or a financial market crash. Such events occur rarely in isolation, but are rather interconnected, and occur simultaneously across certain instances; for example, floods disseminate through a river network, or extreme losses occur across several financial sectors. Such events make it necessary to not only understand dependencies between rare events, but also their causal structure.

When modeling rare events, one faces by definition a limited amount of data. While extremes in a univariate setting are well studied, multivariate extremes still remain a focus of present research. This is partly due to the augmented dimensionality problem, which affects crucially non-parametric methods (see de Haan and Ferreira (2006), Chapter 7), but also by the lack of a parametric family to characterize interdependencies (see Beirlant et al. (2004), Chapters 8, 9).

Recently, there has been interest in graphical models for modeling dependencies between extreme risks, which brings not only a potential complexity reduction, but also allows for modelling cause and effect in the context of extreme risk analysis. The model we consider in the present paper originates from Gissibl and Klüppelberg (2018), where max-linear structural equation models have been proposed and investigated. The underlying graphical structure of the model is a directed acyclic graph (DAG), also called a Bayesian network. Identifiability and estimation of recursive max-linar models are investigated in Gissibl et al. (2018). We refer to Lauritzen (1996) and Diestel (2010) for details on graphical modeling and graph theory, respectively.

Some other methods for combining graphical modeling with extremes have been proposed recently. In Segers (2019), Markov trees with regularly varying node variables are investigated using the so-called tail chains. In Engelke and Hitz (2018), a new approach using conditional independence relations between node variables is introduced, when considering undirected graphical models for extremes. This work is based on the assumption of a decomposable graph as well as the existence of density, which then leads to a Hammersley-Clifford type factorization of the latter into a lower dimensional setting. The method is applied to the estimation of flood in the Danube river network. A recursive max-linear model has been fitted to data from the EURO STOXX 50 Index in Einmahl et al. (2018), where the structure of the DAG is assumed to be known.

High dimensions are a serious challenge of dependence modeling of extreme events, and as a consequence most of the applications so far have focused on a lower dimensional setting. An exception is Cooley and Thibaud (2019), who present a new approach to extract the dependence structure from a regularly varying random vector. The authors propose the use of a dependence summary similar to the extreme dependence measure from Larsson and Resnick (2012), which can be considered an analogue to the covariance. Aiming at reducing the complexity, the authors propose a decomposition technique alike that of the Principal Component Decomposition for normal distributions. Other attempts aiming at dimension reduction of extremes involve Chautru (2015), and Janssen and Wan (2019), who present clustering approaches, or Haug et al. (2015), who propose a factor analysis for extremes.

In the present paper we develop a new structure learning and estimation algorithm for the recursive max-linear model in Gissibl and Klüppelberg (2018). Our approach is motivated by Cooley and Thibaud (2019) and applies to regularly varying node variables, which is a common assumption for extreme risk modelling. Multivariate node distributions have heavy-tailed marginals and are eponymous to those which lie in the domain of attraction of multivariate Fréchet distributions; see Resnick (1987), Section 5.4.2 (Proposition 5.15). We refer to Resnick (1987, 2007) for further details on regular variation.

Our multivariate regular variation setting is similar to that in Gissibl et al. (2018), which investigates the use of the tail dependence coefficients matrix towards the recovery of a causal order and identifiability of the max-linear coefficient matrix. Their method has the severe drawback that the initial nodes have to be known. For instance, in a DAG with two nodes and one edge there can only be one initial node, which can not be determined by the tail dependence coefficient as it is symmetric. This problem is to be encountered also in a DAG of larger size with several initial nodes, thus being a serious disadvantage to find a complete causal order. In contrast, other than regular variation itself, our methodology is free of assumptions. More recently, in a heavy-tailed setting, Gnecco et al. (2019) propose a method for identifying a causal order from the estimated conditional means of the integral transforms of pairs of nodes.

For arbitrary recursive max-linear models, a different identification and estimation method based on a generalized MLE can be found in Gissibl et al. (2018) and Klüppelberg and Lauritzen (2020). An extension of this method to models with observational noise has been investigated in Buck and Klüppelberg (2019). In these papers all innovations have to be independent and identically distributed, which is stronger than the tail assumptions imposed by regular variation.

We develop a new non-parametric methodology aimed at applying recursive max-linear models to extreme phenomena in a multivariate regular variation setting. First, targeting the problem of recovering a causal structure as a graphical model on a DAG, we propose a new technique via the scaling parameters of multivariate marginal distributions. This scaling technique allows for the manipulation of the dependence structure between extremes by simple scalar multiplication. These manipulations then uncover specific parts of the spectral measure, which fully characterize the dependence structure of interest to pave the way for estimating the causal dependence structure of the model. Second, we estimate the spectral measure empirically, where we focus on the relevant parts for the estimation of the required scaling. Asymptotic properties of the empirical spectral measure proven as an extension of a result of Larsson and Resnick (2012) lead to consistent and asymptotically normal estimates of all dependence parameters.

The application of the proposed methodology to financial data and to food dietary data shows that the recursive max-linear model can model multivariate extremes from real-life data with the goal of inferring causality for high risks.

Our paper is structured as follows. Preliminaries on graph theoretical terminology and regular variation, including the scaling parameter, are introduced in Section 2. Section 3 provides relevant properties of recursive max-linear models with regularly varying node variables. In Section 4 we show how the dependence structure of a recursive max-linear model can be identified from the scaling parameters of the model. Section 5 prepares for the causal inference by applying the scaling technique to find initial nodes as well as to reorder all other components into generations. Section 6 deals with statistical inference of the model. We propose non-parametric estimators of the relevant scalings, which also yield estimators of the dependence parameters. This allows us to estimate a partial order of the nodes and, in particular, a well-ordered graphical model on a DAG. We also show the asymptotic normality of the model dependence parameters. Finally, Section 7 is dedicated to two applications, namely to a real world financial data set of industry portfolio returns, as well as food dietary interview data.

Preliminaries

Let $\mathcal{D}=(V,E)$ be a directed acyclic graph (DAG) with nodes $V=\{1,\dots,{d}\}$ and edges $E=\{(j,i):i\in V\mbox{ and }j\in{\rm pa}(i)\}$, where ${\rm pa}(i)$ are the parents of node $i$. Each node of $\mathcal{D}$ is associated with a random variable, and dependence between two random variables can be represented via an edge connecting the corresponding nodes; for background see Lauritzen (1996).

Throughout we use the following notation. A path $p_{ji}\coloneqq[l_{0}=j\to l_{1}\to\cdots\to l_{m}=i]$ from node $j$ to $i$ has length $|p_{ji}|=m$, and we summarize all paths from $j$ to $i$ in the set $P_{ji}$.

For a node $i$ with parents ${\rm pa}(i)$ we set ${\rm Pa}{(i)}={\rm pa}(i)\cup\{i\}$, likewise, we denote by ${\rm an}(i)$ the ancestors of $i$ and set ${\rm An}(i)={\rm an}(i)\cup\{i\}$. The ancestral set of some subset $C\subset V$ of nodes is denoted by ${\rm an}(C)$ or ${\rm An}(C)={\rm an}(C)\cup C$. We also work with the following two notions throughout.

(i) We call $i\in V$ an initial node, if ${\rm pa}(i)=\emptyset$, and denote by $V_{0}$ the set of all initial nodes. (ii) In a DAG $\mathcal{D}$, a generation of nodes is the set of all nodes that have a longest path of same length from any initial node. Let $G_{0}=V_{0}$, then the $i$-th generation of nodes is defined by:

The following two auxiliary results provide some properties of this concept.

In a DAG $\mathcal{D}$ there is no path between two nodes of the same generation.

Suppose that there exists a path $p_{ij}$ in some generation $G_{k}\subset V$, $k\geq 1$ for nodes $i,j\in G_{k}$ on $\mathcal{D}$. A longest path from $V_{0}$ to $i$ would be of length $k$, say $p_{ti}$ for some $t\in V_{0}$. Extend now the same path along $p_{ij}$ to get $p_{tj}=[t\to\cdots\to i\to\cdots\to j]$. Clearly $p_{tj}$ is longer than $p_{ti}$, giving a contradiction to $j\in G_{k}$. ∎

The next result proves useful; its proof is not difficult and can be found in Krali (2018), Lemma 3.3.

Consider a DAG $\mathcal{D}=(V,E)$ with $|V|={d}$, and the set $V_{0}$ of initial nodes. Suppose that ${\mathcal{D}}$ has $l$ generations. Then for $i\in\{1,\dots,l\}$, $1\leq l\leq d$ and $k\notin G_{e}$ for $e<i$, we have $k\in G_{i}$ if and only if for all $j\in\underset{m\geq i}{\cup}G_{m}$ it holds that $j\notin\emph{an}(k)$.

A directed graph $\mathcal{D}=(V,E)$ is well-ordered, if for all $i\in V$ we have $i<j$ for all $j\in{\rm pa}(i)$. We refer to such an order as a causal order. $\Box$

Note that we employ a reverse ordering than in Gissibl and Klüppelberg (2018).

2 Multivariate Regular Variation

Multivariate regular variation can be defined in various ways, and we shall work with the following two equivalent definitions (cf. Resnick (2007), Theorem 6.1).

in $M_{+}((0,\infty]\times\Theta_{+}^{{{d}}-1})$, $d\nu_{\alpha}(x)=\alpha x^{-\alpha-1}dx$ for some $\alpha>0$, and for Borel subsets $C\subseteq\Theta_{+}^{{{d}}-1}$,

The measure $H_{\boldsymbol{X}}$ is called the spectral measure. (c) If $\boldsymbol{X}$ satisfies the above definition, we write $\boldsymbol{X}\in RV^{{d}}_{+}(\alpha)$, and $\alpha$ is called the index of regular variation. $\Box$

As explained in Theorem 6.5 of Resnick (2007), starting with an arbitrary vector $\boldsymbol{X}$ with positive components, we can always standardize all marginals to $\boldsymbol{X}\in RV^{{d}}_{+}(2)$ with normalizing sequence as in (b) chosen as $b_{n}=\sqrt{n}$. This implies that all scaling information is pushed into $H_{\boldsymbol{X}}$.

Let $\boldsymbol{X}\in RV^{{d}}_{+}(2)$ and consider its polar representation $(R,\boldsymbol{\omega})$ as in Definition 3(b) such that $\omega_{i}=\frac{X_{i}}{R}$ for $i=1,\dots,{d}$. For every $1\leq i,j\leq{{d}}$ define

We abbreviate $\sigma_{i}=\sigma_{{X}_{i}}=\sigma_{\boldsymbol{X}_{ii}}$ and call it the scaling or scaling parameter of $X_{i}$. $\Box$

The following auxiliary results are well-known and simple consequences of the definitions of regular variation. For the sake of completeness, we provide short proofs.

(a) From the homogeneity of the exponent measure and its polar representation in Definition 3(b) we obtain

(b) We simply compute the total mass of the ${d}$-dimensional unit simplex $\Theta_{+}^{{d}-1}$:

Immediately from Definition 4 and Lemma 3 above we find for $i\in\{1,\dots,{d}\}$ that, if $X_{i}$ has scaling $\sigma_{i}$, then $cX_{i}$ has scaling $c\sigma_{i}$ for every $c>0$.

(i) As $H_{\boldsymbol{X}}$ is a finite measure, it can be normalised to a probability measure by defining

For a simple assessment of the dependence structure of the components of a random vector, various summary measures have been introduced; see e.g. Sections 8.2.7 and 9.5.1 of Beirlant et al. (2004). We note that Definition 4 is a non-normalized version of the extreme dependence measure (EDM), which is a bivariate dependence measure on the positive unit sphere $\Theta_{+}^{{d}-1}$ and measures the limit of conditional cross-moments in the radial components of two random variables. The EDM has been introduced in Section 3 of Resnick (2004). A more refined version can be found in Propositions 3 and 4 in Larsson and Resnick (2012), where also more details on the EDM are given.

[Extreme dependence measure (EDM)] Let $\boldsymbol{X}\in RV^{d}_{+}(\alpha)$. Then for any two components $X_{i},X_{j}$ of $\boldsymbol{X}$, setting (\omega_{i},\omega_{j}):=\big{(}\frac{X_{i}}{R},\frac{X_{j}}{R}\big{)}, the EDM is given by

Recursive Max-linear Models

Recursive max-linear models were introduced in Gissibl and Klüppelberg (2018) and estimated with different methods in \al@gkl,GKO; \al@gkl,GKO; Klüppelberg and Lauritzen (2020). We summarize notation and results needed in the present paper.

A max-linear structural equation model $\boldsymbol{X}$ on a DAG $\mathcal{D}$ is defined as

From Theorem 2.2 of Gissibl and Klüppelberg (2018) we know that a max-linear structural equation model $\boldsymbol{X}$ from (3.1) has a solution in terms of its innovations $\boldsymbol{Z}$, which can be found by a path analysis. For each path $p_{ji}=[j\to k_{1}\to\cdots\to k_{l}=i]$ of length $l\geq 1$ from $j$ to $i$ define the path weights $d(p_{ji})\coloneqq c_{jj}c_{k_{1}j}\cdots c_{ik_{l-1}}$. Furthermore, define for $i=1,\dots,{d},$

Then $\boldsymbol{X}$ can be written as the recursive max-linear (ML) vector:

The matrix $A=(a_{ij})_{i,j=1,\dots,d}$ is called the ML coefficient matrix. Furthermore, a path $p_{ji}$ from $j$ to $i$ such that $a_{ij}=d(p_{ji})$ is called max-weighted.

If the innovations vector $\boldsymbol{Z}\in{RV}_{+}^{d}(\alpha)$, then by simple calculations given e.g. in Proposition A.2 of Gissibl et al. (2018), see also Proposition 4.1 of Krali (2018), $\boldsymbol{X}\in{RV}_{+}^{d}(\alpha)$ with discrete spectral measure

where $a_{k}=(a_{1k},\dots,a_{{d}k})^{\top}$ is the $k$-th column of $A$. Obviously, the entries of $A$ are the dependence parameters of $\boldsymbol{X}$.

Using the representation of $H_{\boldsymbol{X}}$ in (3.4) together with Remark 1 (ii) we obtain the following lemma.

For ${\alpha}=2$ and the Euclidean norm, the scalings of the recursive ML random vector (3.3) can be expressed by the matrix $A$ as follows.

From Definition 4 and (3.4) we find for $i\neq j$

The calculation of squared scalings $\sigma^{2}_{i}$ is analogous.

Finally, we consider the standardized recursive ML random vector $\boldsymbol{X}$ from (3.3) by standardizing the ML coefficient matrix.

[Standardized ML coefficient matrix] Define

Then $\bar{A}$ is referred to as standardized ML coefficient matrix. $\Box$

Since the innovations vector $\boldsymbol{Z}$ is standardized, all scaling information is in $\bar{A}$: Proposition 1 entails that the recursive ML vector $\boldsymbol{X}=\bar{A}\times_{\max}\boldsymbol{Z}$ has components with squared scalings $\sigma^{2}_{i}=(\bar{A}\bar{A}^{T})_{ii}=1$ for $i=1,\dots,{d}$. $\Box$

The following result has been proven in Lemma 2.1 of Gissibl et al. (2018).

Assume that the DAG corresponding to the recursive ML vector $\boldsymbol{X}$ is well-ordered. Then

We summarize all model assumptions used throughout the rest of the paper.

The innovations vector $\boldsymbol{Z}\in RV^{d}_{+}(2)$ has independent and standardized components.

We work with the Euclidean norm $\|\cdot\|_{2}$.

The ML coefficient matrix $A$ is standardized as in eq. (3.5), such that the components of $\boldsymbol{X}$ are standardized.

Identification of the ML Coefficient Matrix From Scalings

In this section we consider a recursive ML vector $\boldsymbol{X}=A\times_{\max}\boldsymbol{Z}$ such that (A1)-(A3) are satisfied. We show how to identify $A$ from $\boldsymbol{X}$, when $\boldsymbol{X}$ is a recursive ML vector on a well-ordered DAG; i.e.,

We identify $A$ from the squared scalings of maxima over combinations of components of $\boldsymbol{X}$. For a set $\boldsymbol{h}\subseteq\{1,\dots,{d}\}$ we define

We first compute the relevant squared scalings.

The random variable $M_{\boldsymbol{h}}$ is again max-linear, in particular $M_{\boldsymbol{h}}\in RV^{1}_{+}(2)$ with squared scalings as follows: (a) Let $\boldsymbol{h}\subseteq\{1,\dots,{d}\}$, then

(b) If $\boldsymbol{h}=\{1,\dots,{d}\}$, then $\sigma_{M_{\boldsymbol{h}}}^{2}=\sum_{k=1}^{d}a_{kk}^{2}.$

(a) Starting with $X_{i}=\underset{j=1,\dots,{d}}{\bigvee}a_{ij}Z_{j}$ for $i=1,\dots,{d}$, we calculate:

By eq. (3.4) $M_{\boldsymbol{h}}$ is regularly varying. In order to compute the squared scaling $\sigma_{M_{\boldsymbol{h}}}^{2}$, we use the same arguments as in Proposition 1 which yields (4.3). (b) This follows directly from (a) in combination with Lemma 5. ∎

We now illustrate the identification of the ML coefficient matrix by the following example.

Let $\boldsymbol{X}$ be a recursive ML vector on a well-ordered DAG satisfying (A1)-(A3) such that

Note that by standardization every row must have norm 1.

We first compute the diagonal entries. By standardization, $a^{2}_{33}=\sigma_{3}^{2}=1$. By Lemma 5 we know that $a_{ii}>a_{ki}$ for $k<i$. Let $M_{ij}$ for $1\leq i,j\leq 3$ and $M_{123}$ be defined as in (4.2). From Lemma 6 we obtain

From this we first find $a_{11}^{2}=\sigma_{M_{123}}^{2}-\sigma_{M_{23}}^{2}$. Similarly $a_{22}^{2}=\sigma_{M_{23}}^{2}-\sigma_{3}^{2}$.

The next step is to find the remaining entries in the first row of $A$, namely $a_{12}$ and $a_{13}$. Proceeding with $a_{12}$ we find from Lemma 6 for $M_{13}$ first $\sigma_{M_{13}}^{2}=a_{11}^{2}+a_{12}^{2}+a_{33}^{2}=a_{11}^{2}+a_{12}^{2}+\sigma_{3}^{2},$ which yields

Finally, we find $a_{13},a_{23}$, since the rows of $A$ have norm 1. $\Box$

We now proceed by proving the correctness of the above recursion, which gives rise to Algorithm 1 below.

Let $\boldsymbol{X}$ be a recursive ML vector on a well-ordered DAG satisfying (A1)-(A3). Then the following recursion yields the standardized ML coefficient matrix $A$:

We first show (4.5). From Lemma 6 we find for $i=1,\dots,{d}-1$:

which implies that $a_{ii}^{2}=\sigma_{M_{i,\dots,{d}}}^{2}-\sigma_{M_{i+1,\dots,{d}}}^{2}$. For $i=p$, by standardization of $A$ we have $a_{{d}{d}}^{2}={\sigma_{{d}}^{2}}=1$. In order to prove (4.6) we compute first:

Fix now $i\in\{1,\dots,{d}-1\}$. We proceed by induction over $j$. We start with the initial index $j=i+1.$ By (4.5) we know all $a_{ii}$ for $i=1,\dots,{d}$, and by (4.8),

By the induction hypothesis, suppose that we have found $a_{ij}$ for all $j\in\{i+1,...,l-1\}$, where $i+2<l<{d}$. Let $j=l$. Then, it is straightforward to see that

Equation (4.7) follows from the fact that $A$ is standardized, hence, all rows have norm 1 (by Remark 2, $\sigma_{i}^{2}=\sum_{k=i}^{d}a_{ik}^{2}=1.$) ∎

The Algorithm corresponding to Proposition 2 reads as follows.

In Proposition 2 we have shown that we can compute the diagonal entries of $A$ from the squared scalings $\sigma_{M_{1,2,\dots,{d}}}^{2},\sigma_{M_{2,3,\dots,{d}}}^{2},\dots,$ $\sigma_{M_{{d}-1,{d}}}^{2},\sigma_{d}^{2}$ by a recursion algorithm. Furthermore, we have identified the non-diagonal entries of the $i$-th row of $A$ from

Consider the row-wise vectorized version of the squared entries of the upper triangular matrix $A$, where we use $A^{2}$ for the matrix with squared entries of $A$ and its vectorized version

Note that both vectors ${A^{2}}$ and ${S}_{M}$ show a similar structure, built from row vectors with ${d},{d}-1,\ldots,1$ components, respectively; so both have ${d}({d}+1)/2$ components. By means of Proposition 2 we show that ${A^{2}}$ can be written as a linear transformation of $S_{M}$.

Let $S_{M}$ and ${A^{2}}$ be as in (4.9) and (4.10), respectively. Then

$a^{2}_{ii}:\quad t_{l_{ii},l_{ii}}=1$, $t_{l_{ii},l_{i+1,i+1}}=-1$ for $i=1,\dots,{d}-1$;

$a^{2}_{{d}{d}}:\quad t_{l_{ii},l_{ii}}=1$ for $i={d}$;

$a^{2}_{ij}:\quad t_{l_{ij},l_{ij}}=1,t_{l_{ij},l_{j+1,j+1}}=-1,t_{l_{ij},l_{i,j-1}}=-1,t_{l_{ij},l_{jj}}=1$ for $i<j\leq{d}-1$;

$a^{2}_{i{d}}:\quad t_{l_{i{d}},l_{i{d}}}=1,t_{l_{i{d}},l_{i,{d}-1}}=-1,t_{l_{i{d}},l_{{d}{d}}}=1$ for $i=1,\dots,{d}-1$,

where $l_{ij}=(j-{d})+\sum_{k=0}^{i-1}({d}-k)$ for $i=1,...,d$ and $j\geq i$. All other entries of $T$ are equal to zero.

(ii) Starting from (4.6) we show by induction that for $i=1,\dots,{d}-2$ and $j=i+1,\dots,{d}-1$,

For $j=i+1$ we clearly have that $a_{i,i+1}^{2}=(\sigma_{M_{i,i+2,\dots,{d}}}^{2}-\sigma_{M_{i+2,\dots,{d}}}^{2})-(\sigma_{M_{i,i+1,\dots,{d}}}^{2}-\sigma_{M_{i+1,\dots,{d}}}^{2}).$ Suppose now that this holds for all $i<j<{d}-2$. We show now that it holds for $j+1$. More specifically,

where the last equality is due to the telescoping sum after noting that $\sigma_{M_{i,i,\dots,{d}}}^{2}=\sigma_{M_{i,\dots,{d}}}^{2}$.

(iii) Similar to (ii), for $a_{i{d}}$ with $i<{d}$,

while for $i={d}$ we obtain again $a_{{d}{d}}^{2}=\sigma_{{d}}^{2}$.

(iv) The results in (i)-(iii) show already the linearity between the vectors ${A^{2}}$ and $S_{M}$. It remains to construct the matrix $T=(t_{uv})_{k\times k}$ for $k={d}({d}+1)/2$ such that ${A^{2}}=TS_{M}$. We start by renumbering the vector components in ${A^{2}}$ and replacing the double indices $ij$ for $i=1,\dots,{d}$ and $j\geq i$ by

Then the vector in (4.10) becomes $(a^{2}_{1},a^{2}_{2},\dots,a^{2}_{{d}({d}+1)/2})$. Moreover, (4.14) maps $ii$ into $l_{ii}$ and $l_{i,i+k}=l_{ii}+k$ for $i=1,\dots,{d}$ and $1\leq i+k\leq{d}$.

Also notice that for all $i=1,\dots,{d}$, by the structure of $S_{M}$, its $l_{ij}$-th component is $S_{l_{ij}}=\sigma_{M_{i,j+1,\dots,{d}}}^{2}$ for $i\leq j<{d}$, and $S_{l_{i{d}}}=\sigma_{i}^{2}$.

(v) We construct now $T$, where by (i)-(iii) $T$ contains many zeros, and we focus on the non-zero entries.

Since $a^{2}_{ii}$ becomes $a^{2}_{l_{ii}}$ for $i=1,\dots,{d}$ and, by the structure of $S_{M}$, the $l_{11},\dots,l_{{d}{d}}$-th components of $S_{M}$ are $\sigma_{M_{1,2,\dots,{d}}}^{2},\sigma_{M_{2,3,\dots,{d}}}^{2},$ $\dots,\sigma_{M_{{d}-1,{d}}}^{2},\sigma_{M_{d}}^{2}=\sigma_{d}^{2}$, respectively, in each $l_{ii}$-th row of $T$ there must be a 1 on the diagonal; i.e. $t_{l_{ii},l_{ii}}=1$. Furthermore, $t_{l_{ii},l_{i+1,i+1}}=-1$, and all other entries in this row are 0.

Similarly, we find the other non-zero entries by representation (4.12) and (4.13). ∎

We illustrate the linear transformation (4.11) for ${d}=4$, which clarifies the structure also for higher dimensions. For a recursive ML vector with 4 nodes, by (4.12) and (4.13) the identity ${A}^{2}=TS_{M}$ becomes

Reordering the Vector Components

In Section 4 we have assumed that the DAG underlying the recursive ML vector $\boldsymbol{X}$ is well-ordered. In a real life situation this will rarely be the case, and the components of $\boldsymbol{X}$ have to be reordered. In this section we use again the scalings for finding a causal order of the components of $\boldsymbol{X}$. This is achieved by first identifying the initial nodes, which can be ordered arbitrarily within all initial nodes. The same applies for every following generation: within each generation the order is arbitrary. All such obtained partial orders correspond to equivalent well-ordered DAGs and we construct one representative DAG by the method as follows.

We start with an auxiliary result which ensures that the recursive ML vector $\boldsymbol{X}=A\times_{\max}\boldsymbol{Z}$ is invariant with respect to column permutations of the ML coefficient matrix $A$.

Denote by $\pi:\{1,\dots,{d}\}\to\{1,\dots,{d}\}$ an arbitrary permutation of the columns of $A$, and notice that an arbitrary component $i\in\{1,\dots,{d}\}$ of $\boldsymbol{X}$ is given by

and, therefore, $\boldsymbol{X}^{\pi}=\boldsymbol{X}$. ∎

Since by Lemma 7 the distribution of $\boldsymbol{X}$ is invariant with respect to column permutations, we can assume that an arbitrarily ordered recursive ML vector $\boldsymbol{X^{*}}=(X_{1^{*}},\dots,X_{{d}^{*}})=A^{*}\times_{\max}\boldsymbol{Z}$, needs only row permutations in $A$, denoted by $\nu:(1^{*},\dots,{d}^{*})\to(1,\dots,{d})$, to become well-ordered:

We refer to entries of the matrix $A^{*}$ as $a_{i^{*}k}$ and to entries from the row-permuted matrix $A_{\nu}$ as $a_{ik}\coloneqq a_{\nu(i^{*})k}$, corresponding to a reordered vector $\boldsymbol{X}_{\nu}=\boldsymbol{X}$ in distribution on a well-ordered DAG.

In order to find an initial node, we fix one node which we want to investigate and extend the notation from (4.2) to maxima over a ${d}$-tuple of partly scaled random variables: for $a>0$ we define for $m\in\{1,\dots,{d}\}$,

By Lemma 6, also $M_{{-m},am}\in RV_{+}(2)$. The following theorem provides necessary and sufficient conditions for the identification of initial nodes.

Let $\boldsymbol{X^{*}}=(X_{1^{*}},\dots,X_{d^{*}})$ be an arbitrarily ordered recursive ML vector with ML coefficient matrix $A^{*}$ satisfying (A1)-(A3). Then the following holds. (a) If ${m^{*}\in\{1^{*},\dots,{d}^{*}\}}$ is an initial node of the recursive ML vector $\boldsymbol{X}_{\nu}$ in a well ordered DAG, then for all scalars $a>1$ it holds that

(b) If there exists a scalar $a>1$, such that for ${m^{*}\in\{1^{*},\dots,{d}^{*}\}}$ eq. (5.2) holds, then $m^{*}$ is an initial node of the recursive ML vector $\boldsymbol{X}_{\nu}$ in a well ordered DAG.

$(a)$ Let $X_{m^{*}}$ be the component of $\boldsymbol{X^{*}}$ such that $m^{*}$ is an initial node. W.l.o.g. we may set $X_{\nu(m^{*})}=X_{{d}}$. By the representation (4.1), and given the standardized scalings, we know that $a_{\nu(m^{*}),1}=\dots=a_{\nu(m^{*}),{d}-1}=0$, and $a_{\nu(m^{*}),{d}}=1$. By Lemma 6(b), for some $a>1$, we compute

Taking the difference yields (5.2). $(b)$ We prove this by contradiction. Let $\nu$ be a row permutation that transforms $\boldsymbol{X}^{*}$ into a recursive ML vector $\boldsymbol{X}_{\nu}$ on a well-ordered DAG, and suppose that for some non-initial node, say $k^{*}\in\{1^{*},\dots,{d}^{*}\}$, there exists some $a>1$ such that

Given that a recursive ML vector can have more than one initial node, w.l.o.g. assume that there are $0<l<{d}$ initial nodes. Since $k^{*}$ is not an initial node, we know by Definition 2 that $\nu({k^{*}})\leq{d}-l$. Furthermore, as ${k^{*}}$ must have an ancestor, there exists some node $j>\nu({k^{*}})$ for $j\in\{1,...,{d}\}\setminus{\{\nu(k^{*})\}}$, such that $j\in\textrm{an}(\nu({k^{*}}))$ and, hence, $a_{\nu({k^{*}})j}>0$. Since $\nu$ is a row permutation, which permutes $\boldsymbol{X^{*}}$ into a recursive ML vector on a well ordered DAG, we know that there exists some $j^{*}\in\{1^{*},...,{d}^{*}\}\setminus{\{k^{*}\}}$, such that $\nu(j^{*})=j$. By Lemma 5, and since $j>\nu({k^{*}})$, and $j\in\textrm{an}(\nu({k^{*}}))$, it follows that $a_{jj}>a_{\nu(k^{*})j}>0$. This implies that

Next, for the summands in the sum on the right-hand side, following Lemma 5, we obtain

which, when combined with eq. (5.4), yields the inequality

However, this is a contradiction to eq. (5.3). ∎

2 Reordering the Vector Components: Finding the Descendants

Once we have identified the initial nodes of the recursive ML vector $\boldsymbol{X}$, we provide a necessary and sufficient criterion for identifying the causal order of the descendants.

We proceed iteratively by identifying every new generation in the DAG. Suppose we have found all nodes which belong to a certain number of generations, and that there are $h\leq{d}-1$ such nodes which we have ordered as ${d},{d}-1,\dots,{d}-h+1$. Let $X_{\nu^{-1}({d})},\dots,X_{\nu^{-1}({d}-h+1)}$ be the corresponding components in the arbitrarily ordered recursive ML vector $\boldsymbol{X^{*}}.$

The next logical step is to investigate, whether ${m^{*}}\in\{1^{*},\dots,{d}^{*}\}\setminus\{\nu^{-1}({d}),\dots,$ $\nu^{-1}({d}-h+1)\}$, belongs to the next generation of nodes in the causal order. Define $\boldsymbol{h}\coloneqq\{\nu^{-1}({d}),\dots,\nu^{-1}({d}-h+1)\}$ and let $\boldsymbol{h}^{c}$ contain all other components. Then we take the maximum over a ${d}$-tuple of partly scaled random variables: for $a>0$ define

Let $\boldsymbol{X^{*}}=(X_{1^{*}},\dots,X_{{d}^{*}})$ be an arbitrarily ordered recursive ML vector with ML coefficient matrix $A^{*}$ satisfying (A1)-(A3). Let $\boldsymbol{h}=\{\nu^{-1}({d}),\dots,\nu^{-1}({d}-h+1)\}$ be the first $h$ nodes of the recursive ML vector $\boldsymbol{X}_{\nu}$ of a well-ordered DAG, which have already been ordered. Then the following holds: (a) If ${m^{*}}\notin\boldsymbol{h}$ has no ancestors in $\boldsymbol{h}^{c}$, then for all scalars $a>1$ it holds that

(b) If there exists a scalar $a>1$ such that (5.7) holds, then we identify ${m^{*}}\notin\boldsymbol{h}$ as the $(h+1)$-th node.

$(a)$ W.l.o.g. let ${m^{*}}$ be a node such that $\nu(m^{*})={d}-h$. Consider the squared scaling of ${M}_{\boldsymbol{h}_{a},m^{*}_{a},\{\boldsymbol{h}\cup\{m^{*}\}\}^{c}}$ as in (5.6), and $M_{1,\dots,{d}}$. By representation (4.1) and Lemma 6, and following similar steps as in the proof of Theorem 2(i), we find

Therefore, (5.7) is satisfied. $(b)$ Suppose now that for $X_{m^{*}}$ (${m^{*}}\notin\boldsymbol{h}$) there exists some $a>1$ such that

We have the following system of equalities:

The summands in the last summation are non-negative. For $\nu(m^{*})={d}-h$ we can take $m^{*}$ as the $(h+1)$-th node. Then the difference is equal to $(a^{2}-1)\sigma_{M_{\boldsymbol{h},m^{*}}}^{2}$. Similarly, if $a_{\nu(m^{*})j}=0$ for $\nu(m^{*})+1\leq j\leq{d}-h$, then this implies that ${\nu^{-1}(j)}\notin\textrm{an}(m^{*})$ and, thus, that $m^{*}$ can be chosen as the $(h+1)$-th node. Suppose now that $\nu(m^{*})<{d}-h$, and $a_{\nu(m^{*})j}>0$ for some $\nu(m^{*})+1\leq j\leq{d}-h$. Then, by Lemma 5 and since $a>1$, a bound similar to that in (5.5) gives

Therefore we have that either $\nu(m^{*})={d}-h$, or $a_{\nu(m^{*})j}=0$ for $\nu(m^{*})+1\leq j\leq{d}-h$. In both cases $m^{*}$ can be chosen as the $(h+1)$-th node. ∎

One of the consequences of the proof of Theorem 3 provides a criterion, when two or more components of $\boldsymbol{X}$ are neither descendants nor ancestors of one another in a DAG.

Let $\boldsymbol{X}$ be as in Theorem 3 and suppose that we have found the first $h$ nodes. If $\sigma_{\boldsymbol{M}_{\boldsymbol{h}_{a},m^{*}_{a},\{\boldsymbol{h}\cup\{m^{*}\}\}^{c}}}^{2}-\sigma_{M_{1,\dots,{d}}}^{2}=(a^{2}-1)\sigma_{M_{\boldsymbol{h},m^{*}}}^{2}$ for $m^{*}\in\{i^{*},j^{*}\}\cap\boldsymbol{h}^{c}$ and some $a>1$, then $a_{i^{*}j^{*}}=a_{j^{*}i^{*}}=0.$

As another direct consequence of Theorem 3 we obtain the following corollary.

Let $\boldsymbol{X}$ be as in Theorem 3 and suppose that we have found the first $h$ nodes. Let $X_{i^{*}},X_{j^{*}}$ be such that $X_{i^{*}}$ is the $(h+1)$-th node, while $X_{j^{*}}$ belongs to a different generation. Define for $m^{*}\in\{1^{*},\dots,{d}^{*}\}$ and $a>1$

[The reordering algorithms for a 10 nodes model] We assess the performance of the reordering procedure based on Theorems 2 and 3. We assume that the recursive ML vectors have standard Fréchet(2) components, which allows us to estimate the scalings by the standard MLE given for $M_{\boldsymbol{h}}$ by (see Krali (2018), Section 5.4 for details)

To perform the reordering we turn both Theorem 2 and Theorem 3 into Algorithm 2 and Algorithm 3, respectively. For every node, say ${i^{*}}$, Algorithm 2 checks the criterion in Theorem 2 (see Line 4 below). The difference between $\sigma_{M_{-\nu(i^{*}),a\nu(i^{*})}}^{2}-\sigma_{M_{1,\dots,{d}}}^{2}$ and $a^{2}-1$ is entered into the ${i^{*}}$-th- component of the vector $\Delta$. Finally, the vector $N$ is filled with non-zero components only for those nodes, which satisfy the criterion of Theorem 2. Algorithm 3 follows a similar logic.

When estimating the scalings, then the $\hat{\Delta}_{i^{*}}$-s in Algorithms 2 and 3 are a.s. different from 0. Hence, both algorithms are adapted to allow for small bounds to both quantities, which have to be chosen appropriately. We are interested in checking if the final order identifies the nodes in accordance with their respective generations. To this end, we choose a small bound $\epsilon_{3}$, which enables Algorithm 3 to return more than one node per iteration step; namely the generations. Based on simulation experience, we chose $a=\sqrt{2}$ and $\varepsilon_{1}=0.1,\varepsilon_{2}=0.05,\varepsilon_{3}=0.1$. All simulations and data analysis are done using R, R Core Team (2016).

As the initial order is irrelevant, we set w.l.o.g. $\boldsymbol{X}^{*}_{10}=(X_{1},\dots,X_{10})$. We perform 100 simulation runs for each of the sample sizes $n\in\{2000,3000,5000,$ $10000\}$.

The reorderings are obtained by first applying Algorithm 2 and then Algorithm 3. For some simulations it has happened that the bounds in Line 5 in Algorithms 2 and 3 are not satisfied by any of the components. We indicate this in Table 5.1, where the column “Valid Runs” corresponds to the number of simulation runs (out of $100$) for which the conditions of the bounds in the algorithms are satisfied each time the “if” loop is entered.

The column “Correctly Reordered” gives the number of runs for which Algorithms 2 and 3 return the generations $V_{0}=\{10\},G_{1}=\{8,9\},G_{2}=\{5,6,7\},G_{3}=\{1,2,3,4\}$. The column “Success Ratio” presents the ratio of “Correctly Reordered” over “Valid Runs”.

Statistical Theory for Regularly Varying Innovations

The discrete spectral measure of the ML model poses serious challenges towards the objective of estimation. Einmahl et al. (2016), Einmahl et al. (2012), and Einmahl et al. (2018) develop estimation procedures for the stable tail dependence function. In both, Einmahl et al. (2012) and Einmahl et al. (2018), the methods are also applied to models with discrete spectral measure, whose dependence parameters can be obtained from the stable tail dependence function. Janssen and Wan (2019) provide a new way for estimating the atoms of the spectral measure on the unit sphere by using a clustering approach. For our purposes we resort to the empirical spectral measure.

Finally, we estimate $A$ by the linear transformation as given in Theorem 1.

The following CLT holds for every $\boldsymbol{X}\in RV^{d}_{+}(2)$. Again we use the polar representation (6.1). Let the radial component $R$ of $\boldsymbol{X}$ have distribution function $F$.

2 Asymptotic Normality of the Scalings of Maxima

We show asymptotic multivariate normality of the vector $\hat{S}_{M}$. To this end we use the Cramér-Wold device and a properly chosen continuous function $f$ on $\Theta_{+}^{{d}-1}$ to which we then apply Theorem 4.

where the entries of the covariance matrix $W_{M}$ are given by the right-hand sides of the following two limits. The diagonal entries for $\boldsymbol{h}\subset\{1,\dots,{d}\}$ satisfy

and the non-diagonal entries for two different sets $\boldsymbol{h}_{i}\neq\boldsymbol{h}_{j}\subset\{1,\dots,{d}\}$ are given by

Moreover, the covariance matrix $W_{M}$ is singular.

which is—as a linear function of continuous functions—itself continuous on $\Theta_{+}^{{d}-1}$. The empirical estimator for $f$ is by (6.3) given as

Applying Theorem 4 for the given choice of $f$ it follows that

To show that $W_{M}$ is singular, let $t_{i}=1$ for $\boldsymbol{h}_{i}$ such that $|\boldsymbol{h}_{i}|=1$, and set the remaining components of the vector $\boldsymbol{t}$ to zero. Summarize all these $\boldsymbol{h}_{i}$ into the set $\mathcal{H}^{1}\coloneqq\{i\in\{1,\dots,{d}({d}+1)/2\}:|\boldsymbol{h}_{i}|=1\}$, and note that $|\mathcal{H}^{1}|={d}$ since there are exactly ${d}$ such entries in $S_{M}$, namely $\sigma_{1}^{2},\sigma_{2}^{2},...,\sigma_{d}^{2}$ corresponding to the dimension of $\boldsymbol{X}$. Then we obtain

Next, when computing the covariance we simply use the identity $2\textrm{Cov}(X,Y)=\textrm{Var}(X+Y)-\textrm{Var}(X)-\textrm{Var}(Y)$. Let $\boldsymbol{h}_{i}\neq\boldsymbol{h}_{j}\subseteq\{1,\dots,{d}\}$. Consider f(\boldsymbol{\omega})={d}\big{(}\underset{k\in\boldsymbol{h}_{i}}{\bigvee}\omega_{k}^{2}+\underset{l\in\boldsymbol{h}_{j}}{\bigvee}\omega_{l}^{2}\big{)}. Then using again Theorem 4 we get

The asymptotic covariance matrix $W_{m}$ can also be expressed in terms of the squared entries of the matrix $A$.

Let the assumptions of Theorem 5 hold. Then the entries of $W_{M}$ are given by the right-hand sides of the following two limits: On the diagonal we obtain

where $(\sum_{j=1}^{d}(\underset{i\in\boldsymbol{h}}{\bigvee}a_{ij}^{2}))^{2}=\sigma_{M_{\boldsymbol{h}}}^{4}$. For the non-diagonal entries we obtain

where $\sum_{k=1}^{d}\underset{m\in\boldsymbol{h}_{i}}{\bigvee}a_{mk}^{2}=\sigma_{M_{\boldsymbol{h}_{i}}}^{2}$ and $\sum_{k=1}^{d}\underset{l\in\boldsymbol{h}_{j}}{\bigvee}a_{lk}^{2}=\sigma_{M_{\boldsymbol{h}_{j}}}^{2}.$

With the explicit form of the spectral measure (3.4), expression (6.7) becomes:

Similarly, for the covariance, from (6.8) we obtain:

Summing the variance terms together we obtain

The identities giving $\sigma_{M_{\boldsymbol{h}}}^{2}$, $\sigma_{M_{\boldsymbol{h}_{i}}}^{2}$, and $\sigma_{M_{\boldsymbol{h}_{j}}}^{2}$ follow from (4.3). ∎

Finally we prove asymptotic normality of the estimated ML coefficient matrix $A$ computed via Theorem 1 as $A^{2}=TS_{M}$. As $T$ is a deterministic matrix, we obtain ${\hat{A}}^{2}$ as $T\hat{S}_{M}$. Then Theorem 5 and Corollary 3 gives the asymptotic normality of the estimated ML coefficient matrix $A$.

Let $\boldsymbol{X}=A\times_{\max}\boldsymbol{Z}$ be a recursive ML vector satisfying (A1)-(A3), and let $\boldsymbol{X}_{1},\dots,\boldsymbol{X}_{n}$ be i.i.d. copies of $\boldsymbol{X}$. Let the assumptions of Theorem 5 hold and assume that the ML coefficient matrix $A=({a}_{ij})_{{d}\times{d}}$ satisfies

Estimate ${\hat{A}}^{2}=T\hat{S}_{M}$ with $T$ as in Theorem 1 and $\hat{S}_{M}$ as in (6.5). Then

where $W_{M}$ is the covariance matrix in Corollary 3.

Data Applications

For estimating a causal order as well as for the estimation of the ML coefficient matrix $A$ we need estimates for the scalings in Algorithms 1, 2, and 3, respectively. Notice that in Proposition 2 we estimate the ML coefficient matrix $A$ for a well-ordered DAG, so that we first estimate the order of the nodes and then $A$. The structure learning is based on the scalings of $M_{\boldsymbol{h}_{a},m^{*}_{a},\{\boldsymbol{h}\cup\{m^{*}\}\}^{c}}$ as defined in (5.6), and the estimation of $A$ as in Proposition 2 is based on scalings of $M_{\boldsymbol{h}}$ as in (4.2). In contrast to Example 3 we make no distributional assumptions on $\boldsymbol{X}$, but use the non-parametric estimation method developed in Section 6.

For the non-parametric estimation of all scalings needed in the algorithms, as in Section 6 we denote by $k$ the number of upper order statistics corresponding to the radial threshold used for the estimation of the spectral measure. We recall that it has to be chosen as $k=o(n)$ and we choose $k\approx\sqrt{n}$. We observe that we may have rather few components exceeding the radii computed as in (6.1) based on all $d$ components. If some components of an observation are very large, then other components may not exceed the corresponding threshold. As the choice of the $k$ upper order statistics is based on these radii, there may be rather few exceedances in some component. Hence, we resort to lower dimensional vectors $\boldsymbol{X_{{\boldsymbol{q}}}}=(X_{i}:i\in{\boldsymbol{q}})$ for appropriate sets ${\boldsymbol{q}}\subseteq\{1,...,d\}$ for the estimation of the various scalings.

We want to apply Algorithms 2 and 3 for structure learning by replacing the theoretical scalings by their estimated counterparts. However, we have to modify both algorithms to account for estimation errors.

As the limited number of exceedances of radii in some components is particularly critical for identification of the initial nodes, we modify the estimation procedure, which has been presented in (5.1), and estimate for every $m\in\{1,\dots,d\}$ the initial nodes only based on $\max(X_{i},aX_{m})$ for $a>1$ (corresponding to $\boldsymbol{q}=\{i,m\}$), but then for all $i\in\{1,\dots,d\}\setminus\{m\}$. This means that the identification of the initial nodes is carried out by applying a pairwise version of Theorem 2 for all pairs. The following pairwise version of Algorithm 2 identifies the intial nodes of the DAG. It is based on Theorem 5.6 and Algorithm 3 of Krali (2018), adapted for possible estimation errors. The positive bounds $\varepsilon_{1},\varepsilon_{2}$ have to be chosen appropriately to ensure that the estimates $\hat{\Delta}_{i^{*}m^{*}}$ are close to zero. The scalar $a>1$ has to be chosen in accordance with Theorem 2.

Once the initial nodes are identified, we proceed finding the descendants. When searching for the $(h+1)$-th node, we apply the following modification of Algorithm 3, which has again been adapted for possible estimation errors by an application of Corollary 2.

In contrast to Example 3, where our goal was to identify the generations of the graph, here we are only interested in a causal order of the nodes, and thus modify Algorithm 3 based on Corollary 2 so that it returns a unique node at each step of the if-loop.

1.2 Estimation of the Scalings

According to Line 4 of Algorithm 5 three squared scalings need to be estimated:

-${\sigma}_{{M}_{\boldsymbol{h},i^{*},\{\boldsymbol{h}\cup\{i^{*}\}\}^{c}}}^{2}$ : By (5.6), this estimate involves all $d$ components of $\boldsymbol{X}$. Thus, we set $\boldsymbol{q}\coloneqq\{1,...,d\}$ and proceed as in step (i) below;

-${\sigma}_{M_{\boldsymbol{h},i^{*}}}^{2}\hskip 40.40285pt$: We estimate the spectral measure based on $\boldsymbol{q}\coloneqq\boldsymbol{h}\cup\{i^{*}\}$ and proceed as in step (i) below; such scalings we also need to estimate in Line 5 of Algorithm 4;

-${\sigma}_{{M}_{\boldsymbol{h}_{a},i^{*}_{a},\{\boldsymbol{h}\cup\{i^{*}\}\}^{c}}}^{2}$: By (5.6), this is the estimated scaling of the rescaled vector $\boldsymbol{X}$ and we follow step (ii) below.

For Algorithm 1 we have to estimate the following squared scalings for $i,j\in\{1,...,d\}$ and $i\leq j+1$:

-${\sigma}_{M_{i,j,j+1,...,d}}^{2}$: We estimate the spectral measure based on $\boldsymbol{q}\coloneqq\{i,j,j+1,...,d\}$.

-${\sigma}_{M_{j,j+1,...,d}}^{2}$ : Here we set $\boldsymbol{q}\coloneqq\{j,j+1,...,d\}$.

The following two steps modify the setting of Section 6 and summarize the estimation of the scalings in both Algorithms 4, 5, and Algorithm 1.

Then for $1\leq k\leq n$ we estimate ${\sigma}_{M_{\boldsymbol{q}}}^{2}$ as

When $\boldsymbol{q}=\{i\}$, corresponding to the scaling of a single component, then by (7.1), $\omega=1$, and plugging this in the estimator in (2), we obtain $\hat{\sigma}_{i}^{2}=1$, which is the true scaling parameter $\sigma_{i}=1$.

The numerator, $(a^{2}-1)(|\boldsymbol{h}|+1)+d$ corresponds to the new mass of the spectral measure as a consequence of the scaling by $a$ of the components involved in $X_{\boldsymbol{h}}$ and $X_{m^{*}}$, see for instance Lemma 3(b).

1.3 Estimating the ML Coefficient Matrix

After having estimated also the scalings needed for Algorithm 1, we have to take care of estimation errors. Indeed, it can happen that entries of $A^{2}$ are estimated as being negative. For the two data examples to follow we simply set $\hat{A}=\sqrt{\max(\hat{A}^{2},0)}$, with the square root taken entrywise, and keeping all positive estimates. For larger networks it may be advisable to choose a thresholding or lasso procedure to obtain a sparse graph.

2 Industry Portfolio Data

The data consist of seven time series of value-averaged daily percentage returns, each assigned to one of seven industry portfolios as part of the 30-Industry-Portfolio in the Kenneth French Data Library available at https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/datalibrary.html.

All 30 portfolios have been analysed in Cooley and Thibaud (2019); Janssen and Wan (2019), where in Cooley and Thibaud (2019) it is suggested that the tails of the data are regularly varying. We refer to the introduction for more details on the objectives of these papers.

The seven portfolios we consider are: Chemicals (Ch), Fabricated Products (FP), Electrical Equipment (EE), Healthcare (H), Smoke (S), Utilities (U), and Others (O), which includes products which are not specific to any of the other listed industries. A precise description of the data, in particular of each industry sector can be found on the website above.

The data has been collected over the years 1950-2015. Since the time series over this time period is non-stationary and, in particular, since also the dependence structure changes over this long period, we have selected the time window from $01.06.1989$ to $15.06.1998$ containing $2285$ observations, which show marginal stationarity. To each of the 7 time series we have fitted moving average processes of order 3, with the exception of Others where we have fitted a moving average process of order 4, and performed a Ljung-Box test with 8 lags on the residuals. The test did not reject the independence hypothesis of the residuals, supporting the assumption that the time series are stationary. Figure 2 depicts the time series plots in %-returns for the seven industries.

As we assess dependence in extreme negative returns, we transform the vector of the 7 portfolio time series to $\boldsymbol{X}^{*}=\max(-\boldsymbol{X},0)$. Our aim is to fit a recursive ML model model to $\boldsymbol{X}^{*}$. We transform the data by the empirical integral transform to standardize them to Fréchet$(2)$ margins (see for instance, p. 381 in Beirlant et al. (2004), or Cooley and Thibaud (2019)). We map $(S,H,EE,U,O,FP,Ch)\mapsto(1,2,3,4,5,6,7)$ and define for $i=1,\dots,7$

Running maxima. In order to provide some insight in the data structure, we start our analysis with a time-line of the high risk events of the seven industry portfolios and their association with shocks entering the industry network from the relevant component of the innovation vector $\boldsymbol{Z}$. The horizontal axes of Figure 3 shows every time point, when the maximum over all seven standardized industry returns happens and indicates on the vertical axes the respective innovation component causing this maximum. Since the innovations $Z_{i}$ for $i=1,\dots,d$ are atomfree and by Assumption (1) independent, representation (3.2) ensures that each recursive ML component $X_{i}$ realises its maximum in exactly one innovation. If two components of $\boldsymbol{X}$ are realised by the same innovation, still the realised values of the two components are different by Lemma 5, which implies that $\max\{X_{1},...,X_{d}\}$ is unique. These innovations are indicated in Figure 3.

We find that most of the shocks originate from the initial node Chemicals (7) during the time period 1990-1991 which coincides with the Gulf War, caused by the invasion of Kuwait by Iraq. During the war it was feared that Iraq made use of chemical warfare. Regarding Fabricated Products (6), the larger losses occur close to the end of 1997, which is associated with the slowdown in the Asian economies and which had spillover effects on the U.S economy, eventually leading also to the October 27, 1997 Mini-Crash. The Utilities (4) experience large losses in the years 1994 and 1996 associated with deregulation of the electric energy supply in the US, which was initiated in 1992. The Healthcare sector (2) experiences losses in the period 1992-1993 associated with the Clinton Health Care reform.

Our final goal is to approximate the causal dependence structure via a recursive ML model by means of the learning algorithm presented in the previous sections of this paper. This algorithm is based on all returns above a high threshold, not only on the maximum value.

Bivariate extremes. In a second exploratory analysis we plot the bivariate extremes (real data and simulated ones) in Figure 6 of Appendix A. We also simulate a 7-dimensional random vector $\boldsymbol{X}\in RV^{7}_{+}(2)$ from an innovation vector $\boldsymbol{Z}\in RV^{{7}}_{+}(2)$ with independent standard Fréchet(2) components of dimension $n=2285$ via $\boldsymbol{X}=\hat{A}\times_{\max}\boldsymbol{Z}$, where the estimated matrix $\hat{A}$ is given in (7.3). Two columns always belong together, the left one gives the empirical bivariate extremes of two of the seven components, respectively, which have also been the basis for the estimation procedure of Section 7.1. The right one presents a simulation of the estimated model. We plot only those bivariate observations with the 50 largest radii. Left and right (real and simulated data) look very much alike, indicating that the estimated bivariate models are valid approximations to the biviariate empirical distribution in the tails.

We give an interpretation of Utilities versus Electrical Equipment (line 6, columns 3 and 4): Here we notice that large losses in Utilities do not necessarily occur together with high risk for Electrical Equipment. This can be verified in the plot of realised bivariate extremes (column 3), where large losses for Utilities occurring close to the vertical axis correspond only to negligible losses for Electrical Equipment. This suggests that the common large losses between Utilities and Electrical Equipment are rather caused by their common ancestors. This may be due to idiosyncratic risks associated with one but not the other, which in this case might correspond to the innovation terms $Z_{3}$ $(\hat{a}_{43}=0)$ and $Z_{4}$ $(\hat{a}_{34}=0)$.

Fitting a recursive ML model. Finally, we approximate the extreme dependence structure of $\boldsymbol{X}^{*}$ by a recursive ML model. To this end, according to Section 6, we have to choose a threshold value $k=o(n)$ of the radial components. We choose $k\approx\sqrt{n}$ and set $k=50$. For identification of the initial nodes we employ Algorithm 4 with $a=1.01$ and $\varepsilon_{1}=0.0045,\varepsilon_{2}=0.0045$, and then in order to reorder the remaining nodes Algorithm 5. As output we obtain the ordered vector

Finally, we estimate the ML coefficient matrix by the estimation version of Algorithm 1 and setting $\hat{A}=\sqrt{\max(\hat{A^{2}},0)}$. The estimated standardised ML coefficient matrix is given by

The estimated squared scaling parameters of the components of $\boldsymbol{X}$ are obtained by summing the squared entries of the respective row. We find the estimated vector of scalings $(1.036,1.01,1.01,1,1,1,1)$ and recall that the theoretical ones are all equal to 1. Deviations from scalings of 1 stem from the fact that we set $\hat{A}=\sqrt{\max(\hat{A}^{2},0)}$. In doing so, once $\hat{A}^{2}$ has been computed, we ignore its entries that are close to zero but negative, for instance $\hat{a}_{12}^{2},\hat{a}_{14}^{2},\hat{a}_{24}^{2},\hat{a}_{34}^{2}$. Consequently this can make the sums of the square entries of the respective row be slightly greater than one.

The DAG corresponding to $\hat{A}$ is given in Figure 4. We recall that $a_{ij}=0$ implies no edge from $j$ to $i$.

The estimated DAG should provide insight into the causality structure of the seven industries. We associate the estimated scalings with the standardised risk of each component. Then the estimated entries in $\hat{A}$ indicate the proportions of risk inferred from the causal dependence.

The only initial node is Chemicals, whose high losses impact risk on all other industries as it has out-degree 6. For instance, the line of productions in Fabricated Products, Healthcare (which includes pharmaceutical industry), Electrical Equipment, Utilities (which includes electricity services and supply), and Smoke are highly dependent on the supply of chemical products and chemical processing.

Fabricated Products has out-degree 5 with its high risk affecting Utilities, Others, Smoke, Healthcare, and Electrical Equipment. A reason for this may lie in the industrial fabrication of many products in the affected industries. In particular the impact on Smoke, whose production line depends heavily on machinery, is stronger than that on Utilities, Electrical Equipment, and Others, since $a_{16}>\max(a_{26},\dots,a_{56})$.

Others, whose components include also Cogeneration Power Producers, has out-degree 4 with its high risk affecting Utilities, Smoke and Electrical Equipment, and Healthcare to a lesser extent. On the other hand, Others has in-degree 2, so high risk in Chemicals or Fabricated Products affects Others, which can be seen from $a_{57}$ and $a_{56}$ with higher influence from Chemicals.

Electrical Equipment has out-degree 2 and in-degree 3, so its high risk is caused by Others, Chemicals, and Fabricated Products, where the influence of Chemicals is about twice as large as Others, and the influence of high risk in Fabricated Products is much lower. On the other hand, high risk in Electrical Equipment impacts on Healthcare and Smoke in about equal proportions.

Utilities, Healthcare, and Smoke have out-degree 0, so these portfolios are affected by high risk of other portfolios, but their high risks do not spread elsewhere. This is seen from columns 1,2, and 4 of $\hat{A}$, where the quantities on the diagonal correspond to the idiosyncratic risk. The quantities to the right measure the high risk influencing these three portfolios.

3 Dietary Supplement Data

The data is taken from a dietary interview from the NHANES report for the year 2015-2016, which is available at https://wwwn.cdc.gov/Nchs/Nhanes/2015-2016/DR1TOT_I.XPT; here also more details about the 168 data components can be found. The objective is that of estimating the total intake of calories, nutrients and non-nutrient food components from foods and beverages consumed a day prior to the interview. From the above data 38 components have been investigated in Janssen and Wan (2019) using the clustering approach mentioned already in the introduction.

We focus on four of the components, Vitamin A (DR1TVARA), Beta-Carotene (DR1TBCAR), Lutein+Zeaxanthin (DR1TLZ) and Alpha-Carotene (DR1TACAR). We abreviate them as VA, BC, LZ, AC, respectively. For each component there are $n=9544$ observations, each corresponding to a different individual and generated from survey interviews, thus the data sample can be treated as an i.i.d. sample. A Hill plot (see e.g. Embrechts et al. (1997), Section 6.4) suggests that all data components are regularly varying with some positive index. We map (VA, BC, LZ, AC) $\mapsto(1,2,3,4)$ and apply the empirical integral transform to standardize the data to Fréchet(2) margins as in (7.2). As in Section 7.2 we choose $k\approx\sqrt{n}$ as radial threshold (see Section 6), taking $k=100$ upper order statistics. The bivariate extremes (real data and simulated ones) are plotted in Figure 7 of Appendix B with interpretations as in Section 7.2.

We apply Algorithms 4 and 5 to reorder the nodes and estimate a recursive ML model. In the two algorithms we set $a=1.01,\varepsilon_{1}=0.002,\varepsilon_{2}=0.001$. This results in the causal order (VA, BC, LZ, AC). Following the procedure in Section 7.2, we obtain the ML coefficient matrix

The estimated scaling parameters of the components of $\boldsymbol{X}$ equals up to three digits (1,1,1,1). The DAG corresponding to $\hat{A}$ is presented in Figure 5. We interpret dependence in high amounts of the four given food components. The estimated entries in $\hat{A}$ indicate the proportion of high intake from food consumption.

From the DAG in Figure 5 we observe that the only initial node is Alpha-Carotene. Having out-degree 3 its high intake affects the intake of the other three food components. Alpha-Carotene affects in particular Vitamin A, and Beta-Carotene, where in both cases it behaves as the main contributor of high values amongst the respective ancestral nodes. This can be seen from the relative magnitude of the entries, namely $a_{14}>\max(a_{12},a_{13})$ and $a_{24}>a_{23}$. To a lesser degree Alpha-Carotene also affects high intake of Lutein+Zeaxanthin as is seen from $a_{34}=0.281$.