Causal Discovery from Temporally Aggregated Time Series

Mingming Gong∗† , Kun Zhang† , Bernhard Schölkopf‡ , Clark Glymour† , Dacheng Tao]∗
Centre for Artificial Intelligence, FEIT, University of Technology Sydney, NSW, Australia†
Department of Philosophy, Carnegie Mellon University, Pittsburgh, USA‡
Max Plank Institute for Intelligent Systems, Tübingen, Germany]
School of Information Technologies, FEIT, University of Sydney, NSW, Australia

Abstract
Discovering causal structure of a dynamical system from observed time series is a traditional
and important problem. In many practical applications, observed data are obtained by applying subsampling or temporally aggregation to theoriginal causal processes, making it difficult to
discover the underlying causal relations. Subsampling refers to the procedure that for every k consecutive observations, one is kept, the rest being
skipped, and recently some advances have been
made in causal discovery from such data. With
temporal aggregation, the local averages or sumsof k consecutive, non-overlapping observationsin the causal process are computed as new observations, and causal discovery from such data iseven harder. In this paper, we investigate how torecover causal relations at the original causal frequency from temporally aggregated data when k
is known. Assuming the time series at the causal
frequency follows a vector autoregressive (VAR)
model, we show that the causal structure at thecausal frequency is identifiable from aggregatedtime series if the noise terms are independent andnon-Gaussian and some other technical conditions
hold. We then present an estimation method basedon non-Gaussian state-space modeling and evaluate its performance on both synthetic and real
data.

1

INTRODUCTION

Causal modeling (Spirtes et al., 2001; Pearl, 2000) of timeseries data has been widely applied in many fields such as
econometrics (Ghysels et al., 2016), neuroscience (Zhouet al., 2014), and climate science (Van Nes et al., 2015).
Classical causal discovery approaches, e.g., Granger causality test (Granger, 1969), usually assume that the data measurement frequency matches the true causal frequency of theunderlying physical process. However, since the true causal
frequency is usually unknown, the time series data are often
measured at the frequency lower than the causal frequency.
For example, some econometric indicators such as GDP
and non-farm payroll are usually recorded at quarterly andmonthly scales. Causal interactions between the processes,
however, may take place at the weekly or fortnightly scales(Ghysels et al., 2016). In neuroscience, imaging technologies have relatively low temporal resolutions, while manyhigh frequency neuronal interactions are important for understanding neuronal dynamics (Zhou et al., 2014). In thesesituations, the available observations have a lower resolution
than the underlying causal process.
There are two typical schemes to generate low-resolution or
low-frequency data from high-frequency ones (Silvestrini& Veredas, 2008; Marcellino, 1999). One is by subsampling: for every k consecutive observations, one is kept,
the rest being skipped. The other is temporally aggregation, i.e., taking the local averages or sums of k consecutive,
non-overlapping observations from the underlying causal
process as new observations. For instance, the time series
of interest, money supply, and temperature are usually obtained by subsampling; in contrast, the U.S. nominal GDP
was obtained by aggregation – it refers to a total number ofdollars spent over a time period.
Numerous contributions have been made on analyzing theeffects of the above two schemes to generate low-resolution
data on the properties of the time series such as estimatedcausal relations and exogeneity (Tiao, 1972; Weiss, 1984;
Granger, 1987; Marcellino, 1999; Rajaguru & Abeysinghe,
2008). These studies found that temporal aggregation canlead to errors in the estimated causal relations if not properly addressed. For example, Breitung & Swanson (2002)
examined the impact of temporal aggregation on Granger
causality in vector autoregressive (VAR) models and foundthat the results of Granger causal analysis heavily dependon temporal aggregation.
Recovering the high frequency causal relations from temporally aggregated data is a very hard problem due to information loss in the aggregation process. A classical way todiscover high frequency causal relations from temporally
aggregated data is to first disaggregate the low frequency
time series to high frequency ones and then apply standard
causal discovery methods on the disaggregated data. Temporal disaggregation of low resolution time series has been
extensively studied in the econometric and statistical literature (Boot et al., 1967; Stram & Wei, 1986; Harvey &Chung, 2000; Moauro & Savio, 2005; Proietti, 2006), which
is clearly an even harder problem than discovering causal
relations.
Recently, a set of methods have been proposed to estimate
the causal structure at the causal frequency from subsampled
data without resorting to disaggregation techniques (Hyttinen et al., 2016; Gong et al., 2015; Plis et al., 2015a; Danks& Plis, 2013). Plis et al. (2015a,b) first inferred the causal
structure from the subsampled data, and then searched forthe causal structure at the causal frequency from the causal
structure inferred in the first step. Based on this framework,
Hyttinen et al. (2016) proposed a much faster inference
method using a general purpose Boolean constraint solver.
Gong et al. (2015) proposed a model-based approach andexamined the identifiability of the underlying vector autoregressive model (VAR) at the causal frequency from subsampled time series. They showed that the causal transition
matrix is identifiable if the noise process is non-Gaussian.
This work was recently extended to mixed frequency databy structural VAR modeling (Tank et al., 2017) . However,
how to estimate causal relations from aggregated data still
remains unknown.
Compared to subsampling, temporal aggregation is perhaps
more widely used to produce low-resolution time series,
especially in economics and finance. However, the effectof temporal aggregation is more complex, and accordingly
it is technically more difficult to recover the underlying
causal relations from such data. Specifically, because thenoise terms are generated by a larger number of independent
components and thus the mixing matrix contains a more
complicated structure, the estimation is both statisticallyand computationally harder.
The objective of this paper is to seek a possible solution
to this problem, by studying the theoretical identifiabilityof the underlying causal relations and developing a practical causal discovery algorithm. Following (Gong et al.,2015), we assume that the high-frequency data follow a
VAR model, the error terms are non-Gaussian, and there areno confounders (Geiger et al., 2015). We show that the original causal relation can be estimated from the aggregateddata with known k, under a set of technical conditions.
Moreover, we propose an estimation method based on nonGaussian state-space modeling of the aggregated data. Since
the exact inference in the non-Gaussian state-space modelis intractable, we estimate the model parameters using theparticle stochastic approximation EM (PSAEM) algorithm(Lindsten, 2013; Svensson et al., 2014), which combines
the efficient conditional particle filter with ancestor sampling (CPF-AS) (Lindsten et al., 2014) with the stochastic
approximation EM (SAEM) algorithm (Delyon et al., 1999).
Interestingly, in the extreme case where the aggregation
factor k becomes larger and larger, we show that the observed time series will become independent and identically
distributed (i.i.d.), and we study to what extent the underlying time-delayed causal relations can be recovered from theinstantaneous dependence in the observed data.

2

EFFECT OF TEMPORALAGGREGATIONIn the linear case, Granger causal analysis Granger (1969)
can be done by fitting the following first-order VAR model(Sims, 1980):
xt = Axt−1 + et ,(1)
where xt = (xt,1 , xt,2 , ..., xt,n )| is the vector of the observed data, et = (et,1 , ..., et,n )| is the temporally andcontemporaneously independent noise process, and A isthe causal transition matrix containing temporal causal relations.
2.1

WITH A FINITE kGong et al. (2015) studied causal discovery from subsampled data. With subsampling, the observations x̃s1:T ,(x̃s1 , x̃s2 , ..., x̃sT ) = (x1 , x1+k , ..., x1+(T −1)k ) follow
x̃st = Ax1+(t−1)k−1 + e1+(t−1)k= A(Ax1+(t−1)k−2 + e1+(t−1)k−1 ) + e1+(t−1)k= ...
= Ak x̃st−1 +k−1
XAl e1+(t−1)k−l ,

(2)

l=0

which turns out to be a VAR model with temporally independent and contemporaneously dependent noise process.
They demonstrated that it is possible to identify the highresolution causal relation A from the low-resolution observations x̃s1:T if the noise terms are non-Gaussian.
In this paper, we are concerned with the temporally aggregated data x̃1:T , (x̃1 , x̃2 , ..., x̃T ), which are obtained bytaking the average (or sum) of every non-overlapping k
Pkpoints, i.e., x̃t = k1 i=1 xi+(t−1)k , wherexi+(t−1)k = Ak xi+(t−2)k +k−1
XAl ei+(t−1)k−l .

l=0Taking the average of the above equation over i =1, 2, . . . , k, we have
x̃t =

k
k
k−1

Ak X
1 XX lxi+(t−2)k +A ei+(t−1)k−l
k i=1
k i=1
l=0

k−1

= Ak x̃t−1 ++m1 X X n(
A )etk−m
k m=0 n=0

k−1
X k−1
X

(



An )e(t−1)k−m+1 ,(3)

Denote by ēt the error term above, i.e., ēt =Pk
1i=1 ei+(t−1)k+1 . Note that ēt has contemporaneousk
independent components. Sincem=1 n=mwhich is a vector autoregressive-moving-average (VARMA)
model with one autoregressive term and two movingaverage terms:

1(x1+tk − x1+(t−1)k ) → 0, as k → ∞,kwe havex̃t = Ak x̃t−1 + ~etx̃t = Ax̃t + ēt ,= Ak x̃t−1 + M0 t + M1 t−1 ,(4)where t = k1 [etk , etk−1 , . . . , e(t−1)k+1 ]| , M0 = [I, I +
Pk−1
Pk−1
A, . . . , n=0 An ], and M1 = [ n=1 An , . . . , Ak−1 , 0].
Here, I represents the n×n identity matrix, and 0 represents
the n × n zero matrix. We call (A, e, k) the representationof the k-th order aggregated time series x̃. Clearly A cannot
be recovered by simply fitting a VAR model on x̃t , as done
by Granger causal analysis. Even if we using VARMA
modeling, we are only guaranteed to identify Ak insteadof the original A. In Section 3, we will show under what
conditions can we identify the causal relation A at the causal
frequency from the aggregated time series x̃1:T .
2.2

WHEN k → ∞Interestingly, causal discovery from aggregated data witha large aggregation factor k seems to have a wide range ofapplications. For instance, in the stock market, the causal
influences between stocks take place very quickly (as indicated by the efficient market hypothesis), but we usuallywork with low-frequency data such as daily returns. The
daily return is the sum of high-frequency returns within
the same day. Discovering the causal interactions between
stocks from their daily returns then become a problem of
causal discovery from aggregated data with a large k.
When the aggregation factor k is very large, ~et becomes a
mixture of numerous independent components. Fortunately,
we can use a simple model to approximate the generatingprocess of the aggregated data. From (1), we havekX

xi+(t−1)k+1 = A

i=1

k
X

xi+(t−1)k +

i=1

k
X

ei+(t−1)k+1 ,

i=1

(5)

as k → ∞. This is a linear instantaneous causal modelfor the components of x̃t because the components of thetotal error term, ēt , are still contemporaneously independent.
When the error terms are non-Gaussian, it has the same form
as the Linear, Non-Gaussian Model (LiNG) (Lacerda et al.,2008); when the causal relations are further assumed to beacyclic, it follows the form of the Linear, Non-Gaussian
Acyclic Model (LiNGAM) (Shimizu et al., 2006). The
difference is that in LiNG or LiNGAM, the self-loop influences, Aii , are assumed to be zero. We will also investigate
the identifiably of A in this case in Section 3.

3

IDENTIFIABILITY OF CAUSALRELATIONS IN A

We investigate the identifiability of the high-resolution
causal transition matrix A from the aggregated time timeseries x̃1:T . In other words, suppose x̃ also admits anotherrepresentation (A0 , e0 , k), we aim to see whether it is always the case that A = A0 as the sample size T → ∞.
If the noise terms follow the Gaussian distribution, A is
usually not identifiable (Palm & Nijman, 1984). Recently,
it has been shown that A is identifiable from subsampled
time series if the noise terms are non-Gaussian (Gong et al.,2015). However, this does not give rise to the identifiabilityof A from aggregated time series – the latter is much moredifficult to see as the aggregated model described in (3) has
a more complicated structure. Here, we show that, in theexact model (3), A is identifiable from the aggregated data
under appropriate conditions; furthermore, as k → ∞, theapproximate model (5) holds, and A is partially identifiablefrom the aggregated data, but the identification procedure iscomputationally much more efficient.
First, we will show that Ak can be identified by fitting the
VARMA model (4). We make the following assumption.

that is,
k

1X
x̃t =
xi+(t−1)k+1
k i=1
=A

k
h1 X

k i=1

xi+(t−1)k+1 −A1. At least one of the τ -step (τ ≥ 2) delayed cross covariance matrices of x̃t , E[x̃t−1 x̃|t−τ ], is invertible.
i1(x1+tk − x1+(t−1)k ) +kk1Xei+(t−1)k+1 .
k i=1Since t is both temporally and contemporaneously independent, t and t−1 are independent of x̃t−τ , which implies
thatE[t x̃|t−τ ] = 0, E[t−1 x̃|t−τ ] = 0.

Multiplying both sides of (4) from the right side by x̃|t−τ
and taking the expectation, we haveE[x̃t x̃|t−τ ]= Ak E[x̃t−1 x̃|t−τ ] + M0 E[t x̃|t−τ ] + M1 E[t−1 x̃|t−τ ]= Ak E[x̃t−1 x̃|t−τ ].

(6)Under the assumption A1 , we can first see that Ak is identifiable:
Ak = E[x̃t x̃|t−τ ] · E[x̃t−1 x̃|t−τ ]−1 .
3.1

(7)

IDENTIFIABILITY WITH FINITE k

Substituting the above equation into (3), one can then find
~et , which is defined to be linear mixtures of (2k − 1) noise
terms, i.e., etk , etk−1 , . . ., and e(t−2)k+2 . In the following,
we will concentrate on the identifiablity of A from ~e.
Let

The following identifiability result on A states that in various situations, A for the original high frequency data isfully identifiable.
Theorem 1. Suppose all of eit are non-Gaussian, and thatthe data x̃t are generated by (3) and that it also admits
another kth order aggregation representation (A0 , e0 , k).
Let assumptions A1 and A2 hold. When the number ofobserved data points T → ∞, the following statements aretrue.
(i) A0 can be represented as A0 − I = (A − I)D, where
D is a diagonal matrix with 1 or −1 on its diagonal.
If we constrain all the self influences, represented bydiagonal entries of A and A0 , to be no greater than 1,
then A0 = A.
(ii) If each pei is asymmetric, we have A0 = A.
A complete proof of Theorem 1 can be found in Section 6.

H , [I, I + A, . . . ,

k−1
X

Al ,

l=0

k−1
X
l=1

Al , . . . ,k−1
XAl , Ak−1 ].

(9)

(0)
(1)
(2k−2)(1)ẽ = (e1 , ..., e(0), ..., e(2k−2))| .
n , e1 , ..., en , ..., e1n(l)

IDENTIFIABILITY AS k → ∞

l=k−2

(8)
The error terms in (3) correspond to the following mixing
procedure of random vectors:~e = Hẽ, where3.2

We have shown that A is identifiable from aggregated data(3) when k is finite. However, when k becomes larger,
estimating A will encounter more difficulty because more
independent components in (9) are involved. When k = ∞,it is not necessary for Proposition 1, as well as Theorem 1,
to hold, because ~e in (9) is the mixture of an infinite numberof independent components.

(l)Here, el = (e1 , ..., en )| together with the time index t
represent etk−l . The components of ẽ are independent, and(l)for each i, ei , l = 0, ..., 2k − 2, have the same distribution
pei . Under the condition that pei is non-Gaussian for eachi, H can be estimated up to the permutation and scaling
indeterminacies (including the sign indeterminacy) of thecolumns, as given in the following proposition.
Proposition 1. Suppose that all pei are non-Gaussian.
Given k and x̃1:T generated according to (3), H can be
determined up to permutation and scaling of columns.
For the proof of Proposition 1, please refer to (Gong et al.,2015).

Interestingly, as k → ∞, x̃t follows an instantaneous causal
model in the i.i.d. case, as shown in (5). We will thenanswer the following two questions. In this case, can westill estimate A from aggregated data? If we can, is there
an efficient procedure to do so?
Equation (5) implies (I − A)x̃t = ēt . That is, applying
the linear transformation (I − A) on x̃t produces independent components, as components of ēt . This can beachieved by the independent component analysis (ICA) procedure (Hyvärinen et al., 2001), and (I − A) can be estimated up to row scaling and permutation indeterminacies.
We then have the following observations.

We make the following assumptions on the underlying dynamic process (1) and the distributions pei , and then we
have the identifiability result for the causal transition matrix
A.

First, the diagonal entries of A, Aii , which represent theself influences or “self-loops" of the time-delayed causal
relations, cannot be determined (Lacerda et al., 2008). (Herewe have assumed Aii 6= 1.) This is because the scale ofeach row of (I − A) is unknown, and so is (1 − Aii ).

A2. The system is stable, in that all eigenvalues of A havemodulus smaller than one.

Let DA be the diagonal matrix with A11 , A22 , ..., Ann on
its diagonal. Equation (5) is equivalent to

A3. The distributions pei are different for different i after
re-scaling by any non-zero scale factor, their characteristic functions are all analytic (or they are all nonvanishing), and none of them has an exponent factor
with a polynomial of degree at least 2.

(I − DA )x̃t = (A − DA )x̃t + ēt⇒x̃t = (I − DA )−1 (A − DA ) x̃t + (I − DA )−1 ēt .
|{z},AN oSelf Loop

(10)

Secondly, suppose there is no feedback loop between
the processes after removing the self-loops, meaning that(A − DA ) can be permuted to a strictly lower-triangularmatrix by equal row and column permutations. According to the LiNGAM model, which assumes these is noself-loop, (I − DA )−1 (A − DA ) in (10) can be uniquely
estimated (Shimizu et al., 2006). In other words, if one
applies LiNGAM analysis on x̃t , the estimated causal coefficients from th ith variable to the jth variable is actually(1 − Ajj )−1 Aji . From this we can see whether Aji is zeroor not; furthermore, if the self-loops Ajj are given by priorknowledge, then A is fully identifiable.
Thirdly, suppose there exist feedback loops between the processes after removing the self-loops. In this case, (A − DA )cannot be permuted to a strictly lower-triangular matrix byequal row and column permutations. The identifiability of
A in (5) has been studied by Lacerda et al. (2008): suppose
the feedback loops are disjoint, although in theory there are
multiple solutions to AN oSelf Loop , the most stable solution(the product of the coefficients in AN oSelf Loop along eachloop is minimized) is unique.

4

ESTIMATING THE CAUSALRELATIONS FROM AGGREGATEDDATAIn this section, we present the algorithm to estimate A fromaggregated data with finite k. Clearly, the larger k, the moredifficult it is to estimate A from aggregated data. Therefore,
when k is relatively large (say, larger than 6), we advocate
the methods given in Section 3.2 to partially estimate A.
Since the identifiability of A from aggregated data relieson the non-Gaussianity of the error terms, we use Gaussian
mixtures to represent their distributions. It is natural to do
parameter estimation with the Expectation-Maximization
algorithm, which, unfortunately, involves a large number ofGaussian components. To avoid this issue, we propose to
use the Stochastic Approximation EM (SAEM) algorithm,
as a variant of EM, and further resort to conditional particle filtering with ancestor sampling (CPF-AS) to achieve
computational efficiency.

4.1

STATE-SPACE MODELING

We can consider (3) as a special state-space model:x̃t = Ak x̃t−1 + Hẽth P
1l
= Ak x̃t−1 + I
l=0 A

...

Pk−1

l=0 A

l


Pk−1

l=1 A

l

...

Pk−1

l=k−2 A

l

A

k−1

(0)

ẽt
(1)
ẽt
..
.




i



,




(2k−2)
ẽt
(11)

where

 
(0)
etk
ẽt−1
etk−1 
 (1)   ẽt−1  




..
 = F .  + ,.
 .  
 .  e(t−1)k+1 (2k−2)(2k−2)ẽt−1
ẽt0(nk−n)×1

0nk×(nk−n)0nk×(nk)F=,
I(nk−n)×(nk−n) 0(nk−n)×(nk)


(0)

ẽt
(1)
ẽt
..
.





(l)

ẽt = etk−l , and the noise terms etk , etk−1 , ..., e(t−1)k+1
share the same distribution for the same channel and are
mutually independent. Since the non-Gaussianity is essential to the identifiability of A, we use a Gaussian mixture modelof noise term e, i.e.,
Pmto represent each channel
2
p ei =
wi,c N (ei |µi,c , σi,c), where wi,c ≥ 0, and
c=1Pm
c=1 πi,c = 1, for i = 1, ..., n. Correspondingly, eachchannel of ẽ is also represented by a Gaussian mixture
model.
We aim to estimate the parameters A and the noise terms (ifnecessary) in the above state-space model. We introduce theadditional latent variable z̃t = (z̃t,1 , ..., z̃t,nk )| , in whichz̃t,j ∈ {1, . . . , m}, to model the distribution of noise termsẽt by Gaussian mixture models. The joint distribution of thestate-space model (11) over both observed and unobserved
variables is given byp(x̃1:T , ẽ1:T , z̃1:T ) = p(z̃1 )p(ẽ1 |z̃1 )p(x̃1 |ẽ1 )TYp(z̃t )p(ẽt |ẽt−1 , z̃t )p(x̃t |x̃t−1 , ẽt ).

(12)

t=2

The distributions in (12) are specified as follows:p(z̃t ) =nkYp(z̃t,j ) =nkYπ̃j,z̃t,j ,(13a)p(ẽt |ẽt−1 , z̃t ) = N (ẽt |Fẽt−1 + µ̃t , Σ̃t ),

(13b)

j=1

j=1

k

p(x̃t |ẽt , x̃t−1 ) = N (x̃t |A x̃t−1 + Hẽt , Λ).

(13c)

Since there are no additional additive noise terms in themodel, we fix Λ to a small value in our estimation algorithm
for regularization. µ̃t is the conditional mean of ẽt , i.e.,µ̃t = [µ̃1,z̃t,1 , . . . , µ̃nk,z̃t,nk , 01×n(k−1) ]| . Σ̃t is a diagonal matrix containing the conditional variance parameters22of ẽt , i.e., Σ̃t = diag([σ̃1,z̃, . . . , σ̃nk,z̃, 01×n(k−1) ]).
t,1t,nk
According to the structure of ẽ, the parameters π̃j,z̃t,j ,µ̃j,z̃t,j , and σ̃j,z̃t,j are controlled by the parameters of e,i.e., π̃i+nl,c = πi,c , µ̃i+nl,c = µi,c , and σ̃i+nl,c = σi,c , fori = 1, ..., n, l = 0, ..., k − 1, and c = 1, ..., m.

4.2

STOCHASTIC APPROXIMATION EM

4.3

The expectation-maximization (EM) algorithm is usuallyadopted to find the maximum likelihood estimation of the parameters in a probabilistic model with unobserved variables.
We can estimate the parameters θ = (A, wi,c , µi,c , σi,c )in (12) using the EM algorithm that iteratively maximizes the lower boundRof the marginal log-likelihoodPpθ (x̃1:T ) = log z̃1:T pθ (x̃1:T , ẽ1:T , z̃1:T )dẽ1:T . In
the E-step, at the k-th iteration, given the parameters θk−1 estimated from the (k − 1)-th iteration, theEM algorithm firstly computes the posterior distribution pθk−1 (z̃1:T , ẽ1:T |x̃P
1:T ) andR then computes the lower
bound Q(θ, θk−1 ) =pθk−1 (z̃1:T , ẽ1:T |x̃1:T ) logz̃1:Tpθ (x̃1:T , ẽ1:T , z̃1:T )dẽ1:T . In the M-step, the parameters are
updated as θk = arg maxθ Q(θ, θk−1 ).
However, we note that the number of Gaussian mixturesin the posterior distribution grows exponentially with thedimension of the time series, n, the number of aggregation factor, k, and the duration of time series T . Therefore, computing the exact posterior pθk−1 (z̃1:T , ẽ1:T |x̃1:T )and Q(θ, θk−1 ) is intractable in this situation. A possible
solution is to adopt the monte carlo EM (MCEM) algorithm (Wei & Tanner, 1990), which approximately calculates Q(θ, θk−1 ) using samples drawn from the posterior distribution pθk−1 (z̃1:T , ẽ1:T |x̃1:T ). However, MCEM makes
inefficient use of generated samples, as it discards samples generated in the previous EM iterations. Therefore, alarge number of sample points are required for each iteration, which is computationally expensive when the sampling
method is complex.

CONDITIONAL PARTICLE FILTER WITHANCESTOR SAMPLING

In our model, sampling from the posterior pθk−1 (z̃1:T |x̃1:T )is usually performed using forward filter/backward simulator particle smoother, which typically requires a large number of particles to generate a smooth backward trajectoryz̃1:T [k]. To reduce the number of required particles, we use
the Markovian version of SAEM (Kuhn & Lavielle, 2004),
which samples from a Markov kernel Mθt−1 , leaving theposterior distribution invariant. Specifically, let z̃1:T [k − 1]be the previous draw from the Markov kernel, the current
state is sampled by z̃1:T [k] ∼ Mθk−1 (·|z̃1:T [k − 1]). Following (Lindsten, 2013; Svensson et al., 2014), we construct
the Markov kernel using Rao-Blackwellized conditional particle filter with ancestor sampling (RB-CPF-AS) (Lindsten
et al., 2014), which was originally proposed for Gibbs sampling.
The machinery inside RB-CPF-AS resembles a standard
particle filter, with two main differences: one particle tra0
jectory is deterministically set to a reference trajectory z̃1:T ,and the ancestors of the reference trajectory are randomly
chosen and stored during the algorithm execution. Algorithm 1 gives a brief description of the RB-CPF-AS ali
gorithm. Let {z̃i1:t−1 , wt−1}N
i=1 be the approximation ofpθ (z̃1:t−1 |x̃1:t−1 ), RB-CPF-AS propagates this sample totime t by introducing the auxiliary variables {ait }N
i=1 , referred to as ancestor indices. To generate z̃it for the firstN − 1 particle trajectories, we first sample the ancestor inj
dex according to P (ait = j) ∝ wt−1
, and then sample z̃it,jQnk
ait
according to pθ (z̃t |z̃t−1 ) = j=1 pθ (z̃t,j ). The first N − 1ai

To reduce the number of simulated sample points, we propose to use the stochastic approximation EM (SAEM) algorithm (Delyon et al., 1999), which only requires of a single
realization of the unobserved variables at each iteration. At
the k-th iteration, the E-step and M-step are replaced by thefollowing:t, z̃it }. The
trajectories are then augmented as z̃i1:t = {z̃1:t−1
0
N -th particle is set to the reference particle, z̃Nt = z̃t , andthe ancestor index aNt is sampled according toE-step: Generate a single sample point z̃1:T [k] fromthe posterior pθk−1 (z̃1:T |x̃1:T ), and compute∝ pθ (x̃t:T , z̃t:T |x̃1:t−1 , z̃i1:t−1 )pθ (z̃i1:t−1 |x̃1:t−1 ), (15)0∝ pθ (z̃i1:t−1 |z̃t:T , x̃1:T )0iwhere pθ (z̃i1:t−1 |x̃1:t−1 ) = wt−1
andQ̂k (θ) = (1 − γk )Q̂k−1 (θ)
Z+ γk pθk−1 (ẽ1:T |z̃1:T [k], x̃1:T )log pθ (x̃1:T , ẽ1:T , z̃1:T [k])dẽ1:T .

P (aNt = i)

0(14)

M-step: Update parameters by θk = arg maxθ Q̂k (θ).
In (14), P{γk }∞of decreasing step sizes satk=1 is a sequenceP∞∞isfying k=1 γk = ∞ and k=1 γk2 < ∞. Here we use
Rao-Blackwellization (Svensson et al., 2014) to avoid sampling ẽ1:T because it is analytically integrable. It has beenshown in (Delyon et al., 1999) that the resulting sequence{θk }k≥1 will converge to a stationary point of pθ (x̃1:T )under weak assumptions.

pθ (x̃t:T , z̃t:T |x̃1:t−1 , z̃i1:t−1 )1 i∝ |Mit−1 |−1/2 exp(− ηt−1
).
2

(16)

Conditioned on z̃i1:T , we can calculate p(ẽ1:T |z̃i1:T , x̃1:T )using Kalman filter and Rauch-Tung-Striebel (RTS)
smoother. The filtering, prediction, and smoothing PDFs
areipθ (ẽt |z̃i1:t , x̃1:t ) = N (ẽt |µ̂if,t , Σ̂f,t ),ipθ (ẽt+1 |z̃i1:t , x̃1:t ) = N (ẽt+1 |µ̂ip,t , Σ̂p,t ),ipθ (ẽt |z̃i1:T , x̃1:T ) = N (ẽt |µ̂is,t , Σ̂s,t ),

(17a)
(17b)
(17c)

i

respectively. Let kek2Ω = e| Ωe and Σ̂f,t = Γif,t Γi,|f,t . In(16),i
Mit = Γi,|f,t Ωt Γf,t + I,

(18a)

i
2
ηti = kµ̂if,t k2Ωt − 2λ|t µ̂if,t − kΓi,|f,t (λt − Ωt µ̂f,t )kM−1 ,
t

(18b)
where
0

0

1/2

1/2

Ωt = F| (I − Ω̂t+1 Σ̃t+1 Υ−1
t+1 Σ̃t+1 )Ω̂t+1 F,

(19a)

0

0
0
1/2
1/2
λt = F| (I − Ω̂t+1 Σ̃t+1 Υ−1
t+1 Σ̃t+1 )(λ̂t+1 − Ω̂t+1 µ̃t+1 ),

(19b)
|

Ω̂t = Ωt + H Λ−1−1

λ̂t = λt + H| Λ
0

1/2

Υt = Σ̃t

0

H,

(19c)

(x̃t − Ak x̃t−1 ),

(19d)

1/2

Ω̂t Σ̃t+ I,

0Input: z̃1:T = z̃1:T [k − 1], θ = θk−1
Output: z̃1:T [k] ∼ Mθk−1 (·|z̃1:T [k − 1])
Compute {Ωt , λt }Tt=1 according to (19a)-(19g)
Draw z̃i1 with z̃i1 ∼ pθ (z̃1 ) for i = 1, . . . , N
Compute µi1 , Σi1 , and w1i for i = 1, . . . , N
for t=2 to T doj
Draw ait with P (ait = j) ∝ wt−1
for i = 1, . . . , N// Resampling and ancestor sampling
Draw z̃it with z̃it ∼ pθ (z̃t ) for i = 1, . . . , Ni} according to (18a) and (18b)Compute {Mit−1 , ηt−1iDraw aN according to (15) for i = 1, . . . , N// Particle propagationait
Set z̃i1:t = {z̃1:t−1, z̃it } for i = 1, . . . , N
ai

0

t,nk

, 01×n(k−1) ]| ,

i

t,1

t,nk

(19f)

, 01×n(k−1) ]).

(19g)

With ΩT = 0 and λT = 0, {Ωt , λt }Tt=1 can be computed
0recursively for z̃1:T using (19a)-(19g). Once all the ancesi N
tors {at }i=1 have been sampled, we can calculate the newparticle weights as follows
wti ∝ pθ (x̃t |z̃i1:t , x̃1:t−1 )i= N (x̃t |Fµ̂ip,t+1 + Ak x̃t−1 , FΣ̂p,t+1 F| + Λ).

(20)

After all the particle trajectories have been generated, we
obtain z̃1:T [k] by sampling from these trajectories according
to the weights {WTi }N
i=1 at time T .
4.4

T
X

log p(z̃t )
T

1X
−
2 t=1
T

1X
2 t=1

Z

Z

0

0

0

0

−1

q(ẽt )(ẽt − µ̃t )| Σ̃t

0

0

0

(ẽt − µ̃t )dẽt + log |Σ̃t |q(ẽt )(ỹt − Hẽt )| Λ−1 (ỹt − Hẽt )dẽt+ const,(21)
0

where z̃t = z̃t [k], ẽt = [e|tk , e|tk−1 , . . . , ẽ|(t−1)k+1 ]| ,
0

µ̃t

ait

Draw J with P (J = j) ∝ WTj and set z̃1:T [k] = z̃J1:T .

2, . . . , σ̃nk,z̃]), ỹt = x̃t − Ak x̃t−1 , and q(ẽt ) =t,nkpθk−1 (ẽt |z̃1:T , x̃1:T ).

It
R can be seen that
R we only need sufficient statisticsẽt q(ẽt )dẽt and ẽt ẽ|t q(ẽt )dẽt to maximize (21). Denoting a sufficient statistics at the k-th iteration as S k , weuse Sk = (1 − γk )Sk + γk S k to maximize Q̂k (θ). To maximize Q̂k (θ) with respect to A, we compute the gradientof A in terms involving Ak and H and apply a conjugate
gradient descent method as done in Gong et al. (2015).

5

0

=

[µ̃1,z̃t,1 , . . . , µ̃nk,z̃t,nk ]| , Σ̃t=2diag([σ̃1,z̃t,1EXPERIMENTS

In this section, we conduct empirical studies of the two
estimation methods presented in Section 3.2 and Section 4
on both synthetic and real data to show their effectiveness.
5.1

t=1

−

i

PARAMETER UPDATEAt the k-th M step, given the sample z̃1:T [k] drawn by
RB-CPF-AS, we can obtain pθk−1 (ẽt |z̃1:T [k], x̃1:T ) =N (ẽt |µ̂s,t , Σ̂s,t ) using the RTS smoother. Then we haveZpθk−1 (ẽ1:T |z̃1:T , x̃1:T ) log pθ (x̃1:T , ẽ1:T , z̃1:T )dẽ1:T=

ait

Σ̂f,1:t−1 = Σ̂f,t−1 , Σ̂p,1:t−1 = Σ̂p,1:t−1// WeightingiiCompute µ̂if,t , µ̂ip,t , Σ̂f,t , and Σ̂p,tCompute weights wti according to (20)022Σ̃t = diag([σ̃1,z̃
0 , . . . , σ̃
0
nk,z̃

ai

t
t
Set µ̂if,1:t−1 = µ̂f,1:t−1, µ̂ip,1:t−1 = µ̂p,1:t−1
,

(19e)

µ̃t = [µ̃1,z̃0 , . . . , µ̃nk,z̃0t,1Algorithm 1: RB-CPF-AS

SIMULATED DATA

We conduct a series of simulations to investigate the effectiveness of the proposed estimation methods. Following(Gong et al., 2015), we generated the data at the casualfrequency using the VAR model (1) with randomly generated matrix A and independent Gaussian mixture noises
et . The elements in A were drawn from a uniform distribution U(−0.5, 0.5). The Gaussian mixture model contains
two components for each channel. The parameters were
w1,1 = 0.2, w1,2 = 0.8, w2,1 = 0.3, w2,2 = 0.7 ,µi,1 = 0,222µi,2 = 0, σi,1= 1e − 4, σ1,2
= .1, and σ2,2
= 0.2. Lowresolution observations were obtained by aggregating thehigh-resolution data using aggregation factor k. Similarly,

we also generated data with Gaussian noise (by setting
22
σ1,c= 0.01, σ2,c
= 0.02) for comparison of different
methods. We tested data with dimension n = 2, aggregation factor k = 2 and 3, and sample size T = 150 and
300, respectively. For comparison,we replaced the Gaussian
mixture models in our method with Gaussian noise models,leading to a method based on Gaussian noises. We denotethe method proposed in Section 4 that performs causal discovery from temporally aggregated data as CDTAfinite andthe corresponding Gaussian counterpart as CDTAGauss . We
also compare with the NG-EM method (Gong et al., 2015)
on the aggregated data with non-Gaussian noises. The experiments are repeated for 10 replications .
Table 1 shows the mean squared error (MSE) of the estimated causal transition matrix A. It can be seen that as the
sample size T increases, both the proposed CDTAfinite andthe baseline method CDTAGauss obtain smaller estimation errors. On the non-Gaussian data, CDTAGauss produces muchhigher errors than CDTAfinite . On the Gaussian data, neither
CDTAfinite nor CDTAGauss can obtain accurate estimations.
This is because the estimation algorithms can converge tomany possible solutions that have the same marginal likelihood, if the data noises are Gaussian or the estimation
algorithms assume a Gaussian noise model. The results
are consistent with the theoretical results that the causal
relations might not be uniquely determined using Gaussian
noise models. It can also be seen that the NG-EM method
fails on the aggregated data, because NG-EM is proposed
for subsampled rather than aggregated data.
Further, we examined the performance of the method described in Section 3.2, denoted as CDTAinfty , with finite
k values. To achieve so, we generated aggregated data
0.1
with A = [ 0.8
0.6 0.7 ], aggregation factor k = 2, 3, 4, 10,
and the same Gaussian mixture noise parameters describedabove. The true ANoSelfLoop in this case can be calculated as ANoSelfLoop = [ 02 0.5
0 ]. Using the linear instantaneous non-Gaussian model, we can obtain the estimations of ANoSelfLoop on the aggregated data. The results fork = 2, 3, 4, 10 are given as follows:
NoSelfLoop
00.52
0 0.47 ],ÂNoSelfLoop= [ −13.5= [ 1.430 ], Â3
2
0
0 0.44 ], ÂNoSelfLoop = [ 0 0.51 ]. (22)ÂNoSelfLoop= [ 1.61
4
10
0
2.03 0It seems that when k ≥ 4, the linear instantaneous nonGaussian causal model (10), which assumes that there isno self-loop, can estimate the corresponding ANoSelfLoop
accurately and very efficiently, at the cost of losing the selfloops in the original process. However, the self-loops canbe estimated with CDTAfinite when k is reasonably large. As
a cautionary notice, researchers should carefully interpretthe estimated parameters produced by linear instantaneous
causal models, which assume there is no self-loop; as aconsequence, the linear instantaneous non-Gaussian causal
model produces (1 − Ajj )−1 Aji , whose magnitude can bevery different from the true causal parameter Aji .

5.2

REAL DATA

We conducted experiments on the Temperature Ozone data
(Mooij et al., 2016) and the macroeconomic data used in(Moneta, 2008). These two time series are collected byaveraging the records during specified time intervals. For
example, the Temperature Ozone data contain daily meanvalues of ozone and temperature of year 2009 in Chaumontand Switzerland. The macroeconomic data contain quarterly
US macro variables for the period 1947:2 to 1994:1.
Temperature Ozone. The Temperature Ozone datais the 50th causal-effect pair from the website https://webdav.tuebingen.mpg.de/cause-effect/.
The data have records of ozone density X and dailymean temperatureY . The ground truth is causal relation is
Y → X. We first applied CDTAinfty on the data, resulting
0 0.65 ]. From this result, we can findin ANoSelfLoop = [ 0.65
0that instantaneous effects exist in both directions. This
could possibly caused by aggregation with a small k as theestimated ANoSelfLoop is likely to be inaccurate. We thenestimated the causal matrix by CDTAfinite . The estimated0.0945transition matrix A for k = 1, 2, 3 is [ 0.8418
0.0165 0.9785 ],0.84260.11510.83600.1150[ 0.0220 0.9702 ], [ 0.0172 0.9729 ], respectively. It seems thatthe estimated matrices sensibly captured the self-influences
and cross-influences between the ozone and temperatureprocesses.
Macroeconomic Data The data are quarterly U.S. observational on real aggregated macroeconomic variables. Herewe consider the causal relations between two variables, including real balances X and price inflation Y . X denotes
the logarithm of per captita M2 minus the logarithm of theimplicit price deflator. Y is the log of the implicit pricedeflator at the time t minus log of the implicit price deflator
at the time t − 1. Again, we first applied the CDTAinfty
to find the rough estimation of causal relations excluding
00
self-loops. The estimated ANoSelfLoop is [ 0.3074
0 ]. This indicates that no influence from effect to cause can be estimatedfrom the instantaneous dependencies, which is consistent
with the ground truth. We also employed CDTAfinite to obtain the estimation of a complete causal transition matrix
A. The estimated transition matrix A for k = 1, 2, 3, 4 is1.009 −0.0013
0.9007 −0.0011
0.8503 −0.0005[ 0.1718
0.5796 ], [ 0.2024 0.6420 ], [ 0.1378 0.7091 ], respectively. We can see that A gives weaker responses from
effect Y to cause X as k increases. If we consider k = 4 as
the aggregation factor, then we can calculate from the esti00.033 ], which is close tomated A that ANoSelfLoop = [ 0.4737
0
the results estimated by CDTAinfty .

6

CONCLUSION

In this paper, we have investigate the problem of discoveringhigh frequency causal relations from temporally aggregatedtime series. When the aggregation factor is finite, we proved

non-Gaussian noise

Data
Methodsk=2
T=150
2.10e-4
1.28e-2
8.75e-2

CDTAfinite
CDTAGauss
NG-EM

Gaussian noise
k=3

T=300
1.19e-44.49e-3
8.51e-2

T=150
8.17e-4
1.20e-2
5.27e-1

k=2
T=300
7.36e-4
7.22e-3
1.88e-1

T=150
1.42e-2
1.13e-2
-

k=3
T=300
3.67e-3
3.08e-3-

T=150
7.63e-3
6.26e-2-

T=300
9.69e-3
9.07e-3-

Table 1: MSE of different methods on simulated non-Gaussian and Gaussian data. The results are shown for subsampling
factors (k = 2, 3) and lengths of data (T = 150, 300).
0
0
0
Hkn+j
= Hm+j+ Hm+k+j. Hence, Hkn+i is propor00
tional to Hkn+j , i.e., Hkn+j= λki Hkn+i . Assume thatP1
Pk−1 0l Pk−1 0l
Pk−1
0l
0l
l=0 A , . . .,l=0 A ,
l=1 A , . . .,l=k−2 A ,0k−1
0A
are non-diagonal matrices, we have Hj = λ0i Hi .
Since Hi and Hj0 must be columns of I, as implied by
the structure of H and H0 , we can see that λ0i = 1 andthat i = j. Consequently, λki must be 1 or -1. Let
B = I + A + . . . + Ak−1 and B0 = I + A0 + . . . + A0k−1 ,we thus have B = B0 D, where D is a diagonal matrix
with 1 or -1 as its diagonal entries. Moreover, because
AB − B = Ak − I and A0 B0 − B0 = A0k − I, and
Ak = A0k , we havethat the causal relations are fully identifiable if the underlying causal relations are linear and the noise process is
non-Gaussian. We also show that the causal matrix thatremoves self-loops is identifiable from instantaneous dependencies, when the aggregation factor goes to infinity.
Based on these results, we propose an algorithm to recover
the complete causal matrix when the aggregation factor isrelatively small an an very efficient algorithm to partially
recover the matrix when the aggregation factor is relatively
large. Future work will focus on automatically estimating
the aggregation factor k from data.

APPENDIX: PROOF OF THEOREM 1

A0 − I = (A − I)D.

If both A0 and A have diagonal entries which are smaller
than 1, D must be the identity matrix, i.e., A0 = A. Therefore statement (i) is true.

Proof. Here we consider the limit when T → ∞. According to the identifiability results of Ak (7), we have
Ak = A0k .

(23)

If each pei is asymmetric, ei and −ei have different distributions. Consequently, the representation (24) does not holdany more if one changes the signs of a subset of, but not
0all, non-zero elements of {Hln+j| l = 0, ..., 2k − 2}. Thisimplies that for non-zero Hln+i , λli , including λ0i , have thesame sign, and they are therefore 1 since λ0i = 1. λki = 1
leads to D = I and thus gives A0 = A. That is, (ii) is true.

We then consider the remaining error term ~et . The corresponding random vector ~e follows both the representation(9) andH0 , [I, . . . ,

k−1
X
l=0

A0l ,

k−1
X

A0l , . . . ,

l=1

k−1
X

(26)

A0l , A0k−1 ].

l=k−2

(24)
(0)

(0)

(1)

(1)

(2k−2)ẽ0 = (e0 1 , ..., e0 n , e0 1 , ..., e0 n , ..., e0 1(2k−2) |, ..., e0 n

(25)
(l)
with e0 i , l = 0, ..., 2k − 2, having the same distribution
pe0i .
According to Proposition 1, each column of H0 is a scaled
version of a column of H. Denote by Hln+i , l = 0, ..., 2k −2; i = 1, ..., n, the (ln + i)th column of H, and similarly for0
Hln+i. According to the Uniqueness Theorem in Eriksson &Koivunen (2004), we know that under condition A2, for eachi, there exists one and only one j such that the distribution(l)of ei , l = 0, ..., 2k − 2 (which have the same distribution),(l)is the same as the distribution of e0 j , l = 0, ..., 2k − 2,
up to changes of location and scale. As a consequence,0
the columns {Hln+j| l = 0, ..., 2k − 2} correspond to{Hln+i | l = 0, ..., 2k − 2} up to the permutation andscaling arbitrariness.
According to the structure of H, ∀ m ≤ k − 1, Hkn+i =Hm+i + Hm+k+i , and similarly we have ∀ m ≤ k − 1,)

Acknowledgements
The authors would like to thank Dr. Tongliang Liu for helpful discussions. DT and MG would like to acknowledge
the support from DP-140102164, FT-130101457, and LP150100671. CG and KZ would like to acknowledge the support from NIH-1R01EB022858-01 FAIN-R01EB022858,
NIH-1R01LM012087, and NIH-5U54HG008540-02 FAINU54HG008540.

References
Boot, J.CG., Feibes, W., and Lisman, J. Hubertus C. Further
methods of derivation of quarterly figures from annual data.
Applied Statistics, pp. 65–75, 1967.
Breitung, J. and Swanson, N. R. Temporal aggregation and spurious instantaneous causality in multiple time series models.
Journal of Time Series Analysis, 23:651–665, 2002.
Danks, D. and Plis, S. Learning causal structure from undersampled time series. In JMLR: Workshop and Conference Proceedings, 2013.

Delyon, B., Lavielle, M., and Moulines, E. Convergence of astochastic approximation version of the em algorithm. Annalsof statistics, pp. 94–128, 1999.
Eriksson, J. and Koivunen, V. Identifiability, separability, anduniqueness of linear ICA models. IEEE Signal Processing
Letters, 11(7):601–604, 2004.
Geiger, P., Zhang, K., Gong, M., Schölkopf, B., and Janzing,
D. Causal inference by identification of vector autoregressive processes with hidden components. In 32nd International
Conference on Machine Learning, pp. 1917–1925. Microtome
Publishing, 2015.
Ghysels, E., Hill, J.B., and Motegi, K. Testing for granger causality
with mixed frequency data. Journal of Econometrics, 192(1):207–230, 2016.
Gong, M., Zhang, K., Schölkopf, B., Tao, D., and Geiger, P. Discovering temporal causal relations from subsampled data. In
ICML, pp. 1898–1906, 2015.
Granger, Clive WJ. Investigating causal relations by econometricmodels and cross-spectral methods. Econometrica: Journal of
the Econometric Society, pp. 424–438, 1969.
Granger, Clive WJ. Implications of aggregation with common
factors. Econometric Theory, 3(02):208–222, 1987.
Harvey, A. C. and Chung, C. H. Estimating the underlying change
in unemployment in the uk. Journal of the Royal Statistics
Society, Series A, 163:303–309, 2000.
Hyttinen, A., Plis, S., Järvisalo, M., Eberhardt, F., and Danks, D.
Causal discovery from subsampled time series data by constraint
optimization. In International Conference on Probabilistic
Graphical Models, pp. 216–227, 2016.
Hyvärinen, A., Karhunen, J., and Oja, E. Independent Component
Analysis. John Wiley & Sons, Inc, 2001.
Kuhn, E. and Lavielle, M. Coupling a stochastic approximation
version of em with an mcmc procedure. ESAIM: Probability
and Statistics, 8:115–131, 2004.
Lacerda, G., Spirtes, P., Ramsey, J., and Hoyer, P. O. Discovering
cyclic causal models by independent components analysis. In
Proceedings of the 24th Conference on Uncertainty in Artificial
Intelligence (UAI2008), Helsinki, Finland, 2008.
Lindsten, F. An efficient stochastic approximation em algorithm
using conditional particle filters. In Acoustics, Speech and Signal Processing (ICASSP), 2013 IEEE International Conferenceon, pp. 6274–6278. IEEE, 2013.
Lindsten, F., Jordan, M. I, and Schön, T.B. Particle gibbs withancestor sampling. Journal of Machine Learning Research, 15(1):2145–2184, 2014.
Marcellino, M. Some consequences of temporal aggregation in
empirical analysis. Journal of Business and Economic Statistics,
17:129–136, 1999.
Moauro, F. and Savio, G. Temporal disaggregation using multivariate structural time series models. Journal of Econometrics,
8:210–234, 2005.
Moneta, A. Graphical causal models and vars: an empirical assessment of the real business cycles hypothesis. Empirical
Economics, 35(2):275–300, 2008.
Mooij, J. M, Peters, J., Janzing, D., Zscheischler, J., and Schölkopf,
B. Distinguishing cause from effect using observational data:
methods and benchmarks. Journal of Machine Learning Research, 17(32):1–102, 2016.
Palm, F. C. and Nijman, T. E. Missing observations in the dynamic
regression model. Econometrica, 52:1415–1435, 1984.

Pearl, J. Causality: Models, Reasoning, and Inference. Cambridge
University Press, Cambridge, 2000.
Plis, S., Danks, D., Freeman, C., and Calhoun, V. Rate-agnostic(causal) structure learning. In Advances in neural information
processing systems, pp. 3303–3311, 2015a.
Plis, S., Danks, D., and Yang, J. Mesochronal structure learning. In
Proceedings of the 31st Conference on Uncertainty in Artificial
Intelligence, volume 31. NIH Public Access, 2015b.
Proietti, T. Temporal disaggregation by state space methods: Dynamic regression methods revisited. The Econometrics Journal,
9:357–372, 2006.
Rajaguru, G. and Abeysinghe, T. Temporal aggregation, cointegration and causality inference. Economics Letters, 101:223–226,
2008.
Shimizu, S., Hoyer, P.O., Hyvärinen, A., and Kerminen, A.J. A
linear non-Gaussian acyclic model for causal discovery. Journalof Machine Learning Research, 7:2003–2030, 2006.
Silvestrini, A. and Veredas, D. Temporal aggregation of univariate
and multivariate time series models: A survey. Journal of
Economic Surveys, 22:458–497, 2008.
Sims, C. A. Macroeconomics and reality. Econometrica, 48:1–48,
1980.
Spirtes, P., Glymour, C., and Scheines, R. Causation, Prediction,and Search. MIT Press, Cambridge, MA, 2nd edition, 2001.
Stram, Daniel O and Wei, William WS. A methodological note on
the disaggregation of time series totals. Journal of Time Series
Analysis, 7(4):293–302, 1986.
Svensson, A., Schön, T.B., and Lindsten, F. Identification of jump
markov linear models using particle filters. In Decision andControl (CDC), 2014 IEEE 53rd Annual Conference on, pp.
6504–6509. IEEE, 2014.
Tank, A., Fox, E.B., and Shojaie, A. Identifiability and estimationof structural vector autoregressive models for subsampled andmixed frequency time series. arXiv preprint arXiv:1704.02519,
2017.
Tiao, George C. Asymptotic behaviour of temporal aggregates oftime series. Biometrika, pp. 525–531, 1972.
Van Nes, E.H., Scheffer, M., Brovkin, V., Lenton, T.M., Ye, H.,
Deyle, E., and Sugihara, G. Causal feedbacks in climate change.
Nature Climate Change, 5(5):445–448, 2015.
Wei, G.CG. and Tanner, M.A. A monte carlo implementationof the em algorithm and the poor man’s data augmentation
algorithms. Journal of the American statistical Association, 85(411):699–704, 1990.
Weiss, A. Systematic sampling and temporal aggregation in timeseries models. Journal of Econometrics, 26:271–281, 1984.
Zhou, D., Zhang, Y., Xiao, Y., and Cai, D. Analysis of sampling artifacts on the granger causality analysis for topology extractionof neuronal dynamics. Frontiers in computational neuroscience,
8, 2014.