Identifying Nonlinear 1-Step Causal Influences in Presence of Latent Variables

TL;DR

This paper introduces an information-theoretic method to identify 1-step causal influences in stochastic dynamical systems with latent variables, including a linear regression-based approach for linear dynamics, validated through simulations.

Contribution

It presents a novel approach for causal discovery in systems with latent variables, extending existing methods to nonlinear and linear cases with validation.

Findings

Successfully recovers causal relations among observed variables.

Effective in systems where latent variables lack exogenous noise.

Validated through numerical simulations demonstrating practical applicability.

Abstract

We propose an approach for learning the causal structure in stochastic dynamical systems with a $1$-step functional dependency in the presence of latent variables. We propose an information-theoretic approach that allows us to recover the causal relations among the observed variables as long as the latent variables evolve without exogenous noise. We further propose an efficient learning method based on linear regression for the special sub-case when the dynamics are restricted to be linear. We validate the performance of our approach via numerical simulations.

Equations

$$\exists h>0,\mbox{ $s.t.$ }\Sigma_{X_{1}}(h|\Omega_{t})<\Sigma_{X_{1}}(h|\Omega_{t}\backslash\{X_{2}(s)\}_{s=0}^{t}),$$

$$\left[\begin{array}{c}{X}_1(t)\\ X_2(t)\\ Z(t)\end{array}\right]=\left[\begin{array}{ccc}0& 0& 0.5\\ 0.5& 0.1& 0.9\\ 0.9& 0& 0.5\end{array}\right]\left[\begin{array}{c}{X}_1(t-1)\\ X_2(t-1)\\ Z(t-1)\end{array}\right]+\left[\begin{array}{c}{\omega }_1(t)\\ {\omega }_2(t)\\ 0\end{array}\right],$$

$$\displaystyle\left\{\begin{array}[]{ll}X_{i}(t)=F_{i}(\vec{X}(t-1),\vec{Z}(t-1))+\omega_{i}(t),&1\leq i\leq n,\\ \vec{Z}(t)=G(\vec{X}(t-1),\vec{Z}(t-1)),&\end{array}\right.$$

$${\alpha }_{i,j}(t):=\underset{\begin{array}{c}{v}_j,v_j^{\prime },{\vec{v}}_{-j}\end{array}}{\sup }|F_i(v_j,{\vec{v}}_{-j})-F_i(v_j^{\prime },{\vec{v}}_{-j})|,$$

$$X_1(t)=e^{-|X_1(t-1)+X_2(t-1)|}-e^{-|Z(t-1)|}/5+{\omega }_1(t),$$

$$X_2(t)=\sqrt{|X_2(t-1)|/2}+{\omega }_2(t),$$

$$Z(t)=X_2(t-1)Z(t-1),$$

$$I(V_j\to V_i||{\underline{V}}_{\mathcal{K}}):={\mathbb{E}}_{P_{{\underline{V}}_{\mathcal{K}\cup \{i,j\}}}}\left[\log \frac{d{P}_{V_i||V_j,{\underline{V}}_{\mathcal{K}}}}{d{P}_{V_i||{\underline{V}}_{\mathcal{K}}}}\right],$$

$$I(V_j\to V_i||{\underline{V}}_{-\{i,j\}})>0.$$

$$X_1(t)=X_1(t-1)/3+W_1(t),$$

$$Z(t)=X_1(t-1),$$

$$X_2(t)=Z(t-1)/3+W_2(t),$$

$$I(X_i(t);X_j(t-1)|{\underline{X}}_0^{t-1}\backslash \{X_j(t-1)\})=0.$$

$$H(X_i$$

$$=H(F_i(\vec{X}(t-1),{\Psi }_{t-1}({\underline{X}}_0^{t-2},\vec{Z}(0)))+{\omega }_i(t)|{\underline{X}}_0^{t-1}),$$

$$=H({\omega }_i(t)|{\underline{X}}_0^{t-1})=H({\omega }_i(t)).$$

$$\begin{array}{rl}H(X_i(t)& |{\underline{X}}_0^{t-1}\backslash \{X_j(t-1)\})\\ & =H(F_i^{\prime }(X_j(t-1))+{\omega }_i(t)|{\underline{X}}_0^{t-1}\backslash \{X_j(t-1)\}),\end{array}$$

$$\mathbb{E}\{X_i(t)|{\underline{X}}_0^{t-1}\}=\mathbb{E}\{X_i(t)|{\underline{X}}_0^{t-1}\backslash \{X_j(t-1)\}\}$$

$$\mathbb{E}\{X_i(t)|{\underline{X}}_0^{t-1}={\underline{x}}_0^{t-1}\}=F_i(\vec{x}(t-1),\vec{z}(t-1)).$$

$$\left[\begin{array}{c}\vec{X}(t)\\ \vec{Z}(t)\end{array}\right]=\left[\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right]\left[\begin{array}{c}\vec{X}(t-1)\\ \vec{Z}(t-1)\end{array}\right]+\left[\begin{array}{c}\vec{\omega }(t)\\ {\vec{0}}_{m\times 1}\end{array}\right],$$

$$\begin{array}{rl}\vec{X}(t)& =\sum _{k=0}^{t-1}{A}_k^{\ast }\vec{X}(t-1-k)+A_{12}{A}_{22}^{t-1}\vec{Z}(0)+\vec{\omega }(t),\end{array}$$

$$\begin{array}{rl}& \mathbb{E}\{X_i(t)|{\underline{x}}_0^{t-1}\}-\mathbb{E}\{X_i(t)|{\underline{x}}_0^{t-1}\backslash \{x_j(t-1)\}\}=\\ & =A_{11}(i,j)\left(x_j(t-1)-\mathbb{E}\{X_j(t-1)|{\underline{x}}_0^{t-1}\backslash \{x_j(t-1)\}\}\right)\end{array}$$

$$\begin{array}{rl}\vec{X}(t)=& \sum _{k=0}^l{A}_k^{\ast }\vec{X}(t-1-k)+\vec{\omega }(t),t\ge l+1.\end{array}$$

$$V_i(t)=w_{ii}{V}_i(t-1)+\sum _{j\ne i}{w}_{ij}(V_j(t-1)-V_i(t-1))+B_i{\omega }_i(t),$$

$$A(i,k)=\begin{cases}w_{ik},&\mbox{if $i\neq k$},\\ w_{ii}-\sum_{j\neq i}w_{ij},&\mbox{otherwise,}\end{cases}$$

$$\begin{array}{rl}I(X_i(t);X_j(t-1)& |{\underline{X}}_0^{t-1}\backslash X_j(t-1))=\\ & =\frac{1}{2}\log \left(1+\frac{(A_{11}(i,j){)}^2{\sigma }_{{\omega }_j}^2}{{\sigma }_{{\omega }_i}^2}\right),t>l\end{array}$$

$$\displaystyle\left\{\begin{array}[]{lll}X_{1}(t)&=0.2X_{1}(t-1)+0.4\sqrt{|Z(t-1)|}+\omega_{1}(t),\\ X_{2}(t)&=0.5(X_{1}(t-1))^{2}+0.9Z(t-1)+\omega_{2}(t),\\ Z(t)&=0.9(X_{1}(t-1))^{3}+0.4Z(t-1),\end{array}\right.$$

Full text

1 Introduction

Identifying causal influences in a network of time series is one of fundamental problems in many different fields, including social sciences, economics, computer science, and biology. In macroeconomics, for instance, researchers seek to understand what are the factors contributing to economic fluctuations and how these factors interact with each other [12]. In neuroscience, extensive body of research focuses on learning the interactions between different regions of brain by analyzing neural spike trains [16].

In 1960’s, Granger proposed a definition of causality between random processes [8]. The key idea of his definition is that if a process $X_{2}(t)$ causes another process $X_{1}(t)$, then knowing the past of $X_{2}(t)$ up to time $t$ must aid in predicting $X_{1}(t)$. In particular, let $\Sigma_{X_{1}}(h|\Omega_{t})$ be the mean square error (MSE) of the optimal $h$-step predictor of a random process $X_{1}(t)$ at time $t$ given information $\Omega_{t}$. Process $X_{2}(t)$ is said to Granger cause process $X_{1}(t)$ if:

[TABLE]

where the set $\Omega_{t}$ contains all information in the universe related to the past and the present of $X_{1}(t)$. We also say that the process $X_{2}(t)$ has a 1-step cause on $X_{1}(t)$ if the above inequality holds for $h=1$. In other words, considering $X_{2}(t)$ in the set $\Omega_{t}$ improves prediction of $X_{1}(t+1)$.

Granger’s definition of causality is consistent with the belief that a cause cannot come after the effect, but it is not practical in some settings because it requires knowledge of the entire set $\Omega_{t}$. To put it differently, it is hard to identify and account for all parts of universe that are related to a specific process $X_{1}(t)$. Hence, only the available information related to $X_{1}(t)$ is considered in practice [11]. To see what may go wrong in such a situation, consider the following linear model with three state variables:

[TABLE]

where $X_{1}(t)$ and $X_{2}(t)$ are observable but $Z(t)$ is latent. Let $\omega_{1}(t)$ and $\omega_{2}(t)$ be i.i.d random variables with the same variance. If we fit a linear model only on $X_{1}(t)$ and $X_{2}(t)$ without considering $Z(t)$, our estimation of upper left $2\times 2$ submatrix would be: $[0.06,0.32;0.61,0.69]$. This result implies that $X_{2}(t)$ is a 1-step cause of $X_{1}(t)$ with the strength $0.32$ which is wrong. The concept of Granger causality can be generalized to nonlinear setting using an information-theoretic quantity “directed information” [13]. Still the problem caused by latent processes persists in that setting as well.

Identifying causal relations between random variables has been studied in the presence of latent variables to some extent. For instance, Elidan et al. proposed an algorithm based on expectation maximization (EM) to estimate the parameters of their model by fixing the number of latent variables and also the structural relationships between latent and observed variables [3]. Chandrasekaran et al. [1] presented a tractable convex program based on regularized maximum likelihood for recovering causal relations for a model where the latent and observed variables are jointly Gaussian, and the conditional statistics of the observed variables given the latent variables is a sparse graph. A well-known approach for learning latent Markov models uses quartet-based distances to discover the structure [10, 4]. In most of quartet-based solutions, a set of quartets is constructed for all subsets of four observable variables and then quartets are merged to form a tree structure.

In recent years, there has been an increasing interest in inferring causal relations in random processes. Jalali and Sanghavi showed that 1-step causal relations between observed variables can be identified in a Vector Auto-Regressive (VAR) model assuming that connections between observed variables are sparse and each latent variable interacts with many observed variables [9]. In [7], Geiger et al. showed that identifying 1-step causes between observed variables is possible under some algebraic conditions on the transition matrix of VAR model. Recently, Etesami et al. studied a network of processes with polytree structure and introduced an algorithm that can learn latent polytrees using a discrepancy measure [6].

In this paper, we propose an information-theoretic criteria for identifying the causal relations in a general model of stochastic dynamical systems without restricting the mapping functions (say to linear mappings) or the underlying structure (e.g., being a tree) among the observed nodes also when there is no exogenous noise in the latent part. We propose an efficient method to identify functional dependencies for the special case of linear mappings. We further demonstrate the applicability of this criteria though simulation results for both linear and nonlinear cases.

The paper is organized as follows. In Section 2, we provide the preliminary definitions and describe the system model. In Section III, we present the main result and study the special restriction of it to linear models. We provide our simulation results in Section IV. Finally, we conclude in Section V.

2 Problem Definition

In this section, after some notational conventions, the model of stochastic dynamical system is presented. Afterwards, we present our definition of 1-step functional dependency between the processes for this model.

2.1 Notations

Any $n\times 1$ vector with with entries $[V_{1}(t);\cdots;V_{n}(t)]$ is denoted by $\vec{V}(t)$. We denote the $t$-th random variable in the $i$-th process by $V_{i}(t)$. We use underlined characters to represent a collection of processes, for example $\underline{V}_{\mathcal{K},0}^{t}$ is used to denote a set of random processes with index set $\mathcal{K}$ from time [math] up to time $t$. For $\mathcal{K}=\{1,...,n\}$, we denote $\underline{V}_{\mathcal{K},0}^{t}$ by $\underline{V}_{0}^{t}$. We also define: $-j:=\{1,\cdots,n\}\backslash\{j\}$. The identity matrix of size $n$ is shown by $I_{n\times n}$. We denote $(i,j)$ entry of a matrix $A$ by $A(i,j)$.

In a directed graph $\overrightarrow{G}=(V,\overrightarrow{E})$ that is characterized by a set $V$ of vertices (or nodes) and a set of ordered pairs of vertices, called arrows (or edges) $\overrightarrow{E}\subset V\times V$, we denote the set of parents of a node $v$ by $\mathcal{PA}(v)$ and define it as $\mathcal{PA}(v):=\{u\in V:(u,v)\in\overrightarrow{E}\}$.

2.2 System model

Consider a dynamical system described by $n+m$ states $\vec{V}(t)=[V_{1}(t),\cdots,V_{n+m}(t)]$ in which the first $n$ processes, denoted by $\vec{X}=[X_{1},...,X_{n}]$ are observable states and the rest which denoted by $\vec{Z}=[Z_{1},...,Z_{m}]$ are latent. More precisely, the joint dynamic of the system is given by:

[TABLE]

where the exogenous noises $\{\omega_{i}(t)\}$ are i.i.d. with mean zero. $F_{i}:\mathbb{R}^{n+m}\rightarrow\mathbb{R}$, $G:\mathbb{R}^{n+m}\rightarrow\mathbb{R}^{m}$ are mapping functions that belong to appropriately constrained class of functions. Furthermore, we assume that $\vec{Z}(0)$ is a vector of unknown but fixed values. The goal of this work is to identify the causal structure among the observed processes $\vec{X}$ given their realizations. Next, we formally introduce what we mean by a causal structure of a dynamical system.

2.3 Causal Structural Graph

In dynamical systems with functional dependencies, there is a natural notion of influence among the processes, in the sense that process $V_{j}$ causes process $V_{i}$, if $V_{i}$ is a function of $V_{j}$. Such dependencies has been studied in the literature [5]. Adopting the definition of functional dependency in [5], we define the causal structure of the system in (5) as follows.

Random process $V_{i}$ 1-step functionally depends on process $V_{j}$ over the time horizon $[0,T]$, if changing the value of $V_{j}(t-1)$ while keeping all the other variables fixed results in a change in $V_{i}(t)$ for some time $0<t\leq T$. Next, we present our formal definition of functional dependencies in systems whose joint dynamics is described by (5).

Definition 1**.**

We say $V_{j}$ 1-step functionally influences $V_{i}$ if and only if $\alpha_{i,j}:=1/T\sum_{t=0}^{T}\alpha_{i,j}(t)>0$, where

[TABLE]

$(v_{j},\vec{v}_{-j})$* and $(v^{\prime}_{j},\vec{v}_{-j})$ are two realizations of $(V_{j}(t-1),\vec{V}_{-j}(t-1))$.*

In order to visualize the causal structure in (5), we introduce a directed graph whose nodes represent random processes and there is an arrow from node $j$ to nodes $i$, if $V_{j}$ 1-step functionally influences $V_{i}$.

Example 1**.**

Consider a causal system with 3 processes such that their joint dynamic is given by:

[TABLE]

where $\omega_{i}$s are independent exogenous noises. Figure 1 depicts the functional dependency graph of this system.

Directed Information Graphs (DIGs) are another type of graphical models that encode statistical dependencies in dynamical systems [2]. These graphs are defined using an information-theoretic measure, the “conditional directed information” [14, 18]. The relationship between the functional dependencies in a stochastic dynamical system and their corresponding DIG has been studied in [5].

For the sake of completeness, we present the definition of DIG. Consider two random processes $V_{i}$ and $V_{j}$ and a set of indices $\mathcal{K}$ such that $\mathcal{K}\subseteq\{1,...,n\}\setminus\{i,j\}$, then the conditional directed information from $V_{j}$ to $V_{i}$, given $\underline{V}_{\mathcal{K}}$ is defined as:

[TABLE]

where $\frac{dP_{{{\underline{V}}}}}{dQ_{{{\underline{V}}}}}$ is the Radon-Nikodym derivative [19] and $P_{V_{i}||\underline{V}_{\mathcal{K}}}$ denotes the causal conditioning defined as $P_{V_{i}||\underline{V}_{\mathcal{K}}}:=\prod_{t\geq 1}P_{V_{i}(t)|V_{i}^{t-1},\underline{V}^{t-1}_{\mathcal{K}}}.$

Definition 2**.**

[17]** A directed information graph (DIG) is a directed graph, $\overrightarrow{G}_{DI}$, over a set of random processes $\underline{V}$. Node $i$ represents the random process $V_{i}$; there is an arrow from $j$ to $i$ for $i,j\in\{1,...,n\}$ if and only if:

[TABLE]

Note that in the definition of DIG, it is assumed that there are no latent processes. Thus as demonstrated in the example below, when a subset of processes is not observable (as in our model), the corresponding DIG may not encode the 1-step causal relationships accurately.

Example 2**.**

Consider the following joint dynamics:

[TABLE]

where $\{W_{1},W_{2}\}$ are independent exogenous noises. The corresponding DIG of this system when all processes are observed is $X_{1}\rightarrow Z\rightarrow X_{2}$, and when $Z$ is latent, it is $X_{1}\rightarrow X_{2}$. But we know that there is no 1-step functional dependency between $X_{1}$ and $X_{2}$.

Definition 3**.**

A joint distribution $P_{\underline{V}}$ is called positive if there exists a reference measure $\phi$ such that $P_{\underline{V}}\ll\phi$ and $\frac{dP_{{\underline{V}}}}{d\phi}>0$.

Remark 1**.**

In addition to requiring no latent processes, DIGs recover the structure correctly when underlying distribution is positive. This is to avoid degenerate cases that arise with deterministic relationships. For instance, suppose $X$ and $Y$ are two random processes such that $Y=\psi(X)$ for some deterministic function $\psi$. Then $P_{X,Y}$ is not positive since the distribution of $Y$ given $X$ is a point mass.

Note that our model in (5) does not satisfy the non-degeneracy assumption. This is because in this model the hidden variables are a deterministic functions of the other processes. Yet as we will show next, the 1-step causal structure between the observed processes is unique and recoverable as long as the marginal distribution of the observed processes is positive.

3 Main Result

Herein, we introduce our approach for learning the 1-step functional dependencies among the observed variables given their realizations. This approach does not require any prior knowledge about the number of latent process nor functions $\{F_{i}\}$ and $G$.

Theorem 1**.**

Consider the dynamical system in (5) and assume that the marginal distribution of the observed variables $P_{\underline{X}}$ is positive. Then $\alpha_{i,j}(t)=0$ if and only if:

[TABLE]

Proof.

First, we prove that if $\alpha_{i,j}(t)=0$ then (9) holds. Suppose that $X_{i}(t)$ does not 1-step functionally depend on $X_{j}(t-1)$. According to (5), the latent vector $\vec{Z}(t)$ can be determined recursively as a function of $\underline{X}_{0}^{t-1}$ and $\vec{Z}(0)$. We denote this by $\vec{Z}(t)=\Psi_{t}(\underline{X}_{0}^{t-1},\vec{Z}(0))$. Therefore, the entropy of $X_{i}(t)$ given $\underline{X}_{0}^{t-1}$ will be:

[TABLE]

The last equation holds because $F_{i}(\vec{X}(t-1),\Psi_{t-1}(\underline{X}_{0}^{t-2},\vec{Z}(0)))$ is a deterministic function of $\underline{X}_{0}^{t-1},\vec{Z}(0)$, and $\omega_{i}(t)$ is independent of $\underline{X}_{0}^{t-1}$. Furthermore, we have:

[TABLE]

where $F_{i}^{\prime}(X_{j}(t-1))$ is a uni-variate function obtained from $F_{i}(\vec{X}(t-1),\Psi_{t-1}(\underline{X}^{t-2}_{0},\vec{Z}(0)))$ by determining the values of $\underline{X}_{0}^{t-1}\backslash X_{j}(t-1)$. But $F_{i}^{\prime}(X_{j}(t-1))$ does not change by varying $X_{j}(t-1)$ since we assumed $\alpha_{i,j}(t)=0$. Hence, the above equation is equal to $H(\omega_{i}(t))$ and by comparing with (3), we can deduce that (9) holds.

For the converse, note that $X_{i}(t)$ and $X_{j}(t-1)$ are independent given $\underline{X}_{0}^{t-1}\backslash\{X_{j}(t-1)\}$ according to (9). Consequently, we have:

[TABLE]

For any realization of $\underline{X}_{0}^{t-1}$ like $\underline{x}_{0}^{t-1}$, the left hand side of the above equation is equal to:

[TABLE]

where $\vec{z}(t-1)=\Psi_{t-1}(\underline{x}_{0}^{t-2},\vec{Z}(0))$. Since the joint distribution of the observed processes $P_{X}$ is positive, we know that $X_{j}(t-1)$ cannot be written as a deterministic function of $\underline{X}_{0}^{t-1}\backslash\{X_{j}(t-1)\}$. Thus, the right hand of (12) does not depends on $X_{i}(t-1)$. From this fact and (13), we can conclude that $F_{i}(\vec{X}(t-1),\vec{Z}(t-1))$ is not a function of $X_{j}(t-1)$ for any realization of $\underline{X}_{0}^{t-1}$ and thus $\alpha_{i,j}(t)=0$. ∎

This result can be used to recover the 1-step causal structure of the observed processes in (5) given $X_{0}^{t}$. To do so, one can estimate the conditional mutual information in (9) for all $t$. If (9) holds, then we declare that there is no 1-step dependency from $X_{j}$ to $X_{i}$. Next, we propose an efficient method to learn the 1-step causal structure of the observed processes in (5) when $F_{i}$s and $G$ are linear functions.

3.1 The Linear Model

Suppose $\{F_{i}\}$s and $G$ are linear functions, then the equations in (5) can be rewritten as follows:

[TABLE]

where $[A_{11}]_{n\times n}$, $[A_{12}]_{n\times m}$, $[A_{21}]_{m\times n}$, and $[A_{22}]_{m\times m}$ denote the coefficient matrices. We also define $A=[A_{11},A_{12};A_{21},A_{22}]$. The functional dependency of state vector $\vec{X}(t)$ on its history $\underline{X}_{0}^{t-1}$ and also $\vec{\omega}(t)$, and $\vec{Z}(0)$ for $t>1$ can be written as follows:

[TABLE]

where $A_{0}^{*}=A_{11}$ and $A_{k}^{*}=A_{12}A_{22}^{k-1}A_{21},k\geq 1$. Now, suppose that information-theoretic criteria in (9) is zero. By the same arguments in the proof of Theorem 9, we can show that the following term is zero:

[TABLE]

for any realization of $\underline{X}_{0}^{t-1}=\underline{x}_{0}^{t-1}$. Since $P_{X}$ is positive, we can deduce that $A_{11}(i,j)=0$. Consequently, learning the 1-step causal structure among the observed processes reduces to determining the support of $A_{11}$.

Assume that the support of $A_{22}$ corresponds to an acyclic directed graph, i.e. there exists an $l>0$ such that $A_{22}^{l}=0$. Under this condition, the equation (15) can be simplified as:

[TABLE]

The above equation can be interpreted as a VAR model of order $l+1$. Hence, all matrices $\{A_{k}^{*}\}$ can be obtained by doing multivariate least square estimation [11]. Moreover, coefficients in the VAR model can be checked for zero constraints by Wald test [11]. Thus, we can check the information-theoretic criteria merely performing a Wald test.

4 Experimental Results

In this section, we utilize the method described in previous part for network identification problem in consensus protocols [15]. In control systems, a well-known approach for network identification is based on running a series of “node-knockout” experiments in which variables are sequentially forced to be zero without being removed from the network [15, 20]. The main drawback of this approach is that we need to intervene in the system. Here, we will show that the direct edges between observed nodes can be detected just by analyzing the time-series of observed processes.

Consider the weighted consensus protocol within a system with $n+m$ nodes:

[TABLE]

where $V_{i}(t)$ represents the state of node $i$ at time $t$ such as its speed, heading, or position, and the weight $w_{ij}$ denotes the weight on the edge $(i,j)$. The first $n$ state variables correspond to states of observed nodes and the rest belong to hidden nodes. We are trying to find all directed edges (with nonzero weight) between observed nodes by injecting the white noise $\omega_{i}(t)$ into the observed node $i$, i.e. $B_{i}=1$ if $V_{i}$ is an observed node and $B_{i}=0$ otherwise. In fact, this problem can be reformulated to the form in (14) such that:

[TABLE]

where $[V_{1}(t),\cdots,V_{n}(t)]=[X_{1}(t),\cdots,X_{n}(t)]$, and $[V_{n+1}(t),\cdots,V_{n+m}(t)]=[Z_{1}(t),\cdots,Z_{m}(t)]$. Hence, identifying all directed edges with nonzero weight between observed nodes is equivalent to obtaining the support of matrix $A_{11}$.

We generated 1000 instances of the linear system with $n=10$ observed nodes and $m=10$ latent nodes. The weight $w_{ij}$ $(i\neq j)$ was selected randomly from the set $\{-b,0,b\}$ with probability $[q,1-2q,q]$ where $q=0.1$ and $b=0.7$ if $i,j$ were hidden. Otherwise, the weight $w_{ij}$ $(i\neq j)$ was chosen randomly from the set $\{-a,0,a\}$ with probability $[p,1-2p,p]$ where $a=0.2$. Moreover, we set $w_{ii}$ to $\sum_{(i,j)\in\vec{E}}w_{ij}$.

In our simulations, we excluded the generated networks which had cycles in the latent part. Furthermore, the noise process $\vec{\omega}(t)$ was chosen as i.i.d $\mathcal{N}(\vec{0}_{n\times 1},\sigma^{2}I_{n\times n})$ with $\sigma^{2}=0.1$. It can be easily seen that the conditional mutual information in (5) is equal to:

[TABLE]

where $\sigma_{\omega_{i}}^{2}$ and $\sigma_{\omega_{j}}^{2}$ are the variances of $\omega_{i}(t)$ and $\omega_{j}(t)$, respectively. Thus, learning $1$-step functional dependencies, corresponds to finding the support of matrix $A_{11}$, denoted by $Supp(A_{11})$.

In order to obtain nonzero entries in matrix $A_{11}$, we performed a linear regression between $\vec{X}(t+1)$ and data $\underline{X}_{t-l+1}^{t}$ where $l$ is the lag length. Let $\hat{A}_{11}$ be the output of linear regression for time series of length $10000$. According to Wald test [11], for large number of samples, we can obtain $Supp(\hat{A}_{11})$ by setting entry $(i,j)$ to one if $|\hat{A}_{11}(i,j)|>a/2$. In Fig. 2, the error $||Supp(\hat{A}_{11})-Supp(A_{11})||^{2}_{F}$ is averaged over generated random matrices where $||.||_{F}$ is the Frobenius norm of a matrix. As it can be seen, the support of matrix $A_{11}$ can be recovered perfectly as the lag length increases. This trend is expected since the lag length should be at least equal to the order of linear model $l$, in order to have perfect recovery. Moreover, as shown in Fig. 2, for a fixed lag length, the average error is higher for larger $p$. This is because the matrices $\{A_{k}^{*}\}$ become more dense for larger $p$ which leads to higher average error when the right lag length is not selected.

We also examined our proposed criteria in a nonlinear system with three state variables with the following joint dynamics:

[TABLE]

where $\omega_{1}(t)$ and $\omega_{2}(t)$ are i.i.d $\mathcal{N}(0,0.1)$ and $X_{1}(0),X_{2}(0)$ have normal distribution with zero mean and unit variance. The quantity in (5) can be written as a linear combination of some joint entropies. Hence, we can utilize the $K$-nearest neighbor method of [21] to obtain an estimation of the desired quantity. To do so, we generated $1000$ samples of $X_{1},X_{2},Z$ for $t=0,1$. For $K=10$, the numerical results were: $I(X_{1}(1);X_{2}(0)|X_{1}(0))\approx 0.06$, and $I(X_{2}(1);X_{1}(0)|X_{2}(0))\approx 1.68$. From these results, we can infer that $X_{2}(t)$ is 1-step functional dependent on $X_{1}(t-1)$ which is consistent with the system model in (24).

5 Conclusion

We proposed an information-theoretic quantity for identifying causal relations among observed variables for a general $1$-step stochastic dynamical system in the presence of latent variables when there exists no exogenous noise in the latent part. It would be interesting to see if by further imposing some additional constraints on the structure of functional dependencies, it would be possible to recover the inter-connections in the latent sub-graph.