Identifiability and estimation of recursive max-linear models

TL;DR

This paper investigates the identifiability and estimation methods for recursive max-linear models on DAGs, providing new techniques for structure learning and parameter estimation despite non-identifiability issues.

Contribution

It characterizes the class of DAGs and edge weights compatible with observed distributions and proposes a generalized MLE and a structure identification method.

Findings

The class of DAGs and edge weights can be fully identified from observational data.

A generalized maximum likelihood estimator for edge weights is developed.

A simple method reliably identifies DAG structures with high probability.

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 structures of the DAGs. 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.

Tables

Table 1. Table 3.1: Distributional properties of Y j ​ i subscript 𝑌 𝑗 𝑖 Y_{ji} for distinct i , j ∈ V 𝑖 𝑗 𝑉 i,j\in V .

Relationship between $i$ and $j$	$\mathrm{supp}\left(Y_{ji}\right)$	Atoms
$j\in \mathrm{an}\left(i\right)$	$[b_{ji},\infty )$	$\left\{b_{\ell i}/b_{\ell j},\ell \in \mathrm{An}\left(j\right)\right\}$
$i\in \mathrm{an}\left(j\right)$	$(0,1/b_{ij}]$	$\left\{b_{\ell i}/b_{\ell j},\ell \in \mathrm{An}\left(i\right)\right\}$
otherwise:
if $\mathrm{an}\left(i\right)\cap \mathrm{an}\left(j\right)\ne \varnothing$	${ℝ}_{+}$	$\left\{b_{\ell i}/b_{\ell j},\ell \in \mathrm{an}\left(i\right)\cap \mathrm{an}\left(j\right)\right\}$
if $\mathrm{an}\left(i\right)\cap \mathrm{an}\left(j\right)=\varnothing$	${ℝ}_{+}$	$\varnothing$

Equations

$$X_i=\underset{k\in \mathrm{pa}(i)}{\bigvee }{c}_{ki}{X}_k\vee Z_i,i=1,\dots ,d,$$

$$b_{ji}=\underset{p\in P_{ji}}{\bigvee }{d}_{ji}(p)\text{for }j\in \mathrm{an}(i),b_{ii}=1,\text{and}{b}_{ji}=0\text{for }j\in V\setminus \mathrm{An}(i),$$

$$X_i=\underset{j=1}{\overset{d}{\bigvee }}{b}_{ji}{Z}_j=\underset{j\in \mathrm{An}(i)}{\bigvee }{b}_{ji}{Z}_j,i=1,\dots ,d;$$

$$\displaystyle(F=(f_{ij})_{m\times n},G=(g_{ij})_{n\times p})\mapsto F\odot G:=\Big{(}\bigvee\limits_{k=1}^{n}f_{ik}g_{kj}\Big{)}_{m\times p},$$

$$B=(I_d\vee C{)}^{\odot (d-1)}=\underset{k=0}{\overset{d-1}{\bigvee }}{C}^{\odot k},\mathbf{X}=\mathbf{Z}\odot B,$$

$$X_1=Z_1,X_2=c_{12}{X}_1\vee Z_2,\text{and}{X}_3=c_{13}{X}_1\vee c_{23}{X}_2\vee Z_3.$$

$$X_1=Z_1,X_2=c_{12}{Z}_1\vee Z_2,\text{and}{X}_3=(c_{12}{c}_{23}\vee c_{13})Z_1\vee c_{23}{Z}_2\vee Z_3,$$

$$X_3=c_{13}^{\ast }{X}_1\vee c_{23}{X}_2\vee Z_3$$

$$\big{\{}Z_{i}=xZ_{j}\big{\}}$$

$$\displaystyle{\rm{supp}}(Y_{ji})=\begin{cases}[b_{ji},\infty)&\text{if $j\in{\rm{an}}(i)$}\\ \big{(}0,{1}/{b_{ij}}\big{]}&\text{if $j\in{\rm{de}}(i)$}\\ \mathbb{R}_{+}&\text{otherwise},\end{cases}$$

$$\displaystyle\Big{\{}\bigvee_{\ell\in{\rm{An}}(i)\cap{\rm{An}}(j):\atop b_{\ell i}=b_{\ell j}x}b_{\ell i}Z_{\ell}>\bigvee_{\ell\in{\rm{An}}(i)\cap{\rm{An}}(j):\atop b_{\ell i}\neq b_{\ell j}x}(b_{\ell i}\vee xb_{\ell j})Z_{\ell}\vee\bigvee_{\ell\in{\rm{An}}(i)\setminus{\rm{An}}(j)}b_{\ell i}Z_{\ell}\vee\bigvee_{\ell\in{\rm{An}}(j)\setminus{\rm{An}}(i)}xb_{\ell j}Z_{\ell}\Big{\}}$$

$$\mathrm{supp}(Y_{ji})=\left\{\frac{{\bigvee }_{\ell \in \mathrm{An}(i)}{b}_{\ell i}{z}_{\ell }}{{\bigvee }_{\ell \in \mathrm{An}(j)}{b}_{\ell j}{z}_{\ell }}:{\mathbf{z}}_{\mathrm{An}(i)\cup \mathrm{An}(j)}\in {\mathbb{R}}_{+}^{|\mathrm{An}(i)\cup \mathrm{An}(j)|}\right\}.$$

$${\mathbb{R}}_{+}^{|\mathrm{An}(i)\cup \mathrm{An}(j)|}\to {\mathbb{R}}_{+},{\mathbf{z}}_{\mathrm{An}(i)\cup \mathrm{An}(j)}\mapsto \frac{{\bigvee }_{\ell \in \mathrm{An}(i)}{b}_{\ell i}{z}_{\ell }}{{\bigvee }_{\ell \in \mathrm{An}(j)}{b}_{\ell j}{z}_{\ell }}$$

$$\underset{\ell \in \mathrm{An}(i)\cap \mathrm{An}(j)}{\bigvee }a{b}_{\ell j}{z}_{\ell }\vee \underset{\ell \in \mathrm{An}(j)\setminus \mathrm{An}(i)}{\bigvee }a{b}_{\ell j}{z}_{\ell }\le \underset{\ell \in \mathrm{An}(i)}{\bigvee }{b}_{\ell i}{z}_{\ell }$$

$$F_{Z_i}(x)=\frac{\mathbb{P}(X_i\le x)}{{\prod }_{j\in \mathrm{an}(i)}{F}_{Z_j}(x/b_{ji})}.$$

$$[b_{ji}>0⟺j\in \mathrm{An}(i)]\text{ and }B=I_d\vee (B\odot B_0);$$

$$B=\left(\begin{array}{cccc}{b}_{11}& b_{12}& 0& b_{14}\\ 0& b_{22}& 0& b_{24}\\ 0& 0& b_{33}& b_{34}\\ 0& 0& 0& b_{44}\end{array}\right)B_0=\left(\begin{array}{cccc}0& b_{12}& 0& 0\\ 0& 0& 0& b_{24}\\ 0& 0& 0& b_{34}\\ 0& 0& 0& 0\end{array}\right),$$

$$I_4\vee (B\odot B_0)=I_4\vee \left(\begin{array}{cccc}0& b_{11}{b}_{12}& 0& b_{12}{b}_{24}\\ 0& 0& 0& b_{22}{b}_{24}\\ 0& 0& 0& b_{33}{b}_{34}\\ 0& 0& 0& 0\end{array}\right)=\left(\begin{array}{cccc}1& b_{11}{b}_{12}& 0& b_{12}{b}_{24}\\ 0& 1& 0& b_{22}{b}_{24}\\ 0& 0& 1& b_{33}{b}_{34}\\ 0& 0& 0& 1\end{array}\right)$$

$${\breve{b}}_{ii}=1,{\breve{b}}_{ji}=0\text{ for }j\in V\setminus \mathrm{An}(i),\text{and}{\breve{b}}_{ji}=\underset{t=1}{\overset{n}{\bigwedge }}{y}_{ji}^{(t)}=\underset{t=1}{\overset{n}{\bigwedge }}\frac{x_i^{(t)}}{x_j^{(t)}}\text{ for }j\in \mathrm{an}(i).$$

$$[a_{ji}>0⟺j\in \mathrm{An}(i)]\text{ and }A=I_d\vee (A\odot B_0)$$

$$(I_d\vee B_0{)}^{\odot (d-1)}=\underset{k=0}{\overset{d-1}{\bigvee }}{B}_0^{\odot k}=\underset{k=0}{\overset{\infty }{\bigvee }}{B}_0^{\odot k}$$

$$I_d\vee (A\odot B_0)=I_d\vee \{(I_d\vee B_0{)}^{\odot (d-1)}\odot B_0\}=I_d\vee \{\underset{k=1}{\overset{\infty }{\bigvee }}{B}_0^{\odot k}\}=\underset{k=0}{\overset{\infty }{\bigvee }}{B}_0^{\odot k}=A.$$

$$\check{A}$$

$$\hat{B}=(I_d\vee {\breve{B}}_0{)}^{\odot (d-1)}.$$

$$b_{ji}\le {\hat{b}}_{ji}\le {\breve{b}}_{ji}\text{ for }j\in \mathrm{an}(i).$$

$$\displaystyle\big{\{}X_{2}=b_{12}X_{1}\big{\}}\cap\big{\{}X_{3}=b_{13}X_{1}\big{\}}\cap\big{\{}X_{4}=b_{24}X_{2}\big{\}}\cap\big{\{}X_{4}=b_{34}X_{3}\big{\}},$$

$$\mathbb{P}\left(\underset{t=1}{\overset{n}{\bigwedge }}{Y}_{ki}^{(t)}=b_{ki}\right)\ge 1-p\text{ for some p∈(0,1) }$$

$$n\ge \frac{\ln (p)}{\ln (\mathbb{P}(Y_{ki}>b_{ki}))}.$$

$$\displaystyle{\mathbb{P}}\Big{(}\bigwedge_{t=1}^{n}Y_{ki}^{(t)}=b_{ki}\Big{)}=1-{\mathbb{P}}\Big{(}\bigwedge_{t=1}^{n}Y_{ki}^{(t)}>b_{ki}\Big{)}=1-\prod_{t=1}^{n}{\mathbb{P}}({Y_{ki}^{(t)}}>b_{ki})=1-{\mathbb{P}}({Y_{ki}}>b_{ki})^{n}.$$

$$\rho (\mathbf{x},P,Q):=\frac{dP}{d(P+Q)}(\mathbf{x}),$$

Full text

Identifiability and estimation of recursive max-linear models

Nadine Gissibl Center for Mathematical Sciences, Technical University of Munich, 85748 Garching, Boltzmannstrasse 3, Germany; e-mail: n.gissibl@tum.de, cklu@tum.de   

Claudia Klüppelberg11footnotemark: 1

  

Steffen Lauritzen Department of Mathematical Sciences, University of Copenhagen, Universitetsparken 5, 2100 Copenhagen, Denmark; e-mail: lauritzen@math.ku.dk

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.

1 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. [13]), flooding in river networks (Asadi et al. [1], Engelke and Hitz [9]), financial risk (Einmahl et al. [8]), and chemical pollution of rivers (Hoef et al. [15]). 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 [12]. 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 [21], Lauritzen et al. [22]). Both SEMs (see for example Bollen [3], Pearl [23]) and directed graphical models (see for example Koller and Friedman [19], Lauritzen [20], Spirtes et al. [27]) are well-established concepts for the understanding and quantification of causal inference from observational data. We note that Hitz and Evans [14] and Engelke and Hitz [9] 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 [18].

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.

2 Preliminaries — recursive max-linear models

In this section we introduce notation and summarize the most important properties of recursive ML models needed. A recursive ML model for a random vector $\boldsymbol{X}=(X_{1},\ldots,X_{d})$ is specified by an underlying structure in terms of a DAG $\mathcal{D}$ with nodes $V=\{1,\ldots,d\}$, positive edge weights $c_{ki}$ for $i\in V$ and $k\in{\rm{pa}}(i)$, and independent positive random variables $Z_{1},\dots,Z_{d}$ with support ${\mathbb{R}}_{+}:=(0,\infty)$ and atom-free distributions:

[TABLE]

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 [12]. 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. [2], de Haan and Ferreira [7], Resnick [24, 25]).

Throughout the paper we use the following notation. The sets ${\rm{an}}(i)$, ${\rm{pa}}(i)$, and ${\rm{de}}(i)$ contain the ancestors, parents, and descendants of node $i$ in $\mathcal{D}$. We set ${\rm{An}}(i)={\rm{an}}(i)\cup\{i\}$ and ${\rm{Pa}}(i)={\rm{pa}}(i)\cup\{i\}$. For $U\subsetneq V$ we write $\boldsymbol{X}_{U}=(X_{\ell},\ell\in U)$ and accordingly for $\boldsymbol{x}\in{\mathbb{R}}_{+}^{d}$, $\boldsymbol{x}_{U}=(x_{\ell},\ell\in U)$.

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

[TABLE]

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$:

[TABLE]

see Theorem 2.2 of [12].

For two non-negative matrices $F$ and $G$, where the number of columns in $F$ is equal to the number of rows in $G$, we define the matrix product $\odot:\overline{\mathbb{R}}_{+}^{m\times n}\times\overline{\mathbb{R}}_{+}^{n\times p}\rightarrow\overline{\mathbb{R}}_{+}^{m\times p}$ by

[TABLE]

where $\overline{\mathbb{R}}_{+}=[0,\infty)$. The triple $(\overline{\mathbb{R}}_{+},\vee,\cdot)$, is an idempotent semiring with [math] as 0-element and $1$ as 1-element and the operation $\odot$ is therefore a matrix product over this semiring; see for example Butkovič [4]. Denoting by ${\mathcal{M}}$ all $d\times d$ matrices with non-negative entries and by $\vee$ the componentwise maximum between two matrices, $({\mathcal{M}},\vee,\odot)$ is also a semiring with the null matrix as 0-element and the $d\times d$ identity matrix $I_{d}$ as 1-element.

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

[TABLE]

where we have let $A^{\odot 0}=I_{d}$ and $A^{\odot k}=A^{\odot(k-1)}\odot A$ for $A\in\overline{\mathbb{R}}_{+}^{d\times d}$ and $k\in\mathbb{N}$; see Proposition 1.6.15 of Butkovič [4] as well as Theorem 2.4 and Corollary 2.5 of [12].

3 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:

Theorem 3.1**.**

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.

Example 3.2**.**

[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$$2$$3$$\mathcal{D}

According to (2.1), the components of $\boldsymbol{X}$ have the following representations

[TABLE]

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

[TABLE]

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

[TABLE]

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$$2$$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 [12].

Theorem 3.3**.**

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

1. $(a)$

$\mathcal{D}^{B}\subseteq\mathcal{D}^{*}$; 2. $(b)$

$\mathcal{D}^{*}$* and $\mathcal{D}^{B}$ have the same reachability matrix;* 3. $(c)$

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

$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 [11].

The ratios $\boldsymbol{Y}=\{Y_{ij}=X_{j}/X_{i},\;i,j=1\ldots,d\}$ between all pairs of components of $\mathbf{X}$ are the essential quantities used to identify $B$ from ${\mathcal{L}}(\boldsymbol{X})$. We first present distributional properties of these ratios, where we let $(\Omega,\mathcal{F},{\mathbb{P}})$ denote the probability space of $\boldsymbol{Z}$ and, hence, of $\boldsymbol{X}$. In what follows, we use the standard convention and write events such as $\{\omega\in\Omega:X_{i}(\omega)<X_{j}(\omega)\}$ as $\{X_{i}<X_{j}\}$, etc. Unsurprisingly, because of the max-linear representation (2.3) of the components of $\boldsymbol{X}$, the ratios inherit their distributional properties from the innovations. It plays an important role that

[TABLE]

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

Lemma 3.4**.**

Let $i,j\in V$ be distinct.

1. (a)

The ratio $Y_{ji}=X_{i}/X_{j}$ has an atom in $x\in\mathbb{R}_{+}$ if and only if ${\rm{An}}(i)\cap{\rm{An}}(j)\neq\emptyset$ and $x={b_{\ell i}}/{b_{\ell j}}$ for some $\ell\in{\rm{An}}(i)\cap{\rm{An}}(j)$. 2. (b)

We have

[TABLE]

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

Proof.

To establish (a) note that (2.3) and (3.1) imply that the sets \{X_{i}=xX_{j}\}=\big{\{}\bigvee_{\ell\in{\rm{An}}(i)}b_{\ell i}Z_{\ell}=\bigvee_{\ell\in{\rm{An}}(j)}xb_{\ell j}Z_{\ell}\big{\}} and

[TABLE]

differ only by a set of probability zero. Since the innovations are independent and have support ${\mathbb{R}}_{+}$ the conclusion follows.

To establish (b) note that the support ${\mathbb{R}}_{+}$ of the innovations and the representation (2.3) yield

[TABLE]

The continuity of the function

[TABLE]

implies that ${\rm{supp}}(Y_{ji})$ is an interval in ${\mathbb{R}}_{+}$. Since for $j\in\text{an}(i)$ by Corollary 3.13 of [12] $b_{ji}\leq Y_{ji}$ and by (a) $b_{ji}$ is an atom of $Y_{ji}$, it suffices to show that $j\in\text{an}(i)$ if ${\rm{supp}}(Y_{ji})$ has a positive lower bound. For this assume that $j\not\in{\rm{an}}(i)$. Because of the positive lower bound of ${\rm{supp}}(Y_{ji})$, there exists some $a\in{\mathbb{R}}_{+}$ such that

[TABLE]

for all $\boldsymbol{z}_{{\rm{An}}(i)\cup{\rm{An}}(j)}\in{\mathbb{R}}_{+}^{|{\rm{An}}(i)\cup{\rm{An}}(j)|}$. As ${\rm{An}}(j)\setminus{\rm{An}}(i)\neq\emptyset$, for fixed $\mathbf{z}_{\text{An}(i)}\in{\mathbb{R}}_{+}^{|{\rm{An}}(i)|}$, we can choose $z_{\ell}$ for some $\ell\in{\rm{An}}(j)\setminus{\rm{An}}(i)$ so large that $ab_{\ell j}z_{\ell}$ is greater than the maximum on the right-hand side of (3.2). This contradicts (3.2). Hence, $j\in{\rm{an}}(i)$. ∎

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})$.

Algorithm 3.5**.**

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

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

For all $i,j\in V$ with $i\neq j$, find ${\rm{supp}}(Y_{ji})$: 3. if ${\rm{supp}}(Y_{ji})=[a,\infty)$ for some $a\in{\mathbb{R}}_{+}$, then set $b_{ji}=a$;

1. else, set $b_{ji}=0$.

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})$.

Algorithm 3.6**.**

[Find $F_{Z_{1}}(x),\ldots,F_{Z_{d}}(x)$ for $x\in{\mathbb{R}}_{+}$ from $B$ and ${\mathcal{L}}(\boldsymbol{X})$]

For $\nu=0,\ldots,d-1$,

       for $i\in V$ such that $|{\rm{an}}(i)|=|\{j\in V\setminus\{i\}:b_{ji}\neq 0\}|=\nu$, set

[TABLE]

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).

4 Estimation with known directed acyclic graph

In this section we consider independent realizations \boldsymbol{x}^{(t)}=\big{(}x^{(t)}_{1},\ldots,x^{(t)}_{d}\big{)}, $t=1,\ldots,n$, of a random vector $\boldsymbol{X}=(X_{1},\ldots,X_{d})$ following a recursive ML model with its DAG $\mathcal{D}$ given. Further, we consider the distribution of the innovation vector to be fixed; however, we emphasize that our estimates and their validity do not depend on this distribution as long as it prescribes independent, atom-free margins with support ${\mathbb{R}}_{+}$. Our aim is the estimation of the edge weights $c_{ki}$ and the ML coefficient matrix $B$. We recall from Theorem 3.3 that only the ML coefficient matrix $B$ can be directly identified from ${\mathcal{L}}(\boldsymbol{X})$ and hence our focus will be on the estimation of $B$; subsequently all DAGs and systems of edge weights compatible with $B$ can be obtained from Theorem 3.3.

The ML coefficient matrix $\boldsymbol{B}$

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

[TABLE]

see Theorem 4.2 or Corollary 4.3(a) of [12].

Example 4.1**.**

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

1$$2$$3$$4

and a potential ML coefficient matrix $B$ with reduction $B_{0}$, as given below.

[TABLE]

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

[TABLE]

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:

[TABLE]

Davis and Resnick [6] 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:

Example 4.2**.**

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

Consider the DAG

1$$2$$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.

Lemma 4.3**.**

Let $B_{0}\in\overline{\mathbb{R}}_{+}^{d\times d}$ be a matrix with $b_{ji}>0\iff j\to i$. A matrix $A\in\overline{\mathbb{R}}_{+}^{d\times d}$ satisfies

[TABLE]

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

Proof.

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 ([4], Proposition 1.6.10) that

[TABLE]

and hence

[TABLE]

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)$,

[TABLE]

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:

[TABLE]

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

[TABLE]

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.

Example 4.4**.**

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

Consider the DAG

1$$2$$3$$4$$\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

[TABLE]

then $\widehat{B}=B$ so we estimate all ML coefficients exactly. Note that this event has positive probability and occurs ${\mathbb{P}}$-almost surely if and only if $Z_{1}$ realizes all node variables; i.e., if $X_{2}=b_{12}Z_{1}$, $X_{3}=b_{13}Z_{1}$, and $X_{4}=b_{14}Z_{1}$.   $\Box$

Since by Table 3.1 ${\mathbb{P}}(X_{i}=b_{ki}X_{k})>0$ for $k\in{\rm{pa}}(i)$, it follows from the Borel-Cantelli lemma that $\widehat{b}_{ki}$ ${\mathbb{P}}$-almost surely equals the true value for $n$ sufficiently large. Thus, if $n$ is large, $\widehat{B}$ finds, with probability 1, the true $B$. In [6] this is discussed in a time-series framework used there and in Davis and McCormick [5] they show that under suitable assumptions in the time-series framework, this estimator is asymptotically Fréchet distributed. Assuming the probability of $\{X_{i}=b_{ki}X_{k}\}$ is known, we show next how one has to choose $n$ to observe this event with probability greater than $1-p$ for some $p\in(0,1)$. We also prove that the probability for estimating the true $b_{ki}$ converges exponentially fast to 1.

Proposition 4.5**.**

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

[TABLE]

if and only if

[TABLE]

Furthermore, the convergence {\mathbb{P}}\big{(}\bigwedge_{t=1}^{n}Y_{ki}^{(t)}=b_{ki}\big{)}\to 1 as $n\to\infty$ is exponentially fast.

Proof.

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

[TABLE]

Altogether, the statements follow. ∎

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 $\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 [11]) 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 [17], Kiefer and Wolfowitz [18], and Scholz [26] suggested such extensions. We essentially follow the Kiefer–Wolfowitz definition of a GMLE as also done, for example, by Gill et al. [10] and Johansen [16]. 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 [18].

Let ${\mathcal{P}}$ be a family of probability measures on $({\mathbb{R}}_{+}^{d},{\mathbb{B}}({\mathbb{R}}_{+}^{d}))$ where ${\mathbb{B}}({\mathbb{R}}_{+}^{d})$ denotes the Borel $\sigma$-algebra on ${\mathbb{R}}_{+}^{d}$, and $\boldsymbol{x}^{(1)},\ldots,\boldsymbol{x}^{(n)}$ a random sample from some $P_{0}\in{\mathcal{P}}$. For $P,Q\in{\mathcal{P}}$ and $\boldsymbol{x}\in{\mathbb{R}}_{+}^{d}$ we define

[TABLE]

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

[TABLE]

Since $P$ is absolutely continuous with respect to $P+Q$, the density ${dP}/{d(P+Q)}$ always exists according to the Radon-Nikodym theorem. This means that the GMLE is well-defined, save for the usual ambiguity in the method of maximum likelihood that densities are only defined up to null sets and therefore a specific choice of densities must be made. The Kiefer–Wolfowitz definition extends the definition of a MLE in a very natural way as it simply says that for any $Q\in{\mathcal{P}}$, $\widehat{P}$ is the MLE in the smaller family $\{\hat{P},Q\}$, consisting of only two distributions. In [18] only the second condition in (4.4) is required, but the first condition is implicit. The first step in verifying that $\widehat{B}$ is a GMLE of $B$ is to specify densities of $P_{B}$ with respect to $P_{B}+P_{B^{*}}$ for any two $B,B^{*}\in{\mathcal{B}}(\mathcal{D})$. For this purpose we determine a partition \big{\{}A_{0}(B,B^{*}),A_{1/2}(B,B^{*}), A_{1}(B,B^{*})\big{\}} of ${\mathbb{R}}_{+}^{d}$ that satisfies the following three properties,

[TABLE]

Then we choose as density the measurable function from ${\mathbb{R}}_{+}^{d}$ to $\{0,1/2,1\}$ defined as

[TABLE]

This is a valid density because, using the properties (A), (B), (C), we obtain for every $A\in{\mathbb{B}}({\mathbb{R}}_{+}^{d})$,

[TABLE]

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.

Example 4.6**.**

[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

[TABLE]

of ${\mathbb{R}}_{+}^{2}$ satisfies properties (A), (B), (C) of (4.5). Figure 1 shows the corresponding density $\rho(\cdot,B,B^{*})$ from (4.6) for the three possible order relations between $b_{12}$ and $b^{*}_{12}$.

Since by Table 3.1, ${\rm{supp}}({X_{2}}/{X_{1}})=[b_{12},\infty)$ and $b_{12}$ is the only atom of ${X_{2}}/{X_{1}}$, property (A) is true. By reversing the roles of $B$ and $B^{*}$, (C) follows from (A). The condition (B) is obvious if $b_{12}=b^{*}_{12}$. Assume that $b_{12}\neq b_{12}^{*}$. We then have by definition of $\boldsymbol{X}$ that $\{\boldsymbol{X}\in A_{1/2}(B,B^{*})\}=\{X_{2}>(b_{12}\vee b^{*}_{12})X_{1}\}=\{Z_{2}>(b_{12}\vee b^{*}_{12})Z_{1}\}$ and $X_{2}=Z_{2}$ on $\{Z_{2}>(b_{12}\vee b^{*}_{12})Z_{1}\}$. With this, using that $A_{1/2}(B^{*},B)=A_{1/2}(B,B^{*})$, we obtain for $A\in\mathbb{B}({\mathbb{R}}_{+}^{2})$,

[TABLE]

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

[TABLE]

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

[TABLE]

and the second condition if and only if

[TABLE]

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.

1. (a)

$\widetilde{b}_{12}<\widehat{b}_{12}$ is no GMLE:

Set $b_{12}=\widehat{b}_{12}$, and let $\boldsymbol{x}^{(t)}$ be such that $\widehat{b}_{12}x_{1}^{(t)}=x_{2}^{(t)}$. Then \boldsymbol{x}^{(t)}\in\big{\{}\boldsymbol{x}\in{\mathbb{R}}_{+}^{2}:x_{2}=b_{12}x_{1}>\widetilde{b}_{12}x_{2}\big{\}}\subseteq A_{0}(\widetilde{B},B) but no \boldsymbol{x}^{(s)}\in A_{0}(B,\widetilde{B})=\big{\{}\boldsymbol{x}\in{\mathbb{R}}_{+}^{2}:x_{2}<b_{12}x_{1}\big{\}}. This contradicts (4.8); consequently, $\widetilde{b}_{12}$ cannot be a GMLE of $b_{12}$. In Figure 2(a) we illustrate this situation. On the left-hand side a contour plot of the density $\rho(\cdot,\widetilde{B},B)$ is shown, on the right-hand side of $\rho(\cdot,B,\widetilde{B})$. The crosses represent the realizations $\boldsymbol{x}^{(1)},\ldots,\boldsymbol{x}^{(n)}$. In the left plot crosses are in the [math]-area coloured in red, namely, those that realize $\widehat{b}_{12}$, but in the right plot not. So $\widetilde{B}$ cannot be a GMLE of $B$. 2. (b)

$\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). 3. (c)

$\widetilde{b}_{12}=\widehat{b}_{12}$ is a GMLE:

Condition (4.7) holds obviously. To prove (4.8), assume for some $B\in{\mathcal{B}}(\mathcal{D})$ that some $\boldsymbol{x}^{(t)}\in A_{0}(\widehat{B},B)$. By definition of $A_{0}(\widehat{B},B)$, $x_{2}^{(t)}=b_{12}x_{1}^{(t)}>\widehat{b}_{12}x_{1}^{(t)}$, which implies that $b_{12}>\widehat{b}_{12}$. For $\boldsymbol{x}^{(s)}$ such that $\widehat{b}_{12}x_{1}^{(s)}=x_{2}^{(s)}$, we then find that $x_{2}^{(s)}<b_{12}x_{1}^{(s)}$. Hence, $\boldsymbol{x}^{(s)}\in A_{0}(B,\widehat{B})$, and $\widehat{b}_{12}$ is a GMLE of $b_{12}$. We learn this informally from Figure 2(c). The top line shows contour plots of $\rho(\cdot,\widehat{B},B)$ for the three different orders between $b_{12}$ and $\widehat{b}_{12}$, and the bottom line shows the corresponding contour plots of $\rho(\cdot,B,\widehat{B})$. The two plots on the left-hand side correspond to the situation from above: in the upper plot there are realizations in the [math]-area, namely those that are on the line \big{\{}\boldsymbol{x}\in{\mathbb{R}}_{+}^{2}:x_{2}=b_{12}x_{1}\big{\}}, but then there are also realizations in the [math]-area of the lower plot (those that lie below this line). Hence, (4.8) holds. Since there is no realization in the [math]-area of the middle and right plot in the top line, (4.8) is automatically satisfied if $b_{12}\leq\widehat{b}_{12}$.

$\Box$

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$.

Our partition \big{\{}A_{0}(B,B^{*}),A_{1/2}(B,B^{*}),A_{1}(B,B^{*})\big{\}} of ${\mathbb{R}}_{+}^{d}$ is based on the following representation for the components of $\boldsymbol{X}$:

[TABLE]

We begin with the specification of $A_{1/2}(B,B^{*})$ and prove a property needed subsequently to verify property (B). Have in mind that if $b_{ki}>b^{*}_{ki}$ for all $k\in{\rm{pa}}(i)$ or $b_{ki}<b^{*}_{ki}$ for all $k\in{\rm{pa}}(i)$ then \big{\{}\boldsymbol{x}\in{\mathbb{R}}_{+}^{d}:x_{i}=\bigvee_{k\in{\rm{pa}}(i)}b^{*}_{ki}x_{k}=\bigvee_{k\in{\rm{pa}}(i)}b_{ki}x_{k}\big{\}}=\emptyset.

Lemma 4.7**.**

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

[TABLE]

Then for every $F\in{\mathcal{F}}$,

[TABLE]

Proof.

The proof is deferred to the appendix. ∎

As a partition of ${\mathbb{R}}_{+}^{d}$ we now suggest \big{\{}A_{0}(B,B^{*}),A_{1/2}(B,B^{*}),A_{1}(B,B^{*})\big{\}}, where $A_{1/2}(B,B^{*})$ is defined above in Lemma 4.7,

[TABLE]

and A_{1}(B,B^{*})={\mathbb{R}}_{+}^{d}\setminus\big{(}A_{0}(B,B^{*})\cup A_{1/2}(B,B^{*})\big{)}. With this partition we then have:

Theorem 4.8**.**

Let $B,B^{*}\in{\mathcal{B}}(\mathcal{D})$. Then the function $\rho:{\mathbb{R}}_{+}^{d}\to\{0,1/2,1\}$

[TABLE]

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

Proof.

See the appendix. ∎

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

Example 4.9**.**

[Local densities $\rho_{i}$]

Consider the DAGs

[TABLE]

Let $\rho$, $\rho_{2}$, and $\rho_{3}$ be the corresponding densities from (4.11). For the ML coefficient matrix $B$ of a recursive ML model on $\mathcal{D}$, let $B_{2}$ and $B_{3}$ be the ML coefficient matrices of recursive ML models on $\mathcal{D}_{2}$ and $\mathcal{D}_{3}$ with edge weight $c_{12}=b_{12}$ and $c_{23}=b_{23}$, and let starred quantities denote the same for $B^{*}$. We then find for $\boldsymbol{x}=(x_{1},x_{2},x_{3})\in{\mathbb{R}}_{+}^{3}$,

[TABLE]

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:

[TABLE]

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.

Proposition 4.10**.**

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.

1. (a)

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

[TABLE] 2. (b)

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

[TABLE]

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

The function $\rho_{i}:{\mathbb{R}}_{+}^{d}\to\{0,1/2,1\}$ such that $\boldsymbol{x}_{{\rm{Pa}}(i)}\mapsto\rho_{i}(\boldsymbol{x}_{{\rm{Pa}}(i)},B_{i},B_{i}^{*})$ is a density of $P^{i\mid{\rm{pa}}(i)}_{B}$ with respect to $P^{i\mid{\rm{pa}}(i)}_{B}+P^{i\mid{\rm{pa}}(i)}_{B^{*}}$, where $P_{B}^{i\mid{\rm{pa}}(i)}$ is a regular conditional distribution of $X_{i}$ given $\boldsymbol{X}_{{\rm{pa}}(i)}$ and $P_{B^{*}}^{i\mid{\rm{pa}}(i)}$ one of $X^{*}_{i}$ given $\boldsymbol{X}^{*}_{{\rm{pa}}(i)}$.

Proof.

See the appendix. ∎

Next, we show that $\widehat{B}$ is indeed a GMLE in the sense of [18]. 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.

Theorem 4.11**.**

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.

1. (a)

The matrix $\widehat{B}$ from (4.3) is a GMLE of $B$. 2. (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}$. 3. (c)

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}$.

Proof.

(a) First, recall that $\widehat{B}$ is indeed a ML coefficient matrix of a recursive ML model on $\mathcal{D}$. The first condition in the definition of a GMLE in (4.4) is satisfied due to the definition of $\rho(\cdot,\widehat{B},\widehat{B})$ since $A_{1/2}(\widehat{B},\widehat{B})=\mathbb{R}_{+}^{d}$. Since the densities $\rho(\cdot,\widehat{B},B)$ and $\rho(\cdot,B,\widehat{B})$ have the values [math], $1$, $1/2$, and $A_{1/2}(\widehat{B},B)=A_{1/2}(B,\widehat{B})$, to verify the second condition in (4.4), it suffices to show that there is some realization $\boldsymbol{x}^{(t_{1})}\in A_{0}(B,\widehat{B})$ whenever there is some realization $\boldsymbol{x}^{(t_{2})}\in A_{0}(\widehat{B},B)$; cf. Example 4.6, in particular (4.8). So let $\boldsymbol{x}^{(t_{2})}\in A_{0}(\widehat{B},B)$ for some $t_{2}\in\{1,\ldots,n\}$. We find, for some $i\in V$, from the definition of $A_{0}(\widehat{B},B)$ and the fact that $x_{i}^{(t)}\geq\bigvee_{k\in{\rm{pa}}(i)}\widehat{b}_{ki}x_{k}^{(t)}$,

[TABLE]

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.

5 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.

Algorithm 5.1**.**

[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$. 2. 2.

For all $i,j\in V$ with $i\neq j$ ,

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)}$;

2. else,    else, set $\widecheck{b}_{ji}=0$.

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 [12]) 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.

6 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.

We emphasize again that, although our estimates are derived under the assumption that the distribution of the innovation vector $\boldsymbol{Z}$ is fixed, the estimates do not depend on what this distribution is and would therefore also be valid in the situation where the innovations are independent with unkown distributions that are atom-free and have support equal to ${\mathbb{R}}_{+}$. Algorithm 3.6 provides a recursive procedure to obtain the distribution functions $F_{Z_{i}}$ from $B$ and the marginal distribution functions $F_{X_{i}}$ of $X_{i}$. Estimating $B$ by $\widehat{B}$ and the distributions $F_{X_{i}}$, for example, by their empirical versions, we can apply this procedure to find estimators of the distributions $F_{Z_{i}}$ although it will formally violate the assumption of atom-freeness and thus it is both more efficient and formally correct to estimate these parametrically, or under suitable monotonicity restrictions.

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 [5].

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.

Appendix A Appendix: some technical proofs

Proof of Lemma 4.7

Proof.

First, define for $i\in V$

[TABLE]

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

[TABLE]

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

[TABLE]

Noting from the proof of Theorem 4.2 of [12] that

[TABLE]

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

[TABLE]

Thus, again by (2.3),

[TABLE]

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

Proof.

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

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

[TABLE]

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

[TABLE]

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

(B) Recall that $P_{B}$ and $P_{B^{*}}$ share the same innovation vector when represented by a recursive ML model. Furthermore, note that the set $\Omega(B,B^{*})$ from Lemma 4.7 is a subset of \bigcap_{i\in V}\big{\{}X_{i}=\bigvee_{j\in{\rm{An}}(i):b_{ji}=b^{*}_{ji}}b_{ji}Z_{j}\big{\}}. We have $\Omega(B,B^{*})=\Omega(B^{*},B)$ and hence we obtain from (4.10) for $A\in{\mathbb{B}}({\mathbb{R}}_{+}^{d})$,

[TABLE]

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

[TABLE]

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

Proof.

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$,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

is a regular conditional distribution function of $X^{*}_{i}$ given $\boldsymbol{X}^{*}_{{\rm{pa}}(i)}$. Figure 5 depicts the two conditional distribution functions for the three possible orders between $\bigvee_{k\in{\rm{pa}}(i)}b_{ki}x_{k}$ and $\bigvee_{k\in{\rm{pa}}(i)}b^{*}_{ki}x_{k}$. It then suffices to show for all $\boldsymbol{x}_{{\rm{pa}}(i)}\in{\mathbb{R}}_{+}^{|{\rm{pa}}(i)|}$ and $y\in{\mathbb{R}}_{+}$,

[TABLE]

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

[TABLE]

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

∎