Identifiability and estimation of recursive max-linear models

Abstract

We address the identifiablity and estimation of recursive max-linear structural equation models represented by an edge weighted directed acyclic graph (DAG). Such models are generally unidentifiable and we identify the whole class of DAGs and edge weights corresponding to a given observational distribution. For estimation, standard likelihood theory cannot be applied because the corresponding families of distributions are not dominated. Given the underlying DAG, we present an estimator for the class of edge weights and show that it can be considered a generalized maximum likelihood estimator. In addition, we develop a simple method for identifying the structure of the DAG. With probability tending to one at an exponential rate with the number of observations, this method correctly identifies the class of DAGs and, similarly, exactly identifies the possible edge weights.

MSC 2010 subject classifications: Primary 60E15, 62H12; secondary 62G05, 60G70, 62-09

Keywords and phrases: Causal inference, Bayesian network, directed acyclic graph, extreme value theory, generalized maximum likelihood estimation, graphical model, identifiability, max-linear model, structural equation model.

Introduction

Establishing and understanding cause-effect relations is an omnipresent desire in science and daily life. It is especially important when dealing with extreme events, because they are mostly dangerous and very costly; knowing and understanding the causes of such events and their causal relations could help us to deal better with them. Examples include incidents at airplane landings (Gissibl et al. ), flooding in river networks (Asadi et al. , Engelke and Hitz ), financial risk (Einmahl et al. ), and chemical pollution of rivers (Hoef et al. ). Such applications, where extreme risks may propagate through a network, have been the motivation behind the definition of recursive max-linear (ML) models in Gissibl and Klüppelberg . Recursive ML models are structural equation models (SEMs) represented by a directed acyclic graph (DAG) and thereby obey the basic Markov properties associated with directed graphical models (Lauritzen , Lauritzen et al. ). Both SEMs (see for example Bollen , Pearl ) and directed graphical models (see for example Koller and Friedman , Lauritzen , Spirtes et al. ) are well-established concepts for the understanding and quantification of causal inference from observational data. We note that Hitz and Evans and Engelke and Hitz discuss graphical models for extremes that are based on undirected graphs.

Recursive ML models are defined by a DAG, a collection of edge weights, and a vector of independent innovations. Important research problems that are addressed for recursive SEMs are the question of identifiability of the coefficients and the associated DAG from the observational distribution. Although the true DAG and edge weights for a recursive ML model are not identifiable from the observational distribution, the so-called max-linear coefficient matrix is identifiable and determines the possible class of DAGs and edge weights uniquely.

We shall show that estimation and structure learning of recursive ML models can be done in a simple and efficient fashion by exploiting properties of the ratios between observable components of the model. For a sufficiently large number of observations, these ratios identify the true ML coefficient matrix with a probability that converges exponentially fast to 1. For the situation where the DAG is known, we show that our estimator can be considered a maximum likelihood estimator in an extended sense, originally introduced by Kiefer and Wolfowitz .

Our paper is organized as follows. In Section 2 we introduce the model class of recursive ML models and the notation used throughout. In Section 3 we discuss the identifiability of a recursive ML model from its observational distribution. Here we show distributional properties of the ratio between two components. Based on these properties, we suggest an identification method. Section 4 is then devoted to the estimation of recursive ML models where we assume the DAG to be known. We show that the proposed estimates are generalized maximum likelihood estimates (GMLEs) in the sense of Kiefer–Wolfowitz. The main part is here the derivation of a specific Radon-Nikodym derivative. In Section 5 we complement the theoretical findings on the identifiability of recursive ML models with an efficient procedure to learn recursive ML models from observations only, even when the DAG itself is also unknown. Section 6 concludes and suggests further directions of research.

Preliminaries — recursive max-linear models

where ${\rm{pa}}(i)$ are the parents of node $i$ in $\mathcal{D}$. To highlight the DAG $\mathcal{D}$, we say that $\boldsymbol{X}$ follows a recursive ML model on $\mathcal{D}$. Note that this is a slight variation of the original definition in . We shall refer to $\boldsymbol{Z}=(Z_{1},\ldots,Z_{d})$ as the vector of innovations.

In the context of risk analysis, natural candidates for distributions of the innovations are extreme value distributions or distributions in their domain of attraction, resulting in a corresponding multivariate distribution (for details and background on multivariate extreme value models, see for example Beirlant et al. , de Haan and Ferreira , Resnick ).

Instead of $k\in{\rm{pa}}(i)$ we also write $k\to i$. Assigning the weight $d_{ji}(p)=\prod_{\nu=0}^{n-1}c_{k_{\nu}k_{\nu+1}}$ to every path $p=[j=k_{0}\to k_{1}\to\dots\to k_{n}=i]$ and denoting the set of all paths from $j$ to $i$ by $P_{ji}$, the non-negative matrix $B=(b_{ij})_{d\times d}$ with entries

is said to be the ML coefficient matrix of $\boldsymbol{X}$. This means for distinct $i,j\in V$, $b_{ji}$ is positive if and only if there is a path from $j$ to $i$; in that case $b_{ji}$ is the maximum weight of all paths from $j$ to $i$, where the weight of a path is the product of all edge weights $c_{ki}$ along this path. We say that a path from $j$ to $i$ whose weight equals $b_{ji}$ is max-weighted.

The components of $\boldsymbol{X}$ can also be expressed as max-linear functions of their ancestral innovations and an independent one; the corresponding ML coefficients are the entries of $B$:

The matrix product $\odot$ allows us to represent the ML coefficient matrix $B$ of $\boldsymbol{X}$ in terms of the weighted adjacency matrix $(c_{ij}\mathds{1}_{{\rm{pa}}(j)}(i))_{d\times d}$ of $\mathcal{D}$ since (2.2) and (2.3) simply become

Identifiability of a recursive max-linear model

In this section we discuss the question of identifiability of the elements of a recursive ML model from the distribution ${\mathcal{L}}(\boldsymbol{X})$ of $\boldsymbol{X}$. Indeed we shall show the following:

Let ${\mathcal{L}}(\boldsymbol{X})$ be the distribution of $\boldsymbol{X}$ following a recursive ML model. Then its ML coefficient matrix $B$ and the distribution of its innovation vector $\boldsymbol{Z}$ are identifiable from ${\mathcal{L}}(\boldsymbol{X})$. Furthermore, the class of all DAGs and edge weights that could have generated $\boldsymbol{X}$ by (2.1) can be obtained.

The remaining part of this section is devoted to proving Theorem 3.1, but first we shall consider a small example, illustrating the issues.

[The DAG and the edge weights are not necessarily identifiable] Consider a recursive ML model on the DAG $\mathcal{D}$ depicted below with edge weights $c_{12},c_{23},c_{13}$.

$1<span class="katex-display"><span class="katex"><span class="katex-mathml"><math xmlns="http://www.w3.org/1998/Math/MathML" display="block"><semantics><mrow><mn>2</mn></mrow><annotation encoding="application/x-tex">2</annotation></semantics></math></span><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.6444em;"></span><span class="mord">2</span></span></span></span></span>3$$\mathcal{D}$ According to (2.1), the components of $\boldsymbol{X}$ have the following representations

but also representations in terms of the innovations using (2.3) as

If $c_{13}\leq c_{12}c_{23}$ we have for any $c^{*}_{13}\in[0,c_{12}c_{23}]$ that $b_{13}=c_{12}c_{23}\vee c^{*}_{13}=c_{12}c_{23}\vee c_{13}=c_{12}c_{23}$; so we could also write

without changing the distribution ${\mathcal{L}}(\boldsymbol{X})$ of $\boldsymbol{X}$. This implies that if $c_{13}\leq c_{12}c_{23}$, $\boldsymbol{X}$ follows a recursive ML model on $\mathcal{D}$ with edge weights $c_{12},c_{23},c^{*}_{13}$ but it also follows a recursive model on the DAG $\mathcal{D}^{B}$ depicted below with edge weights $c_{12},c_{23}$.

$1<span class="katex-display"><span class="katex"><span class="katex-mathml"><math xmlns="http://www.w3.org/1998/Math/MathML" display="block"><semantics><mrow><mn>2</mn></mrow><annotation encoding="application/x-tex">2</annotation></semantics></math></span><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.6444em;"></span><span class="mord">2</span></span></span></span></span>3$$\mathcal{D}^{B}$ Consequently, we can neither identify $\mathcal{D}$ nor the value $c_{13}$ from the distribution ${\mathcal{L}}(\boldsymbol{X})$ of $\boldsymbol{X}$. However, note that the ML coefficient $b_{13}=c_{12}c_{23}\vee c_{13}$ is uniquely determined. If we however assume that $c_{13}>c_{12}c_{23}$, only $\mathcal{D}$ and the edge weights $c_{12},c_{23},c_{13}$ represent $\boldsymbol{X}$ in the sense of (2.1). Thus in this case the DAG and the edge weights are identifiable from the distribution ${\mathcal{L}}(\boldsymbol{X})$. $\Box$

As conclusion of Example 3.2, it is generally not possible to identify the true DAG $\mathcal{D}$ and the edge weights $c_{ki}$ underlying $\boldsymbol{X}$ in representation (2.1) from ${\mathcal{L}}(\boldsymbol{X})$, since several DAGs and edge weights may exist such that $\boldsymbol{X}$ has this representation. The smallest DAG of this kind is the DAG that has an edge $k\to i$ if and only if $k\to i$ is the only max-weighted path from $k$ to $i$. We call this DAG $\mathcal{D}^{B}$ the minimum ML DAG of $\boldsymbol{X}$ and note that this is uniquely determined from the ML coefficient matrix $B$. All other DAGs representing $\boldsymbol{X}$ are those that include the edges of $\mathcal{D}^{B}$ and whose nodes have the same ancestors. The edge weights $c_{ki}$ in the representation (2.1) of $\boldsymbol{X}$ are only uniquely determined for edges contained in $\mathcal{D}^{B}$; namely, by $b_{ki}$; otherwise, $c_{ki}$ may be any number in $(0,b_{ki}]$. We summarize these findings in the following theorem which is paraphrasing Theorems 5.3 and 5.4 of .

Suppose $\boldsymbol{X}$ follows a recursive ML model with edge weights $C=\{c_{ij}\}$ and ML coefficient matrix $B$. Let $\mathcal{D}^{B}$ be the minimum ML DAG of $\boldsymbol{X}$ as described above. Then a DAG $\mathcal{D}^{*}$ with associated weight matrix $C^{*}$ is a valid representation of $\boldsymbol{X}$ if and only if

$\mathcal{D}^{B}\subseteq\mathcal{D}^{*}$;

$\mathcal{D}^{*}$ and $\mathcal{D}^{B}$ have the same reachability matrix;

$c^{*}_{ij}=c_{ij}$ for $i\in{\rm{pa}}^{B}(j)$;

$c^{*}_{ij}\in(0,b_{ij}]$ for $i\in{\rm{pa}}^{*}(j)\setminus{\rm{pa}}^{B}(j)$,

where ${\rm{pa}}^{B}(j)$ and ${\rm{pa}}^{*}(j)$ denote the parents of $j$ in $\mathcal{D}^{B}$ and $\mathcal{D}^{*}$ respectively.

Based on the above observations, we investigate the identifiability of the whole class of DAGs and edge weights representing the max-linear structural equations (2.1) of $\boldsymbol{X}$ from ${\mathcal{L}}(\boldsymbol{X})$. Since this class can be recovered from $B$, it suffices to clarify whether $B$ is identifiable from ${\mathcal{L}}(\boldsymbol{X})$. There are many ways to prove that this is indeed the case. The way we present in this section suggests a simple procedure to estimate $B$ from independent realizations of $\boldsymbol{X}$ (see Algorithm 5.1 below). An alternative way can be found in Appendix 4.A.1 of .

which follows from the independence of the innovations and the fact that their distributions are atom-free.

where ${\rm{supp}}(Y_{ji})$ denotes the support of $Y_{ji}$.

In Table 3.1 we summarize the results of Lemma 3.4: depending on the relationship between $i$ and $j$ in $\mathcal{D}$, the support and atoms of $Y_{ji}$ are shown.

Table 3.1 and the fact that $b_{ji}=0$ for $j\not\in{\rm{An}}(i)$ (cf. (2.2)) suggest the following algorithm to find $B$ from ${\mathcal{L}}(\boldsymbol{X})$ since we can identify the support of $Y_{ji}$ from ${\mathcal{L}}(\boldsymbol{X})$. This proves the identifiability of $B$ from ${\mathcal{L}}(\boldsymbol{X})$. In fact, it is sufficient to know ${\rm{supp}}(Y_{ji})$ for all $i,j\in V$ with $i\neq j$ rather than the whole distribution ${\mathcal{L}}(\boldsymbol{X})$.

[Find $B$ from ${\mathcal{L}}(\boldsymbol{X})$]

For all $i\in V=\{1,\ldots,d\}$, set $b_{ii}=1$.

For all $i,j\in V$ with $i\neq j$, find ${\rm{supp}}(Y_{ji})$:

So far we have shown that the ML coefficient matrix $B$ of $\boldsymbol{X}$ can be obtained from ${\mathcal{L}}(\boldsymbol{X})$. Since all DAGs and edge weights that represent $\boldsymbol{X}$ in the sense of (2.1) can be determined from $B$, the only quantities we do not know about yet but appear in the definition of $\boldsymbol{X}$ are the innovations. In what follows we show that the distribution of the innovation vector $\boldsymbol{Z}$ is also identifiable from ${\mathcal{L}}(\boldsymbol{X})$. For this, due to the identifiability of $B$ from ${\mathcal{L}}(\boldsymbol{X})$ and the independence of the innovations, it suffices to provide an algorithm that determines the distributions of the innovations from ${\mathcal{L}}(\boldsymbol{X})$ and $B$. Note that $B$ also determines the ancestral relationships between any pair of nodes in that $j\in{\rm{An}}(i)$ for any DAG representing $\boldsymbol{X}$ if and only if $b_{ji}>0$.

We denote by $F_{Z_{i}}$ the distribution function of the innovation $Z_{i}$. For this algorithm, we do not have to know the whole distribution ${\mathcal{L}}(\boldsymbol{X})$; it is enough to know the ML coefficient matrix $B$ and the univariate marginal distribution functions of ${\mathcal{L}}(\boldsymbol{X})$.

Here we have used the convention that $\prod_{j\in\emptyset}a_{j}=1.$ The correctness of Algorithm 3.6 follows from the independence of the innovations and representation (2.3).

Estimation with known directed acyclic graph

In the following we let $\mathcal{B}(\mathcal{D})$ denote the class of possible ML coefficient matrices of all recursive ML models on $\mathcal{D}$. For $B$ being a matrix with non-negative entries and diagonal elements $b_{ii}=1$ we define $B_{0}:=(b_{ij}\mathds{1}_{{\rm{pa}}(j)}(i))_{d\times d}$. Then it holds that $B\in\mathcal{B}(\mathcal{D})$ if and only if $B$ satisfies the following

[Illustration of (4.1)] To illustrate the above, consider the small network below

$1<span class="katex-display"><span class="katex"><span class="katex-mathml"><math xmlns="http://www.w3.org/1998/Math/MathML" display="block"><semantics><mrow><mn>2</mn></mrow><annotation encoding="application/x-tex">2</annotation></semantics></math></span><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.6444em;"></span><span class="mord">2</span></span></span></span></span>3$$4$ and a potential ML coefficient matrix $B$ with reduction $B_{0}$, as given below.

where we have used that $1$ and $2$ are not ancestors of $3$ and $1$ is not a parent of $4$. We wish to check whether $B\in B(\mathcal{D})$ for this particular DAG so we further calculate

Now $B=I_{4}\vee(B\odot B_{0})$ readily implies that $b_{ii}=1,i=1,\ldots,4$ and $b_{14}=b_{12}b_{24}$.

A simple estimate of 𝑩𝑩\boldsymbol{B}

Next we discuss a sensible estimate of $B$. Table 3.1 shows that for $j\in{\rm{an}}(i)$ the minimal value that can be observed for the ratio $Y_{ji}=X_{i}/X_{j}$ is $b_{ji}$, which is an atom of $Y_{ji}$. This suggests the following estimate $\breve{B}$ of the ML coefficient matrix:

Davis and Resnick suggested such minimal observed ratios as estimates for parameters in max-ARMA processes. For $n$ sufficiently large, we can expect to observe the atoms $b_{ji}$ for $j\in{\rm{an}}(i)$ in the sample $\boldsymbol{x}^{(1)},\ldots,\boldsymbol{x}^{(n)}$ and, hence, to estimate the ML coefficients exactly. However, if $n$ is not large we may with positive probability have that $\breve{B}$ is not an ML coefficient matrix of any recursive ML model on $\mathcal{D}$ as the following simple example shows:

[$\breve{B}$ is not necessarily in ${\mathcal{B}}(\mathcal{D})$] Consider the DAG

$1<span class="katex-display"><span class="katex"><span class="katex-mathml"><math xmlns="http://www.w3.org/1998/Math/MathML" display="block"><semantics><mrow><mn>2</mn></mrow><annotation encoding="application/x-tex">2</annotation></semantics></math></span><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.6444em;"></span><span class="mord">2</span></span></span></span></span>3$$\mathcal{D}$ and assume we observe $\breve{b}_{31}>\breve{b}_{32}\breve{b}_{21}$. Then the matrix $\breve{B}$ fails to satisfy (4.1) and hence is not an element of ${\mathcal{B}}(\mathcal{D})$. $\Box$

However, if we only estimate the ML coefficients corresponding to edges in $\mathcal{D}$ and then compute an estimate based on Lemma 4.3 below this phenomenon cannot occur.

if and only if $A=(I_{d}\vee B_{0})^{\odot(d-1)}$.

We first show that $A=(I_{d}\vee B_{0})^{\odot(d-1)}$ satisfies (4.2). It is immediate that $a_{ji}>0\iff j\in{\rm{An}}(i)$. We have (, Proposition 1.6.10) that

It is easy to see directly that $B_{0}^{\odot k}=0$ for $k\geq d$ and hence if $\check{A}$ is a solution to (4.2) we get by iteration, using that $(M\vee N)\odot K=(M\odot K)\vee(N\odot K)$,

and hence the solution to the equation is unique. ∎

Thus we may define the estimate $\widehat{B}$ by first calculating the matrix $\breve{B}_{0}=(\breve{b}_{ij}\mathds{1}_{{\rm{pa}}(j)}(i))_{d\times d}$ and then iterating the $\odot$-matrix product as:

It then follows that $\widehat{B}_{0}=\breve{B}_{0}$ and Lemma 4.3 yields that $\widehat{B}$ is the unique element of $\mathcal{B}(\mathcal{D})$ satisfying (4.3). By Lemma 3.4(b), we also have

Consequently, when using $\widehat{B}$ or $\breve{B}$ as an estimate of $B$, we never underestimate a ML coefficient; furthermore, the matrix $\widehat{B}$ always estimates $B$ more precisely than $\breve{B}$ and since we always have $\widehat{B}\in{\mathcal{B}}(\mathcal{D})$, $\widehat{B}$ seems to be clearly preferable as an estimate of $B$.

The following example shows how effective the estimate $\widehat{B}$ can be; in particular, $n$ does not necessarily need to be large.

[One observation may be enough to estimate $B$ exactly] Consider the DAG

$1<span class="katex-display"><span class="katex"><span class="katex-mathml"><math xmlns="http://www.w3.org/1998/Math/MathML" display="block"><semantics><mrow><mn>2</mn></mrow><annotation encoding="application/x-tex">2</annotation></semantics></math></span><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.6444em;"></span><span class="mord">2</span></span></span></span></span>3<span class="katex-display"><span class="katex"><span class="katex-mathml"><math xmlns="http://www.w3.org/1998/Math/MathML" display="block"><semantics><mrow><mn>4</mn></mrow><annotation encoding="application/x-tex">4</annotation></semantics></math></span><span class="katex-html" aria-hidden="true"><span class="base"><span class="strut" style="height:0.6444em;"></span><span class="mord">4</span></span></span></span></span>\mathcal{D}$ and assume that the paths $[1\to 2\to 4]$ and $[1\to 3\to 4]$ are both max-weighted, which is equivalent to $b_{12}b_{24}=b_{13}b_{34}$. If we observe the event

Let \boldsymbol{X}^{(t)}=\big{(}X_{1}^{(t)},\ldots,X_{n}^{(t)}\big{)} for $t=1,\ldots,n$ be a sample from a recursive ML model on a DAG $\mathcal{D}$ with ML coefficient matrix $B$. Let $i\in V$ and $k\in{\rm{pa}}(i)$. It then holds that

First note that the events $\{X_{i}=b_{ki}X_{k}\}$ and $\{X_{i}>b_{ki}X_{k}\}$ are complementary and both have positive probability. Further, using that $\boldsymbol{X}^{(1)},\ldots,\boldsymbol{X}^{(n)}$ are independent and identically distributed yields

In conclusion, $\widehat{B}$ has the nice property to be ’geometrically consistent’ in the sense that the probability of $\{\widehat{B}=B\}$ converges exponentially fast to one.

The matrix 𝑩^bold-^𝑩\boldsymbol{\widehat{B}} is a generalized maximum likelihood estimate

As we found in the previous section, the estimate $\widehat{B}$ is preferable to the direct estimate $\breve{B}$ as it will always be closer to the true value. In this section we further establish that $\widehat{B}$ is not just an ad hoc estimator, but can indeed be derived from likelihood considerations.

For $B\in{\mathcal{B}}(\mathcal{D})$ and a fixed distribution of the innovation vector we let $P_{B}$ denote the probability measure induced by a recursive ML model on $\mathcal{D}$ with ML coefficient matrix $B$, i.e. the distribution of $\boldsymbol{X}$ where $\boldsymbol{X}=\boldsymbol{Z}\odot B$. We shall denote the family of these probability measures by ${\mathcal{P}}(\mathcal{D})$.

We cannot use standard maximum likelihood methods to estimate $B$, since the family ${\mathcal{P}}(\mathcal{D})$ is not dominated (cf. Example 4.4.1 of ) and hence the standard likelihood function is not well defined. However, there exist generalizations of maximum likelihood estimation (GMLE) that cover the undominated case as well; Kalbfleisch and Prentice , Kiefer and Wolfowitz , and Scholz suggested such extensions. We essentially follow the Kiefer–Wolfowitz definition of a GMLE as also done, for example, by Gill et al. and Johansen . In the following we shall show that $\widehat{B}$ can be seen as a maximum likelihood estimate of $B$ in the extended sense introduced by Kiefer and Wolfowitz in .

where ${dP}/{d(P+Q)}$ denotes a density of $P$ with respect to $P+Q$. Then we call $\widehat{P}$ a generalized maximum likelihood estimate of $P_{0}$ if

We begin with an example that shall help to get an idea and provide insights into the concepts and arguments we shall use in the general case. It is deliberately very detailed and although it deals with a very special case, it illustrates the main issues also for the general case.

[How to find a density and the associated GMLEs] For $B,B^{*}\in{\mathcal{B}}(\mathcal{D})$ where $\mathcal{D}=(\{1,2\},1\to 2)$, we show that the partition

We now use the density found to determine the GMLE of $B$. The only ML coefficient we have to estimate is $b_{12}$. As before we let $\widehat{b}_{12}=\breve{b}_{12}$ be the minimal observed ratio of ${X_{2}}/{X_{1}}$ and let $\widehat{B}$ be the corresponding ML coefficient matrix from (4.3). Defining $n(B,B^{*})=|\{t:\boldsymbol{x}^{(t)}\in A_{1/2}(B,B^{*})\}|$ and using that $n(B,B^{*})=n(B^{*},B)$, we obtain

Let now $\widetilde{B}$ be an arbitrary potential GMLE of $B$. Then $P_{\widetilde{B}}\in{\mathcal{P}}(\mathcal{D})$ satisfies the first condition in (4.4) if and only if

In summary, some $\widetilde{B}\in\mathcal{B}(\mathcal{D})$ is a GMLE of $B$ if and only if (4.7) and (4.8) are satisfied. We discuss the possible GMLEs of $b_{12}$ in detail.

$\widetilde{b}_{12}>\widehat{b}_{12}$ is no GMLE: This follows directly from (4.7). Figure 2(b) shows a situation that contradicts (4.8), similarly to Figure 2(a) in (1).

In what follows we specify, for the general case, one density of $P_{B}$ with respect to $P_{B}+P_{B^{*}}$ that has a representation as in (4.6) and leads to $\widehat{B}$ as a GMLE of $B$.

Let $B,B^{*}\in{\mathcal{B}}(\mathcal{D})$ and define

is a density of $P_{B}$ with respect to $P_{B}+P_{B^{*}}$.

We observe an interesting relation between the density (4.11) for $\mathcal{D}$ and corresponding densities for subgraphs of $\mathcal{D}$.

[Local densities $\rho_{i}$] Consider the DAGs

This can be observed from Figure 3, where the densities are depicted as functions of ${x_{2}}/{x_{1}}$ and/or ${x_{3}}/{x_{2}}$ for all nine different orders between the ML coefficients in $B$ and $B^{*}$.

Conversely, $\rho_{2}$ and $\rho_{3}$ can be derived from $\rho$ as follows:

which we learn from Figure 3 again. $\Box$

We now extend the findings from Example 4.9 to the general case. Furthermore, we show that the densities $\rho_{i}$ are densities of regular conditional distributions.

Let $B,B^{*}\in{\mathcal{B}}(\mathcal{D})$ and let $\boldsymbol{X}=\boldsymbol{Z}\odot B,\boldsymbol{X}^{*}=\boldsymbol{Z}\odot B^{*}$ follow corresponding recursive ML models on $\mathcal{D}$. For $i\in V$, let $\rho_{i}$ be the density given in (4.11) with respect to the DAG $\mathcal{D}_{i}=({\rm{Pa}}(i),\{(k,i):k\in{\rm{pa}}(i)\})$ as well as $B_{i}$ and $B_{i}^{*}$ the ML coefficient matrices of recursive ML models on $\mathcal{D}_{i}$ with edge weights $c_{ki}=b_{ki}$ and $c^{*}_{ki}=b^{*}_{ki}$, respectively.

We have for $\rho(\boldsymbol{x},B,B^{*})$ given in (4.11)

The function $\rho_{i}$ can be computed from $\rho$ by

where we set $\min_{\boldsymbol{y}\in\emptyset}\rho(\boldsymbol{y},B,B^{*})=0$.

Next, we show that $\widehat{B}$ is indeed a GMLE in the sense of . Note also that the GMLE is obtained by piecing together individual GMLEs corresponding to conditional distributions of any variable given its parents. Thus this is similar to what is obtained in cases where the distributions have densities with respect to a product measure, as the maximum of the likelihood function is then obtained by maximizing each conditional likelihood function for the density of a node given its parents.

Let \boldsymbol{x}^{(t)}=\big{(}x_{1}^{(t)},\ldots,x_{n}^{(t)}\big{)} for $t=1,\ldots,n$ be a sample from a recursive ML model on a DAG $\mathcal{D}$ with ML coefficient matrix $B\in{\mathcal{B}}(\mathcal{D})$ unknown.

The matrix $\widehat{B}$ from (4.3) is a GMLE of $B$.

For every $i\in V$, $(\widehat{b}_{ki},k\in{\rm{pa}}(i))$ is a GMLE of the ML coefficients $(b_{ki},k\in{\rm{pa}}(i))$ of a random vector following a recursive ML model on $\mathcal{D}_{i}=({\rm{Pa}}(i),\{(k,i):k\in{\rm{pa}}(i)\})$ with edge weights $c_{ki}=b_{ki}$.

For every $i\in V$ and $k\in{\rm{pa}}(i)$, $\widehat{b}_{ki}$ is the only GMLE of the ML coefficient $b_{ki}$ of a random vector following a recursive ML model on $\mathcal{D}_{ki}=(\{k,i\},\{(k,i)\})$ with edge weight $c_{ki}=b_{ki}$.

Hence, $x^{(t_{2})}_{i}=b_{ki}x^{(t_{2})}_{k}$ for some $k\in{\rm{pa}}(i)$ with $\widehat{b}_{ki}<b_{ki}$. Let now $t_{1}\in\{1,\ldots,n\}$ such that $\bigwedge_{s=1}^{n}y^{(s)}_{ki}=y^{(t_{1})}_{ki}$. As $\widehat{b}_{ki}=\bigwedge_{s=1}^{n}y^{(s)}_{ki}$, we have $x_{i}^{(t_{1})}<b_{ki}x_{k}^{(t_{1})}$ implying that $\boldsymbol{x}^{(t_{1})}\in A_{0}(B,\widehat{B})$. The statement in (b) is a consequence of (a), and (c) has already been shown in Example 4.6. ∎

Figure 4 illustrates the DAGs $\mathcal{D}_{i}$ in Theorem 4.11(b) or Proposition 4.10.

Learning the structure of a recursive max-linear model

In contrast to the assumptions in the previous section, we now assume independent realizations $\boldsymbol{x}^{(1)},\ldots,\boldsymbol{x}^{(n)}$ of $\boldsymbol{X}$ following a recursive ML model but the underlying DAG $\mathcal{D}$ is unknown. We know from previous discussions that it is not possible to recover $\mathcal{D}$ and the true edge weights $c_{ki}$, and we therefore again focus on the estimation of $B$.

Following Algorithm 3.5, it suffices for any pair of distinct $i,j\in V$ to decide whether ${\rm{supp}}(Y_{ji})={\rm{supp}}(X_{i}/X_{j})$ has a positive lower bound, alternatively a finite upper bound, and if so, to estimate the bound. Recall from Table 3.1 that, if there is such a bound, then it is an atom of $Y_{ji}$. Since we can expect to observe atoms more than twice for $n$ sufficiently large, we propose the following estimation method.

[Find an estimate $\widecheck{B}$ of $B$ from $\boldsymbol{x}^{(1)},\ldots,\boldsymbol{x}^{(n)}$]

For all $i\in V=\{1,\ldots,d\}$, set $\widecheck{b}_{ii}=1$.

if $\#\left\{t:\bigwedge_{s=1}^{n}y_{ji}^{(s)}=y_{ji}^{(t)}\right\}\geq 2$, then conclude $j\in{\rm{an}}(i)$, set $\widecheck{b}_{ji}=\bigwedge_{t=1}^{n}y_{ji}^{(t)}$;

The second item summarizes two steps: the first is concerned with estimating the ancestors of the nodes, the second with estimating the ML coefficients.

Note that the estimate $\widecheck{B}$ from Algorithm 5.1 is not necessarily a ML coefficient matrix of a recursive ML model. For example, the property that $b_{ji}>0$ if $b_{jk}b_{ki}>0$ (see, for example, Corollary 3.12 of ) is not guaranteed. Many modifications of $\widecheck{B}$ are possible, and here we shall not discuss this in detail. Rather we notice that the probability that Algorithm 5.1 outputs the true ML coefficient matrix $B$ tends to one as $n\to\infty$. As in the case where the DAG is known — see Proposition 4.5 — this probability converges to one at an exponential rate.

Conclusion and outlook

We studied the identifiability of the elements of a recursive ML model from the distribution ${\mathcal{L}}(\boldsymbol{X})$ of $\boldsymbol{X}$. The associated DAG and the edge weights are not identifiable, however, the ML coefficient matrix $B$ is. In other words, we can identify the representation (2.3) but not (2.1). The class of all DAGs and edge weights that could have generated $\boldsymbol{X}$ via (2.1) and the distribution of the innovation vector are identifiable from ${\mathcal{L}}(\boldsymbol{X})$. As a consequence, we can recover $B$, the class of the DAGs and edge weights, and the innovation distributions from realizations of $\boldsymbol{X}$.

We have shown that $\hat{B}$ is a generalized maximum likelihood estimate. This is primarily of theoretical interest as it shows the estimate is not purely based on an ad hoc procedure. However, it opens up the possibility of going further, using likelihood theory, for example to study issues of likelihood ratio testing of hypothesis for specific values of the coefficients, or even for the presence or absence of edges in the underlying graph.

Parameter estimation and structure learning for recursive ML models seem to be challenging tasks because assumptions usually made in standard methods are not met. However, in both cases, $B$ can be estimated by a simple procedure. The key idea of our approach is to consider the observed ratios between any pair of components, i.e. to perform a transformation on the realizations. The transformed realizations or rather the distributional properties of the corresponding random variables make it possible to identify, with probability 1, the true $B$ whenever the number of observations $n$ is sufficiently large. It would be interesting to investigate the relationship between the performance of our procedures and the number $n$ of observations. Here, one possible question is how many observations are at least necessary to estimate $B$ exactly; see, Example 4.4. In addition it would be interesting to study estimation of the DAG structure for moderate sample sizes, where exact estimation is not guaranteed.

An important goal for future work is to apply the procedures to real-world data. However, it is unreasonable to expect any non-simulated data to follow a recursive ML model exactly, and the model should then be modified by adding appropriate noise terms. In particular we should not expect that we observe a minimal observed ratio more than twice, as we exploit in Algorithm 5.1. It seems to be more reasonable to expect values close to each other. We therefore want to develop methods based on accumulation points. It is hard to imagine noise models that would lead to simple exact likelihood analysis. One should then rather study the asymptotic precision of reasonable estimates and their behaviour under appropriate scaling, for example along the lines of .

Acknowledgements

We thank Justus Hartl for providing a first discussion about the different estimators suggested in this paper in his master’s thesis. NG acknowledges support by Deutsche Forschungsgemeinschaft (DFG) through the TUM International Graduate School of Science and Engineering (IGSSE). All authors benefited from financial support from the Alexander von Humboldt Stiftung.

References

Appendix A Appendix: some technical proofs

The proof is by induction on the number of nodes of $\mathcal{D}$. For $d=1$ the statement is clear. Assume now that $\mathcal{D}=(V,E)$ has $d+1$ nodes and that the assertion holds with respect to DAGs with at most $d$ nodes. Furthermore, assume without loss of generality that $d+1$ is a terminal node (i.e., ${\rm{de}}(d+1)=\emptyset$). Since $(X_{1},\ldots,X_{d})$ follows a recursive ML model on the DAG $(\{1,\ldots,d\},E\cap(\{1,\ldots,d\}\times\{1,\ldots,d\}))$ with ML coefficient matrix $B=(b_{ij})_{d\times d}$ and $B^{*}=(b^{*}_{ij})_{d\times d}$ is the ML coefficient matrix of a recursive ML model on this DAG as well, the induction hypothesis yields that

For every $i\in V$ we have by (2.3) on $\Omega_{i}$ that

Noting from the proof of Theorem 4.2 of that

we obtain from (A.2) on $\bigcap_{i=1}^{d}\Omega_{i}$,

From (3.1) we then finally observe that \bigcap_{i=1}^{d}\Omega_{i}\cap\big{(}\Omega_{1/2}^{1,d+1}\cup\Omega_{1/2}^{2,d+1}\big{)} and $\bigcap_{i=1}^{d}\Omega_{i}\cap\Omega_{d+1}$ only differ by a set of probability zero, and, hence, (4.10) follows from (A.1). ∎

Proof of Theorem 4.8

We must verify properties (A)–(C) of (4.5).

(A) Since $V$ is finite, it suffices to show for every $i\in V$,

The former is immediate by (4.9). By the same argument we have for the latter,

where we have used (2.3) and (3.1) for the last inequality and equality, respectively. Thus we have verified (A).

(C) We observe from the definition of $A_{0}(B,B^{*})$ and $A_{1/2}(B,B^{*})$ that

Since $A_{0}(B^{*},B)$ is a $P_{B^{*}}$-null set by (A), this holds for the subset $A_{1}(B,B^{*})$ as well. ∎

Proof of Proposition 4.10

Denoting by $A^{i}_{0}(B_{i},B^{*}_{i})$, $A^{i}_{1/2}(B_{i},B^{*}_{i})$, $A^{i}_{1}(B_{i},B^{*}_{i})$ the sets defining $\rho_{i}(\cdot,B_{i},B^{*}_{i})$, we have for the corresponding sets of $\rho$,

From this we obtain (a) and (b). Now, to see (c) we reason as follows:

is a regular conditional distribution function of $X_{i}$ given $\boldsymbol{X}_{{\rm{pa}}(i)}$. To see this, use (4.9) and the independence of the innovations to obtain

Since $\boldsymbol{X}$ and $\boldsymbol{X}^{*}$ share the same innovation vector, we have

and for this again by definition of $\rho_{i}$ (cf. (4.6) and the related discussion) that

Since $F_{Z_{i}}$ is atom-free, this can be read directly from Figure 5.