Learning Vector Autoregressive Models with Latent Processes
Saber Salehkaleybar Jalal Etesami Negar Kiyavash Kun Zhang











Learning Vector Autoregressive Models with Latent Processes
    Saber Salehkaleybar^*, Jalal Etesami^*†, Negar Kiyavash^†, Kun Zhang^

^*Coordinated Science Laboratory, University of Illinois at Urbana-Champaign, Urbana, USA.

^†Department of ISE, University of Illinois at Urbana-Champaign
Urbana, USA.

^Department of ECE, University of Illinois at Urbana-Champaign, Urbana, USA.

^Department of Philosophy, Carnegie Mellon University, Pittsburgh, USA.




    December 30, 2023
==============================================================================================================================================================================================================================================================================================================================================================================================================
















We study the problem of learning the support of transition matrix between random processes in a Vector Autoregressive (VAR) model from samples when a subset of the processes are latent. It is well known that ignoring the effect of the latent processes may lead to very different estimates of the influences among observed processes, and we are concerned with identifying the influences among the observed processes, those between the latent ones, and those from the latent to the observed ones. We show that the support of transition matrix among the observed processes and lengths of all latent paths between any two observed processes can be identified successfully under some conditions on the VAR model. From the lengths of latent paths, we reconstruct the latent subgraph (representing the influences among the latent processes) with a minimum number of variables uniquely if its topology is a directed tree. Furthermore, we propose an algorithm that finds all possible minimal latent graphs under some conditions on the lengths of latent paths. Our results apply to both non-Gaussian and Gaussian cases, and experimental results on various synthetic and real-world datasets validate our theoretical results.





§ INTRODUCTION

Identifying causal influences among time series is a problem of interest in many fields. In macroeconomics, for instance, researchers seek to understand what factors contribute to economic fluctuations and how they interact with each other <cit.>. In neuroscience, many researchers focus on learning the interactions between different regions of brain by analyzing neural spike trains <cit.>.


Granger causality <cit.>, transfer entropy <cit.>, and directed information <cit.> are some of the most commonly used measures in the literature to calculate time-delayed dependence structures in time series. 
Measuring the reduction of uncertainty in one variable after observing another variable is the key concept behind such measures. Under certain assumptions, these measures may represent causal relations among the variables <cit.>. In <cit.>, an overview of various definitions of causation is given for time series.

In this work, we study the causal identification problem in VAR models when only a subset of times series is observed.  
More precisely, we assume that the available measurements are a set of random processes X⃗(t)∈ℝ^n which, together with another set of latent random processes Z⃗(t)∈ℝ^m, where m≤ n form a first order VAR model as follows:


    [ X⃗(t+1); Z⃗(t+1) ]
    =[ A_11 A_12; A_21 A_22 ][ X⃗(t); Z⃗(t) ]
    +
    [ ω⃗_X(t+1); ω⃗_Z(t+1) ].


Here we assume that observed data were measured at the right causal frequency of the VAR process; otherwise one may need to consider the effect of the sampling procedure such as subsampling or temporal aggregation <cit.>. Under certain assumptions (e.g., causal sufficiency), the support of the transition matrix corresponds to the causal structure between these processes <cit.>. 
If we ignore the influence of latent processes and just regress X⃗(t+1) on X⃗(t), we may get a wrong estimate of the transition matrix between observed processes (see the example in <cit.>). Hence, it is crucial to consider the presence of latent processes and their influences on the observed processes.






Contributions: The contributions of this paper are as follows: we propose a learning approach that recovers the observed sub-network (support of A_11) by regressing the observed vector X⃗(t+1) on a set of its past observations (not just X⃗(t)) as long as the graph representation of latent sub-network (support of A_22) is a directed acyclic graph (DAG).
We also derive a set of sufficient conditions under which we can uniquely recover the influences from latent to observed processes, (support of A_12) and also the influences among the latent variables, (support of A_22). 
Additionally, we propose a sufficient condition under which the support of the complete transition matrix can be recovered uniquely. 

More specifically, we show that under an assumption on the observed to latent noise power ratio, if neither of the sub-matrices A_12 and A_21 are zero,
it is possible to determine the length of all directed latent paths[A directed path is a latent path if it connects two observed variables and all the intermediate variables on that path are latent.]. 
We refer to this information as linear measurements[This is because it can be inferred from the observational data using linear regression.].
This information reveals important properties of the causal structure among the latent and observed processes, i.e., support of [0, A_12; A_21, A_22].
 We call this sub-network of a VAR model unobserved network. 
 We show that in the case that the unobserved network is a directed tree and each latent variable has at least two parents and two children, a straightforward application of <cit.> can recover the unobserved network uniquely.

Furthermore, we propose Algorithm <ref> that recovers the support of A_22 and A_12 given the linear measurements when only the latent sub-network is a directed tree plus some extra structural assumptions (see Assumption <ref>).
Lastly, we study the causal structures of VAR models in a more general case in which there exists at most one directed latent path of length k≥2 between any two observed processes (see Assumption <ref>).
For such VAR models, we propose Algorithm <ref> that can recover all possible unobserved networks with minimum number of latent processes. Our results apply to both non-Gaussian and Gaussian cases, and experimental results on various synthetic and real-world datasets validate our theoretical results. All proofs can be found in supplemental material.

Related works: 
The problem of recovering latent causal structure for time series has been studied in the literature. Assuming that connections between observed variables are sparse and each latent variable interacts with many observed variables, it has been shown 
 that the transition matrix between observed variables can be identified in a VAR model <cit.>. However, their approach focuses on learning only the observed sub-network.
<cit.> applied a method based on expectation maximization (EM) to infer properties of partially observed Markov processes, without providing theoretical analysis for identifiability. 
<cit.> showed that if the exogenous noises are independent non-Gaussian and additional so-called genericity assumptions hold, then the sub-networks A_11 and a part of A_12 are uniquely identifiable. However, these assumptions may not hold true in a real-world dataset even with three variables <cit.>.
They also presented a result in which they allowed Gaussian noises in their VAR model and obtained a set of conditions under which they can recover  up to 2n  n candidate matrices for A_11.
Their learning approach is also based on EM and approximately maximizes the likelihood of a parametric VAR model with a mixture of Gaussians as noise distribution.
Recently, <cit.> studied a network of processes (not necessary a VAR model) whose underlying structure is a polytree and introduced an algorithm that can learn the entire casual structure (observed and unobserved networks) using a particular discrepancy measure.


Compared to related works, we improve the state of the art for latent recovery by showing the identifiability of a much larger class of structures. Unlike <cit.>, we do not assume the non-Gaussian distribution of the exogenous noises or those genericity assumptions. Moreover, our results do not rely on the assumption that connections between observed variables are sparse or each latent variables interacts with many observed variables as in <cit.>. Furthermore, these works <cit.> can uniquely identify at most a part of transition matrix (A_11 or a part of A_12). 




§ PROBLEM DEFINITION




In this part, we review some basic definitions and our notation. 
Throughout this paper, we use an arrow over the letters to denote vectors.
We assume that the time series are stationary and denote the autocorrelation of X⃗ by γ_X(k):=𝔼[X⃗(t)X⃗(t-k)^T]. 
We denote the support of a matrix A by Supp(A) and use Supp(A)⊆ Supp(B) to indicate [A]_ij=0 whenever [B]_ij=0.
We also denote the Fourier transform of g by ℱ(g) and it is given by ∑_h=-∞^∞ g(h) e^-hΩ j.

In a directed graph G=(V,E) with the node set V and the edge set E, we denote the set of parents of a node v by 𝒫_v:={u: (u,v)∈E} and the set of its children by 𝒞_v:={u: (v,u)∈E}. 
The skeleton of a directed graph G is the undirected graph obtained by removing all the directions in G.






 §.§ System Model


Consider the VAR model in (<ref>). Let ω⃗_X(t)∈ℝ^n and ω⃗_Z(t)∈ℝ^m be i.i.d random vectors with mean zero. For simplicity, we denote the matrix [A_11,A_12;A_21,A_22] by A. 
Our goal is to recover Supp(A) from observational data, i.e., {X⃗(t)}. Rewrite <ref> as follows
 



    X⃗(t+1)=   ∑_k=0^t A^*_kX⃗(t-k) + A_12A_22^t Z⃗(0)+
       ∑_k=0^t-1Ã_kω⃗_Z(t-k) +ω⃗_X(t+1),

where A_0^*:=A_11, A_k^*:=A_12A_22^k-1A_21 for k≥ 1, and Ã_k:= A_12A_22^k. 


We assume that the A_22 is acyclic, i.e., ∃ 0<l≤ m, such that A_22^l=0.



Based on the above assumption, for t≥ l, Equation (<ref>) becomes[Note that the limits of summations in (<ref>) are changed.]
 
 

    X⃗(t+1)=∑_k=0^l A^*_kX⃗(t-k) +∑_k=0^l-1Ã_kω⃗_Z(t-k)+ω⃗_X(t+1).

We are interested in recovering the set {Supp(A_k^*)}_k=0^l because it captures important information about the structure of the VAR model. 
Specifically, Supp(A^*_0)=Supp(A_11); so it represents the direct causal influences between the observed variables and Supp(A^*_k) for k≥1 determines whether at least one directed path of length k+1 exists between any two observed nodes which goes through the latent sub-network.[Herein, we exclude degenerate cases where there is a direct path from an observed node to another one with length k but the corresponding entry in matrix Supp(A^*_k) is zero. In fact, such special cases can be resolved by small perturbation of nonzero entries in matrix A. In the causal discovery literature, this assumption is known as faithfulness <cit.>.] 
We will make use of this information in our recovery algorithm.
We call the set of matrices {Supp(A_k^*)}_k≥0, linear measurements.
In Section 4, we present a set of sufficient conditions under which given the linear measurements, we can recover the entire or most parts of the unobserved network uniquely. 
<ref>


Note that in general, the linear measurements cannot uniquely specify the unobserved network. 
For example, Figure <ref> illustrates two different unobserved networks that both share the same set of linear measurements, A^*_k=0 for k>2 and the only nonzero entries of A^*_1 and A^*_2 are {(3,2)} and {(4,1),(4,2)}, respectively.
































§ IDENTIFIABILITY OF THE LINEAR MEASUREMENTS

As we need the linear measurements for our structure learning, in this section, we study a sufficient condition under which we can recover the linear measurements from the observed processes {X⃗(t)}.

To do so, we start off by rewriting Equation (<ref>) as follows,
 

    X⃗(t+1)=𝒜𝒳⃗_t-l:t+∑_k=0^l-1Ã_kω⃗_Z(t-k)+ω⃗_X(t+1),

 where 𝒜:= [A_0^*,...,A_l^*], and 𝒳⃗_t-l:t:=[X⃗(t);⋯;X⃗(t-l)].
By projecting Ã_kω⃗_Z(t-k) onto the vector space spanned by the observed processes, i.e., {X⃗(t),...,X⃗(t-l)}, we obtain



    Ã_kω⃗_Z(t-k)=∑_r=0^l C_r^s X⃗(t-r) + N⃗_Z(t-k),   0≤ k≤ l-1,

where {N⃗_Z(t-k)} denote the residual terms and {C_r^s} are the corresponding coefficient matrices.
Substituting (<ref>) into (<ref>) implies


    X⃗(t+1)=ℬ𝒳⃗_t-l:t+θ⃗(t+1),

where ℬ:=[B^*_0,...,B^*_l], 

B^*_k:=A_k^*+∑_s=0^l-1 C_k^s, and  θ⃗(t+1):=ω⃗_X(t+1)+∑_k=0^l-1N⃗_Z(t- k).

Note that by this representation, θ⃗(t+1) is orthogonal to 𝒳⃗_t-l:t.

Hence, Equation (<ref>) shows that the minimum mean square error (MMSE) estimator 
can learn the coeffiecient matrix ℬ given the observed processes. More precisely, let Γ_X(l):= 𝔼{𝒳⃗_t-l:t𝒳⃗_t-l:t^T}, then we have 

    ℬ=[γ_X(1),..,γ_X(l+1)]×Γ_X(l)^-1.


Under Assumption <ref>, for the stationary VAR model in (<ref>), we have

    ||B^*_k-A_k^*||_1≤√(n(l-k-1) M/L)||A_12||_2||A_22||_2^k+1,

where 
M:=λ_max(Γ_ω_Z(0)) and L:=λ_min(Γ_X(0)). 
 




This result implies that we can asymptotically recover the support of  {A_k^*}_k=0^l as long as the absolute values of non-zero entries of A_k^* are bounded away from zero by 2√(n(l-k-1)M/L)||A_12||_2||A_22||_2^k+1. Please note that A_11=A_0^*=B_0^* if ||A_12||_2=0.
In Appendix (the second section), we explained how these bounds can be estimated from observational data. 



Let Σ_X=σ^2_X I_n× n and Σ_Z=σ_Z^2 I_m× m be the autocovariance matrices of ω⃗_X(t) and ω⃗_Z(t), respectively. Then, the ratio M/L strictly increases by decreasing σ_X^2/σ_Z^2.




Proposition <ref> implies that when the σ_X^2/σ_Z^2 increases, M/L will decrease, and based on the bound in Proposition <ref>, the estimation error will decrease (it goes to zero as σ_X^2/σ_Z^2 tends to infinity). This shows that recovering the linear measurements is much easier in high σ_X^2/σ_Z^2 regime as illustrated in Figure <ref>. Note that Proposition <ref> stresses a suffiecient condition for recovering the linear measurements. As shown in Figure <ref>, in practice, the actual estimation error is much smaller than the bound in Proposition <ref>. In the next section, we will make use of {Supp(A_k^*)}_k>0 to recover the unobserved network. We assume that the correct linear measurements can be obtained from matrix ℬ.




In order to estimate the support of matrix ℬ from a finite number of samples drawn from the observed processes, say {X⃗(t)}_t=1^T, first we obtain the lag length l in (<ref>) by AIC or FPE criterion (see Chapter 4 in <cit.>). Afterwards, we can estimate the coefficient matrix ℬ, using an empirical estimator for Γ_X(l), {γ_X(h)}_h=1^l+1, and then applying (<ref>). Denote the result of this estimation by ℬ_T. It can be shown that <cit.>, 

√(T)vec(ℬ_T-ℬ)𝒩(0,  Γ^-1_X(l)⊗Σ),

where  denotes convergence in distribution, and Σ is the autocovariance matrix of θ⃗(t). vec(.) transforms a matrix to a vector by stacking its columns and ⊗ is the Kronecker product. 
Having the estimates of Γ_X(l) and Σ, we can test whether the entries of matrix ℬ are greater than the bounds in Proposition <ref>  (see Chapter 3 in <cit.>).














§ LEARNING THE UNOBSERVED NETWORK

Recall that we refer to Supp([0, A_12; A_21, A_22]) as the unobserved network and Supp(A_22) as the latent sub-network. 
We present three algorithms that take the linear measurements {Supp(A^*_k)}_k≥ 0 as their input.
The first algorithm recovers the entire unobserved network uniquely as long as it is a directed tree and  each latent node has at least two parents and two children.
The output of the second algorithm is Supp([0, A_12; A_21, A_22]), where Supp(A_21)⊆ Supp(A_21). 
This is guaranteed whenever the latent sub-network is a directed tree and some extra conditions are satisfied on how the latent and observed nodes are connected.
The third algorithm finds the set of all possible networks with minimum number of latent nodes that are consistent with the measurements.
This algorithm is able to do so when there exists at most one directed latent path of any arbitrarily length between two observed nodes. A directed path is latent if all the intermediate variables on that path are latent.



 §.§ Unobserved Network is a Directed Tree

Authors in <cit.> introduced a necessary and sufficient condition for recovering a weighted directed tree uniquely from a valid distance matrix D defined on the observed nodes,[The skeleton of the recovered tree is the same as the original one but not necessary the weights.]  and also proposed a recovery algorithm. 
The condition is as follows: every latent node must have at least two parents and two children. 
A matrix D, in <cit.>, is a valid distance matrix, when [D]_ij equals the sum of all the weights of those edges that belong to the directed path from i to j, and [D]_ij=0, if there is no directed path.

The algorithm in <cit.> has two phases. In the first phase, it creates a directed graph among the observed nodes with the adjacency matrix Supp(D). 
In the second phase, it recursively finds and removes the circuits by introducing latent nodes for each circuit.[In a directed graph, a circuit is a cycle after removing all the directions.]

For more details, see <cit.>.




































In order to adopt <cit.>'s algorithm for learning the unobserved network, we introduce a valid distance matrix using our linear measurements as follows, D_ij=k+1 if  [Supp(A_k^*)]_ji≠ 0 and 0, otherwise.


   
      
   

Recall that [Supp(A_k^*)]_ji indicates whether there exists a directed latent path from i to j of length k+1 in the unobserved network. From theorem 8 in <cit.>, it is easy to show that the unobserved network can be recovered uniquely from above distance matrix if its topology is a directed tree and every latent node has at least two parents and two children.




 §.§ Latent Sub-network Is a Directed Tree



We denote the subset of observed nodes that are parents of a latent node h by 𝒫^O_h and denote the subset of observed nodes for which h is a parent, by 𝒞^O_h. We further denote the set of all leaves in the latent sub-network by ℒ.


We consider learning an unobserved network G that satisfies the following assumptions.


Assume that the latent sub-network of G is a directed tree. Furthermore, for any latent node h in G, (i) 𝒫^O_h⊈∪_h≠ j𝒫^O_j and, (ii) if h is a leaf of the latent sub-network, then 𝒞^O_h⊈∪_ i∈ℒ,i≠ h𝒞^O_i.



This assumption states that the latent sub-network of G must be a directed tree such that each latent node in G has at least one unique parent in the set of observed nodes. That is, a parent who is not shared with any other latent node. Furthermore, each latent leaf has at least one unique child among the observed nodes.
For instance, when Supp(A_22) represents a directed tree and both Supp(A_12) and Supp(A_21) contain identity matrices, Assumption <ref> holds. As we will see later in Experimental Results (Figure <ref>), a large portion of randomly generated graphs satisfy Assumption <ref>.

Figure <ref> illustrates a simple network that satisfies Assumption <ref> in which the unique parents of latent nodes a, b,c, and d are {1}, {3}, {2}, and {4}, respectively. The unique children of latent leaves c and d are {5} and {2,4}, respectively.


Among all unobserved networks that are consistent with the linear measurements induced from (<ref>), any graph G that satisfies Assumption <ref> has the minimum number of latent nodes.











Note that if Assumption <ref> is violated, one can find many unobserved networks that are consistent with the linear measurements but are not minimum (in terms of the number of latent nodes).
For example, the network in Figure <ref> satisfies Assumption <ref> (ii) but not (i). 
Figure <ref> depicts an alternative network with the same linear measurements as the network in Figure <ref> but it has fewer number of latent nodes.
Similarly, the graph in Figure <ref> satisfies Assumption <ref> (i) but not (ii). Figure <ref> shows an alternative graph with one less latent node.





































Consider an unobserved network G with adjacency matrix Supp([0, A_12; A_21, A_22]). 
If G satisfies Assumption <ref>, then its corresponding linear measurements uniquely identify G upto Supp([0, A_12; A_21, A_22]), where Supp(A_21)⊆ Supp(A_21).


Figure <ref> gives an example of a network satisfying Assumption <ref> and an alternative network, Figure <ref>, with the same linear measurements which departs from the Figure <ref> in the A_21 component.




Next, we propose the directed tree recovery (DTR) algorithm that takes the linear measurements of an unobserved network G satisfying Assumption <ref> and recovers G upto the limitation in Theorem <ref>. This algorithm consists of three main loops. 
Recall that Assumption <ref> implies that each latent node has at least one unique observed parent.
The first loop finds all the unique observed parents for each latent node (lines: 4-11). 
The second loop reconstructs Supp(A_22) and Supp(A_12) (lines: 12-17). And finally, the third loop constructs Supp(A_21) such that Supp(A_21)⊆ Supp(A_21) (lines: 18-22).

The following lemma shows that the first loop of Algorithm <ref> can find all the unique observed parents from each latent node. To present the lemma, we need the following definitions. 



For an observed node i, we define


    l_i:=max{k:[A^*_k-1]_si≠0,  for some s},
        R_i:={j:[A^*_l_i-1]_ji≠0},  M_i:={(j,r) : [A^*_r-1]_ji≠0}.



In the above equations, l_i denotes the length of longest directed latent path that connects node i to any observed node. R_i is the set of all observed nodes that can be reached by i with a directed latent path of length l_i and set M_i consists of all pairs (j,r) such that there exists a directed latent path from i to j with length r.


Under Assumption <ref>, an observed node i is the unique parent of a latent node if and only if for any other observed node j s.t. l_i=l_j, we have 
(R_j⊈R_i) ∨ (R_j=R_i ∧ M_i⊆ M_j).





In the first loop, if there exist multiple unique parents of a latent node (for instance, node 2 and node 3 in Figure <ref>), we pick the one with a minimum index (lines: 7-9). 

 The second loop recovers Supp(A_22) based on the following observation. If a latent node h_k is the parent of latent node h_s, then h_k can reach all the observed nodes in R_s, i.e., R_s⊆ R_k  and l_k=l_s+1 (line: 13).
 Furthermore, Supp(A_12) can be recovered using the fact that an observed node j is a children of a latent node h_s, if a unique parent of h_s, e.g., s, can reach j by a directed latent path of length 2 (line: 16).
 Finally, the third loop reconstructs Supp(A_21) by adding an observed node i to the parent set of latent node h_j, if i can reach all the observed nodes that a unique parent of h_j, e.g., j, reaches (lines: 18-22).
 



Suppose network G satisfies Assumption <ref>. Then given its corresponding linear measurements, Algorithm <ref> recovers G upto the limitation in Theorem <ref>. 








 §.§  Learning More General Unobserved Networks with Minimum Number of Latent Nodes

In general, the latent sub-network may not be a tree or there may not be a unique minimal unobserved network consistent with the linear measurements (see Figure <ref>). 
Hence, we try to find an efficient approach to recovering all possible minimal unobserved networks under some conditions.
In fact, without any extra conditions, finding a minimal unobserved network is NP-hard.


Finding an unobserved network that is both consistent with a given linear measurements and has a minimum number of latent nodes is NP-hard.

Below, after some definitions, we propose the Node-Merging (NM) algorithm that returns all possible unobserved networks with minimum number of latent nodes under the following assumption. 


Assume that there exists at most one directed latent path of each length between any two observed nodes.







For example, the graph in Figure <ref> satisfies this assumption but not the one in Figure <ref>. This is because there are two directed latent paths of length 2 from node 5 to node 4.


(Merging)
We define merging two nodes i' and j' in graph G as follows: remove node j^' and the edges between i' and j', and then give all the parents and children of j'to i^'. 
We denote the resulting graph after merging i' and j' by Merge(G,i',j').
We say that two nodes i^' and j^' are mergeable if Merge(G,i',j') is consistent with the linear measurements of G.



 (Connectedness) Consider an undirected graph G̅ over the observed nodes which is constructed as follows: there is an edge between two nodes i and j in G̅, if there exists k≥ 1 s.t. Supp([A_k^*]_ij)=1 or Supp([A_k^*]_ji)=1;
We say that two observed nodes i and j are “connected" if there exist a path between them in G̅.




It can be seen that if pairs i,j and j,k are connected then node i,k are also connected. 
We then define a connected class as a subset of observed nodes in which any two nodes  are connected.


Initialization:
We first find the set of all connected classes, say S_1,S_2,...,S_C. 
For each class S_c, we create a directed graph G_0,c that is consistent with the linear measurements. 
To do so, for any two observed nodes i,j∈ S_c, if [A_r^*]_ji≠ 0, we construct a directed path with length r+1 from node i to node j by adding r new latent nodes to G_0,c. 


Merger:
In this phase, for any G_0,c from the initialization phase, we merge its latent nodes iteratively until no further latent pairs can be merged. 
Since the order of mergers leads to different networks with minimum number of latent nodes, the output of this phase will be the set of all such networks.
Algorithm <ref> summarizes the steps of NM algorithm.
In this algorithm, subroutine Check(G,i',j') checks whether two nodes i' and j' are mergeable.


Under Assumptions <ref> and <ref>, the NM algorithm returns the set of all networks that are consistent with the linear measurements and have minimum number of latent nodes.








§ EXPERIMENTAL RESULTS






  
Synthetic Data: We considered a directed random graph, denoted by DRG(p,q), such that there exists a directed link between an observed and latent node with probability p,  independently across all pairs, and there is a directed link between two latent nodes with probability q. If there is a link between two nodes, we set the weight of that link uniformly from [-a,a]. 


 
We utilize the method described in Section 3 to estimate linear measurements with a significance level of 0.05.
In order to evaluate how well we can estimate the linear measurements, we generated 1000 instances of DRG(0.4,0.4) with n+m=100, Σ_X=0.1I_n× n,Σ_Z=0.1I_m × m, and a=0.1. 
The length of the time series was set to T=1000. Let Supp(Â_11) be the estimate of support of A_11.
In Figure <ref>, the expected estimation error, i.e. || Supp(Â_11)-Supp(A_11)||^2_F/n^2, is computed, where ||.||_F is the Frobenius norm. One can see that the
estimation error decreases as the number of observed variables increases.





We also studied the effect of the observed to latent noise power ratio (OLNR), σ_X^2/σ_Z^2, on ||B_0^*-A_0^*||_1, and compared it with the bound given in Proposition <ref>. We generated 1000 instances of DRG(0.05,0.05) with n=5, m=5, and a=0.1. As it can be seen in Figure <ref>, the average estimation error decreases as OLNR increases, as expected from Proposition <ref>.

We investigated what percentage of instances of the random graphs satisfy Assumption <ref>. 
We generated 1000 instances of DRG(p,1/n) with n=100, and p∈[0.04,0.2].
 In Figure <ref>, the probability of satisfying Assumption <ref>, P_sat., is depicted versus  p for different numbers of latent variables in the VAR model. 
 For larger m, it is less likely to see a unique observed parent for each latent node and thus P_sat. decreases.
 For a fixed m, the same phenomenon will occur if we increase p when p is relatively large. 
 Furthermore, for small p, there might exist some latent nodes that have no observed parent or no observed children. 


We also evaluated the performance of the NM algorithm in random graphs. 
We generated 1000 instances of DRG(1/2n,1/2n) with n=10,...,100 and m=n/2, and computed the linear measurements. 
To save time, if for a class of connected nodes the number of latent nodes generated in the initial phase exceeds 40, 
we supposed that the corresponding instance cannot be recovered efficiently in time and did not proceed to the merging phase. 

Figures <ref> and <ref> depict the percentage of instances in which the algorithm can recover all possible minimal unobserved networks and the average run time (in seconds) of the algorithm, respectively.[We performed the experiment on a Mac with 2× 2.4 GHz 6-Core Intel Xeon processor and 32 GB of RAM.] 
This plot shows that we can recover all possible minimal unobserved networks for a large portion of instances efficiently even in relatively large networks.
 
 











  
US Macroeconomic Data:
We considered the following set of time series from the quarterly US macroeconomic data for the period from 31-Mar-1947 to 31-Mar-2009 collected from the St. Louis Federal Reserve Economic Database (FRED) <cit.>:
GDP, GDPDEF, COE, HOANBS, TB3MS, PCEC, GPDI.

Assuming that the underlying dynamics is linear (Eq. (<ref>)), we considered the estimated VAR model over all variables as the ground truth.
Then, we selected four arbitrary times series as observed processes and computed Supp(Â_11). 
We divided the 7  4=35 possible selections into two classes: 1) high power, where tr(𝔼{ω_X(t)ω_X(t)^T})>τ for a fixed threshold τ; 
2) low power:  where tr(𝔼{ω_X(t)ω_X(t)^T})<τ.
In this experiment, we set τ=0.02.
 In Figure <ref>, we plotted the histograms of || Supp(Â_11)-Supp(A_11)||^2_F for these two classes. 
As it can be seen, in the high power regime, most of the possible selections have small estimation errors.

We also considered the following six time series of US macroeconomic data during 1-Jun-2009 to 31-Dec-2016 from the same database: GDP, GPDI, PCEC, TBSMS, FEDFUND, and GS10.  We obtained the causal structure among these six time series by fitting a VAR model on all of them and considered the result as our ground truth (see Figure <ref>). 
Then, we removed GPDI from the dataset and considered the remaining five time series as observe processes and checked whether the influences from the “latent" process (GPDI) can be corrected estimated. 
r5cm

 !<0cm,0cm>;<1.3cm,0cm>:<0cm,1.2cm>:: 
!(-.9,.4) *+GDP="gdp"   
  !(0,-.2) *+GPDI∘="gpdi"	
	!(2,.4) *+PCEC="pcec"   	 		
!(1,-.2) *+GS10="gsio" 
!(2,-.7) *+TBSMS="tbsms"  
!(-.9,-.7) *+FED="fed" 
"gdp":@/^/"gsio" "gdp":"fed" "gdp":@/_.9cm/"tbsms"
"gdp":"gpdi" "pcec":@/_/"gpdi" "pcec":"gsio" "pcec":"tbsms" "gpdi":"fed" "gpdi":"tbsms" "pcec":@/_-.7cm/"fed"
 
US macroeconomic data.


 
We estimated the linear measurements and gave them as an input  to Algorithm <ref>, which successfully recovered the ground truth (the estimated structure, in which the latent process is denoted by a circle, is identical to that in Figure <ref>).






  
Dairy Prices:
A collection of three US dairy prices has been observed monthly from January 1986 to December 2016 <cit.>: milk, butter, and cheese prices. 
r4.7cm

 !<0cm,0cm>;<1.2cm,0cm>:<0cm,.8cm>:: 
!(-.5,-.6) *+Milk ="x1"   
  !(1,.3) *+Butter ="z"	
   !(2.3,-.6) *+Cheese="x2"  
"x1":"x2" "x1":"z"  "z":"x2"  "x2":@/_-.4cm/"x1"
 
Dairy prices
	

 
We estimated the VAR model on all the time series with lag length l=1  and considered the resulting graph as our ground truth (see Figure <ref>).
 Next, we omitted  the butter prices from the dataset and considered  
 the milk and cheese prices as observed processes.
 The estimated linear measurements were: Supp(A_0^*)=Supp(A_11)=[1,1;1,0] and Supp(A_1^*)=[0,0;1,0]. 
Algorithm <ref> correctly recovered  the true causal graph using this linear measurements. Note that the genericity assumptions in <cit.> do not hold true for this data set (see Experiments section).






  
West German Macroeconomic Data:We considered the quarterly West German consumption expenditures X_1, fixed investment X_2, and disposable income X_3, during 1960-1982 <cit.>. 
r4.7cm



 !<0cm,0cm>;<1.2cm,0cm>:<0cm,.8cm>:: 
!(-.5,-.6) *+Expend="x1"   
  !(1,.5) *+Income ="z"	
   !(2.3,-.6) *+Invest="x2"  
"x1":"x2" "x1":"z"  "z":"x2"  "z":@/_.4cm/"x1"
 
West German macroeconomic data.
	

 
Similar to the previous experiment with dairy prices, we first obtained the entire transition matrix among all the process. Figure <ref> depicts the resulting graph.
Next, we considered X_3 to be latent and used {X_1,X_2} to estimate the linear measurements Supp(A_0^*)=Supp(A_11)=[0,0;1,1] and Supp(A_1^*)=[1,0;1,0]. Using this linear measurements, Algorithm <ref> recovered the true network in Figure <ref> correctly.





§ CONCLUSION AND FUTURE WORK


We considered the problem of estimating time-delayed influence structure from partially observed time series data. 
Our approach consisted of two parts: First, we studied sufficient conditions under which certain aspects of the influence structure of the underlying system are identifiable. Second, we proposed two algorithms that recover the influence structures satisfying the sufficient conditions given in the first part. The proposed algorithms can construct the observed sub-network (support of A_11), the causal influences from latent to observed processes (support of A_12), and also the causal influences among the latent variables (support of A_22), uniquely under a set of sufficient conditions. As a future direction, we plan to extend our results to the case that A_22 might have cycles. In the paper, we have seen examples showing that unique recovery is not possible if any conditions of Assumption <ref> are violated. These conditions are a good starting point for the case that we have cycles in A_22.






aaai












§ PROOF OF PROPOSITION <REF>

We project the vector A_r+1:l-1[ω⃗_Z(t-r-1);⋯;ω⃗_Z(t-l+1)] onto X⃗(t-r) as follows: 


    A_r+1:l-1[ ω⃗_Z(t-r-1);           ⋮; ω⃗_Z(t-l+1) ]=C_r
    X⃗(t-r)
    +[ N⃗_Z(t-r-1);           ⋮; N⃗_Z(t-l+1) ],


where A_r+1:l-1=diag(Ã_r+1,...,Ã_l-1), and C_r is a block matrix with C_r^s as its sth block for s=0,...,l-r-2. Please note that ω⃗_Z(t-r) is orthogonal to X⃗(t-k) for k≥ r. 
Since N⃗_Z and X⃗(t-r) are orthogonal, we can see 


    ||A_r+1:l-1Γ_ω_Z(l-r-2)A_r+1:l-1^T||_2≥ ||C_rΓ_X(0)C_r^T||_2.


Using (<ref>) and the relationship between ℓ_2 and ℓ_1  norms of a matrix, we obtain

    λ_max(Γ_ω_Z(0))||A   _r+1:l-1||_2^2≥
       λ_min(Γ_X(0))||C_r||_1^2/(n(l-r-1)),

 where λ_min(·) and λ_max(·) denote the minimum and maximum eigenvalues of a given matrix, respectively. Please note that ω⃗_Z(t) is white noise and thus we have: λ_max(Γ_ω_Z(l-r-2))=λ_max(Γ_ω_Z(0)). Using the fact that A_r+1:l-1 is diagonal and ||A_22||_2<1, we obtain


    ||C_r||_1    ≤√(n(l-r-1)M/L)||A_12||_2max_r+1≤ k≤ l-1||A_22||^k_2
       ≤√(n(l-r-1)M/L)||A_12||_2||A_22||_2^r+1.


where M:=λ_max(Γ_ω_Z(0)) and L:=λ_min(Γ_X(0)).


From (<ref>), we have B^*_r-A_r^*= ∑_s=0^l-r-2 C_r^s. This implies that 
||B^*_r-A_r^*||_1≤ ||C_r||_1. Combining this inequality and the bound in (<ref>) concludes the result.



§ ESTIMATING THE BOUNDS IN PROPOSITION <REF>

The bound √(n(l-r-1)M/L)||A_12||_2||A_22||_2^k+1 can be estimated as follows:


  * The lag length l in (<ref>) can be obtained from AIC or FPE criterion (see chapter 4 in <cit.>).

  * We can estimate L by observation vector X⃗(t). We also consider a bound σ^2_max,Z on the maximum variance of exogenous noises in latent part.

  * We assume a bound on ||A_12||_2≤ρ_12 and ||A_22||_2≤ρ_22<1. 


In summary, an upper bound would be: √(n(l-r-1)σ^2_max,Z/L)ρ_12ρ_22^k+1. Suppose that absolute values of nonzero entries of A_k^* are greater than a_min,k. We can recover the support of matrix A_k^* successfully if

    4n(l-r-1)ρ_12^2/a_min,k^2(ρ_22)^2(k+1)≤L/σ_max,Z^2.



 
 

§ PROOF OF PROPOSITION <REF>

The spectral density of matrix γ_X(h) can be computed as follows:

    ℱ(γ_X)= σ_X^2 F_X(Ω) F_X(Ω)^H+ σ_Z^2 F_Z(Ω) F_Z(Ω)^H

where F_X(Ω)= [e^jΩI_n× n -A_11 -∑_k=0^l-1 A_k^*e^-kjΩ]^-1, F_Z(Ω)=F_X(Ω)(A_12∑_k=0^l-1 A_22^k × e^-kjΩ), and H denotes Hermitian of a matrix. Thus, we have:

    Γ_X(0)=1/2π∫_0^2 πℱ(γ_X) dΩ= σ_X^2 F_X^0+σ_Z^2 F_Z^0,

where F_X^0=1/(2π)∫_0^2π F_X(Ω) F_X(Ω)^H dΩ and F_Z^0=1/(2π)∫_0^2π F_Z(Ω) F_Z(Ω)^HdΩ.

We define the function ψ_σ_X/σ_Z(v):=v⃗^TΓ_X(0)v⃗/σ_Z^2= (σ_X^2/σ_Z^2)F_X^0+ F_Z^0 where v⃗ is a unit vector. Suppose that v⃗^* minimizes the function ψ_σ_X/σ_Z(.). By the definition of L and M, the ratio M/L is equal to 1/ψ_σ_X/σ_Z(v⃗^*). Now if we decrease σ_X/σ_Z to σ'_X/σ'_Z, then we have: ψ_σ'_X/σ'_Z(v⃗^*)<ψ_σ_X/σ_Z(v⃗^*). Moreover, for the optimal solution v⃗'^* of ψ_σ'_X/σ'_Z(.), we know that: ψ_σ'_X/σ'_Z(v⃗'^*)≤ψ_σ'_X/σ'_Z(v⃗^*). Thus, we can conclude that: 1/ψ_σ'_X/σ'_Z(v⃗'^*)>1/ψ_σ_X/σ_Z(v⃗^*). 



§ PROOF OF THEOREM <REF>

First, we show such G has a minimum number of latent nodes. We do this by means of contradiction. 
But first observe that since the latent subnetwork of G is a directed tree, we can assign a non-negative number l_h to latent node h that represents the length of longest directed path from h to its latent descendants. Clearly, all such descendants are leaves which we denote them by L̃_h.
For instance, if the latent subnetwork of G is a→ b→ c, then l_a=2 and L̃_a={c}. 

Suppose that G contains m latent nodes {h_1,...,h_m} and there exists another network G_1 (not necessary with tree-structure induced latent subgraph), with m_1<m number of latent nodes that it is also consistent with the same linear measurements as G.
Due to assumption (i), there is at least m distinct observed nodes that have out-going edges to the latent subnetwork. More precisely, each h_i has at least a unique observed node as its parent. We denote a unique observed parent of node h_i by o_i.

Because m_1<m, there exists at least one observed node in O̅:={o_1,...,o_m} that has shared its latent children with some other latent nodes in G_1. Among all such observed nodes, let o_i^* to be the one whose corresponding latent node in G, (h_i^*), has maximum l_h_i^*.[If there are several such observed node, let o_i^* to be one of them.]
Furthermore, let Ĩ_i^*⊂{1,...,m}∖{i^*} to be the index-set of those observed nodes that o_i^* has shared a latent child with them in G_1.

By the choice of o_i^*, we know that l_h_j≤ l_h_i^* for all j∈Ĩ_i^* and if for some 1≤ k≤ m, l_h_k>l_h_i^*, then o_k has not shared its latent child in G_1 with any other observed nodes in O̅. Moreover, there should be at least a latent node h_j^* where j^*∈Ĩ_i^* such that l_h_j^*= l_h_i^*. Otherwise, G_1 will not be consistent with the linear measurements of G.
Let Ĩ_**:={j: l_h_j= l_h_i^*}∩Ĩ_i^*. Because o_i^* shares its latent children with ∪_j∈Ĩ_**o_j in G_1 and both G and G_1 consistent with the same linear measurements,  the following holds in graph G,

    𝒞^O_L̃_h_i^*(G)⊆∪_j∈Ĩ_**𝒞^O_L̃_h_j(G),

where 𝒞^O_L̃_h_j(G) indicates the set of observed children of the set L̃_h_j. 
This indeed contradicts assumption (ii).



§ PROOF OF THEOREM <REF>

First, we require the following definition. 
For a network G with corresponding latent sub-network that is a tree, we define U_k(G):={h∈ G:  l_h=k}.
To prove the equivalency, suppose there exists another network G_2 such that its latent sub-network is a tree and has a minimum number of latent nodes.
Let {h_1,...,h_m} to denote the latent nodes in G. 
Since G satisfies Assumption (i), for every latent node h_i there exists a unique observed node o_i such that o_i∈𝒫^O_h_i(G) and o_j∉𝒫^O_h_i(G) for all j≠ i.

Since both G and G_2 are consistent with the same linear measurement, it is easy to observe that if h_i∈ U_k(G), then o_i must have at least a latent child in G_2, say h'_i, such that l_h_i=l_h'_i. Note that l_h_i is computed in G and  l_h'_i in G_2.
Moreover, we must have: 

    𝒞_L̃_h_i^O(G)=⋃_h'∈ H'(o_i)∩ U_l_h_i(G_2)𝒞^O_L̃_h'(G_2),

where H'(o_i) denotes the set of latent nodes in G_2 that have o_i as their observed parent.
 In other words, observed nodes that can be reached by a directed path of length l_h_i+2 from o_i should be the same in both graph G and G_2.
This results plus the fact that G satisfies Assumption (ii), imply: 


I) For every h_i∈ U_k(G), there exists a unique latent node h'_i∈ U_k(G_2), such that o_i∈𝒫^O_h'_i(G_2) and o_j∉𝒫^O_h'_i(G_2) for all j≠ i, and

    𝒞^O_L̃_h(G)=𝒞_L̃_h'_i^O(G_2).

Using I) and knowing that both G and G_2 have the same number of latent nodes, we obtain:

II) |U_k(G)|=|U_k(G_2)|, for all k.


Using I) and II), we can define a bijection ϕ between the latent subnetworks of G and G_2 as follows ϕ(h_i)=h'_i. 
Using this bijection and Assumption (ii) of G conclude that if h∈ U_k(G) is the common parent of {h_j_1,...,h_j_s}⊆ U_k-1(G), then  ϕ(h)∈ U_k(G_2) should be the common parent of {ϕ(h_j_1),...,ϕ(h_j_s)}⊆ U_k-1(G_2) and the proof is complete.



§ PROOF OF LEMMA <REF>

Suppose that o_i is the unique observed node of a latent node h_i. Then, for any o_j such that l_i=l_j, if h_i is not a child of o_j, then from assumption ii we have R_j⊈R_i. If h_i is a child of o_j, then since we know that l_i=l_j, M_i⊆ M_j and R_i=R_j. 

Now, suppose that the observed node o_i satisfies conditions but it is not unique parent of any latent node. Let h_i and h_i^' be children of o_i. At least one of them, say node h_i, can reach an observed node by a path of length l_i-1. If h_i^' has the same property, then consider the unique observed parent of h_i^', say node o_j. Based on Assumption (ii), we have R_j⊆ R_i, which is in contradiction with the assumption that node o_i satisfies conditions of Lemma. Moreover, if h_i^' does not have a path to observed node with a length of l_i-1, then for any observed parent of h_i, one of the conditions in the Lemma is not satisfied. Thus, the proof is complete.



§ PROOF OF PROPOSITION <REF>

Notice that the first loop in Algorithm <ref> uses the result of Lemma <ref> and finds all the latent nodes and their corresponding unique observed parents.
The next loop uses the fact that the latent sub-network is a tree and also it satisfies Assumption <ref>. Hence, if there exist two latent nodes h and h', one with depth l and the other one with depth l+1, such that R_h⊆ R_h', then h' must be the parent of h in the latent sub-network. 

Moreover, since each latent node has a unique observed parent, using A^*_1, Algorithm <ref> can identify all the observed children of a latent node.
Finally, the last loop in this algorithm locates the rest of observed nodes as the input of the right latent nodes. The algorithm does it by using the fact that if an observed node i shares a latent child with another observed node j∈ U, then M_j⊆ M_i.
Clearly, if the true unobserved network satisfies Assumption <ref>, the output of this algorithm will have a latent sub-network that is a tree and consistent with the linear measurement. Thus, by the result of Theorem <ref>, it will be the same as the true unobserved network up to some permutations in Supp(A_21).




§ PROOF OF THEOREM <REF>

Consider the instance of the problem where A_22=0_m× m. Without loss of generality, we can assume that entries of A_12 and A_21 are just zero or one. Thus, we need to find [A_12]_n× k and [A_21]_k× n such that Supp(A_12A_21)=Supp(A_1^*) and k is minimum. We will show that the set basis problem <cit.> can be reduced to the decision version of finding the minimal unobserved network which we call it the latent recovery problem. But before that, we define the set basis problem:

The Set Basis Problem <cit.>: given a collection 𝒞 of subsets of a finite set U={1,⋯,n} and an integer k, decide whether or not there is a collection ℬ⊆ 2^U of at most k sets such that for every set C∈𝒞, there exists a collection ℬ_𝒞⊆ℬ where ⋃_B∈ℬ_C B=C.

Any instance of the basis problem can be reduced to an instance of latent recovery problem. To do so, we encode any set C in collection 𝒞 to a row of A_1^*=A_12A_21 where i-th entry is equal to one if i∈ C, and otherwise zero. It is easy to verify that the rows of matrix A_21 correspond to sets in collection ℬ if there exist a solution for the basis problem. Since the basis problem is NP-complete, we can conclude that finding the minimal unobserved network is NP-hard.





§ PROOF OF THEOREM <REF>

Consider a minimal unobserved network G_min. Pick any latent node i^' which its in-degree or out-degree is greater than one. Let V^-_i^' and V^+_i^' be the sets of nodes that are going to and incoming from node i^', respectively. We omit the node i^' and create |V_i^'^-|×|V_i^'^+| latent nodes {i^'_j^'k^'| j^'∈ V_i^'^-, k^'∈ V_i^'^+}. We also add a direct link from node j^'∈ V_i^'^- to i^'_j^'k^' and from i^'_j^'k^' to k^'∈ V_i^'^+ in order to be consistent with measurements. We continue this process until there is no latent node with in-degree or out-degree greater than one. Since there exists at most one path with length k from any observed node to another observed node, the resulted graph is exactly equal to graph G_0. Hence we can construct the minimal graph G_min just by reversing the process of generating latent nodes from G_min to merging latent nodes from G_0. But the NM algorithm consider all the sequence of merging operations. Thus, G_min would be in the set 𝒢_out and the proof is complete.


