# Graphical estimation of multivariate count time series

Vurukonda Sathish<EMAIL_ADDRESS>
Department of Electrical Engineering
Indian Institute of Technology Bombay,
Powai, 400076, India Debraj Chakraborty<EMAIL_ADDRESS>
Department of Electrical Engineering
Indian Institute of Technology Bombay,
Powai, 400076, India Siuli Mukhopadhyay<EMAIL_ADDRESS>
Department of Mathematics
Indian Institute of Technology Bombay,
Powai, 400076, India

###### Abstract

The problems of selecting partial correlation and causality graphs for count
data are considered. A parameter driven generalized linear model is used to
describe the observed multivariate time series of counts. Partial correlation
and causality graphs corresponding to this model explain the dependencies
between each time series of the multivariate count data. In order to estimate
these graphs with tunable sparsity, an appropriate likelihood function
maximization is regularized with an $\ell_{1}$-type constraint. A novel MCEM
algorithm is proposed to iteratively solve this regularized MLE. Asymptotic
convergence results are proved for the sequence generated by the proposed MCEM
algorithm with $\ell_{1}$-type regularization. The algorithm is first
successfully tested on simulated data. Thereafter, it is applied to observed
weekly dengue disease counts from each ward of Greater Mumbai city. The
interdependence of various wards in the proliferation of the disease is
characterized by the edges of the inferred partial correlation graph. On the
other hand, the relative roles of various wards as sources and sinks of dengue
spread is quantified by the number and weights of the directed edges
originating from and incident upon each ward. From these estimated graphs, it
is observed that some special wards act as epicentres of dengue spread even
though their disease counts are relatively low.

Keywords: graphical models, count data, partial correlation, Monte Carlo
expectation and maximization, $\ell_{1}$-type regularization.

##  1 Introduction

The dependencies between multiple interdependent time series data, such as the
temperatures recorded from neighbouring geographical regions (Bach and Jordan,
2004), stock prices from related markets (David A. Bessler, 2003), biological
signals that measure activities in the human brain (Qian et al., 2015) etc,
are studied using graphical models. In a graphical model, each data series is
represented as a vertex of a graph and the dependencies between the different
time series are represented by edges in the graph. In many applications such
as finance, insurance, biomedical, public health, etc, time series data is
frequently measured in the form of counts (Mehdi Jalalpour and Levin, 2015;
R.K.Freeland and McCabe, 2004). In this paper, we propose a novel algorithm
for inferring the dependencies and the strength of the dependencies (i.e., the
graphical model) for multiple inter-related count time series data. The
algorithm is applied to discover the spreading pattern of dengue infections in
a large Indian city.

In graphical models, the dependencies between multiple time series are
represented either by undirected or directed graphs (e.g., see W. J. Granger,
1969; Eichler, 2012; Dahlhaus and Eichler, 2003; Jung et al., 2015 and
references therein). First we review the undirected graphical models in brief
and thereafter review the directed graphical models. Graphical models were
introduced initially for multiple random variables in (Dempster, 1972) using
the inverse covariance matrix of the associated random vectors. (Brillinger,
1996) introduced such models for multivariate time series, where the
dependencies between the nodes were represented using partial correlations. In
(Songsiri et al., 2009) and (Dahlhaus, 2000), graphical description of vector
autoregressive (AR) models were considered for stationary time series and
defined in terms of the inverse spectral density matrix. These results were
extended to graphical models with vector autoregressive moving average (ARMA)
models in (Avventi et al., 2013). In (Songsiri and Vandenberghe, 2010), a
method to control the sparsity of the estimated graph was introduced through
an $\ell_{1}$-type regularization. Such a regularization (often called Least
Absolute Shrinkage and Selection Operator (LASSO)) on the standard maximum
likelihood estimation (MLE) method, reduces overfitting and provides control
over the sparsity of the estimated graphical model (Tibshirani, 1996).

Directed graphical models were introduced in (Wright, 1921, 1934) as path
diagrams for the study of hereditary properties, linking parents and children
graphically by arrows. For multivariate time series, one series is said to be
causal for another series if the prediction of second series using all
available information except the first series, can be improved by adding the
available information about the first series. This definition of the so-called
Granger causality was introduced in (W. J. Granger, 1969). Using this notion,
the causalities between multiple time series were represented as directed
Granger causality graphs, where each time series was represented as a vertex
and the causality between two different time series was represented by a
directed edge (Sims, 1972; Pierce and Haugh, 1977). The multivariate
autoregressive (AR) model of second order stationary time series was used to
obtain Granger causality graph (TjØstheim, 1981; Hsiao, 1982) where the
causalities between multiple time series were calculated directly from the
autoregressive coefficients. The Granger causality graphs were studied by
(Eichler, 2006, 2007) through multivariate autoregressive (AR) models of
second order stationary time series. Later, sparse Granger causality graphs
were discussed in (Songsiri, 2013, 2015; Zorzi and Sepulchre, 2016; Alpago et
al., 2018; Ciccone et al., 2018) through multivariate AR time series.

In all the above papers, the graphical models were discussed for time series
of Gaussian data. In contrast to the graphical models for Gaussian data,
literature on graphical models for count data is still sparse. The undirected
graphical models were introduced for multivariate count data in (Park and
Park, 2019) using a Bayesian approach to inference. Later, the undirected
graphical models were considered for multivariate count data in (Roy and
Dunson, 2020) to obtain the complex interactions between genes using a pseudo-
likelihood based algorithm. Unlike these directed models, Poisson directed
acyclic graphical (DAG) models, also referred to as Bayesian Networks (Park
and Raskutti, 2015, 2017), were obtained for multivariate count data in (Hue
and Chiogna, 2021). In all these papers, the graphical models are obtained for
multiple random variables. In contrast, our aim in this paper is to estimate
undirected and directed graphical model for multivariate count time series
where the dependencies between these time series are represented using
causalities and partial correlations.

To model count time series, generalized linear models (GLM) were discussed in
(Zeger and Qaqish, 1988; Li, 1994; Fahrmeir and Tutz, 2001; Fahrmeir and
Wagenpfeil, 1997; Durbin and Koopman., 2000). In (Utazi et al., 2018),
parameter driven models were used for modeling multivariate time series of
counts, using a latent process to account for the correlation between the
observed counts. In this paper, we choose this model for our study. In this
model, the latent process is unknown. While the unknown parameters in such
models can, in principle, be estimated using the well known expectation and
maximization (EM) algorithm (Dempster et al., 1977), associated difficulties
in calculating the conditional density functions requires us to use the the
Monte Carlo expectation and maximization (MCEM) algorithm (Wei and Tanner,
1990; Chan and Ledolter, 1995). In this paper, we extend a combination of the
methods from (Dempster et al., 1977; Wei and Tanner, 1990; Chan and Ledolter,
1995) to infer graphical models with tunable sparsity, using LASSO.

First we introduce certain choices about the model, formulation and
methodology. Following (W. J. Granger, 1969; Songsiri and Vandenberghe, 2010),
the dependencies between the observed multivariate time series of counts are
represented using partial correlations and causalities. However partial
correlations and causalities between multiple time series can only be
calculated easily for second order stationary time series. In our case, the
observed count data is usually non stationary, e.g., the data might be
seasonal in nature or follow various trends (Fahrmeir and Tutz, 1994). To
address this issue, we introduce an unobserved multivariate second order
stationary latent process to calculate the partial correlation and causality
between the observed multiple count time series. We define the partial
correlation and causality between elements of the observed time series in
terms of the partial correlation and causality between elements of the
stationary latent processes. As in (Chan and Ledolter, 1995), the observed
count data given the latent process, follows a Poisson distribution, while the
multivariate latent process follows an autoregressive (AR) Gaussian model.

It is known that the partial correlation between multiple time series can be
computed from the inverse spectral density matrix of a multivariate stationary
time series (Brillinger, 1981). Thus the partial correlations for the
stationary latent process in our model are obtained in terms of the parameters
of the latent AR model (Brockwell and Davis, 2002) using the inverse spectral
density matrix. Further, to obtain the partial correlations and causalities
for the observed multivariable count data from the stationary latent process,
we need to estimate unknown parameters of the AR model. Since the number of
possible edges in the partial correlation graph can be large, overfitting can
be a potential issue in this estimation problem (Cawley and Talbot, 2010).
Thus a maximum likelihood estimation (MLE), with $\ell_{1}$-type
regularization term involving the partial correlation constraints, is
formulated to get the desired sparsity in the estimated partial correlation
and causality graph. Since the latent process is unknown, a direct
maximization of the log-likelihood function with $\ell_{1}$-type
regularization is not possible. Following (Chan and Ledolter, 1995), an MCEM
algorithm with $\ell_{1}$-type regularization is proposed to estimate the
model parameters. In this algorithm, samples from the conditional density
function of the latent process given observed data, are generated using the
Metropolis-Hastings algorithm (Chan and Ledolter, 1995).

First, the proposed MCEM algorithm with $\ell_{1}$-type regularization is
tested on a simulated multivariable count dataset. The regularization
parameter is estimated using the trade-off curve (Songsiri, 2010) between log-
likelihood and $\ell_{1}$-type regularization function. We collect all the
graphs along the trade-off curve for different values of regularization
parameter and then choose the partial correlation and causality graph based on
the Bayes information criterion (BIC)(Songsiri, 2010).

Next, we use the developed algorithm to estimate the partial correlation and
causality graphs for the observed dengue disease data. The dengue data under
investigation consist of weekly dengue counts in a six year period, January
$2010$ to December $2015$, collected from each ward of Greater Mumbai city,
India. From these estimated partial correlation and causality graphs, we infer
the number of undirected, incoming and outgoing edges for each ward. These
edges and their weights give a quantitative measure of the interdependence and
directionality of the disease prevalence and development in the various wards.
Surprisingly, we observe that some wards act as the epicentres of disease
spread even though their absolute disease counts are relatively low.

This manuscript is an extended version of (Sathish et al., 2019) which was
presented in IEEE 58th Conference on Decision and Control (CDC), December
2019. It differs from (Sathish et al., 2019) as follows:

  1. 1.

In this manuscript, in addition to the problem of estimating partial
correlation graphs studied in (Sathish et al., 2019), the problem of
estimating causality graphs is also considered. Causality graphs were not
studied in (Sathish et al., 2019).

  2. 2.

No real world application was presented in (Sathish et al., 2019). Here, the
developed algorithms are applied on observed weekly dengue disease counts in
each of the 24 wards of Greater Mumbai city from 2010- 2015, to learn directed
and undirected models for the spread of dengue disease.

Our contributions in this paper are as follows.

  1. 1.

We formulate an optimization problem of maximum likelihood estimation (MLE) of
observed multivariate count data with $\ell_{1}$-type regularization to
estimate the partial correlation and causality graphs.

  2. 2.

A Monte Carlo expectation and maximization (MCEM) algorithm with
$\ell_{1}$-type regularization on multivariate count data is proposed to solve
the formulated problem.

  3. 3.

The proposed MCEM algorithm with $\ell_{1}$-type regularization is tested
successfully using simulated data.

  4. 4.

The partial correlation and causality graphs are estimated for the observed
weekly dengue disease counts from each ward of Greater Mumbai city.

##  2 Preliminaries and problem formulation

###  2.1 Parameter driven model

Let {$\mathbf{Y}(t)\in\mathbb{R}^{n},t=1\dots N\\}$ be an observed
multivariate time series of counts. An unobserved multivariate second order
stationary latent process {$\mathbf{X}(t)\in\mathbb{R}^{n},t=1\dots N\\}$ is
used to introduce the correlation between observations measured at successive
time points and individual time series. To model $\mathbf{Y}(t)$ using
$\mathbf{X}(t)$, a parameter driven model introduced in (Utazi et al., 2018)
is considered. Let $\textbf{z}_{t,i}$ be a $q\times 1$ vector of covariates at
time $t$ for the $i$th time series and $\bm{\beta}_{i}$ be a $q\times 1$
vector of regression coefficients corresponding to the $i$th covariate vector
$\textbf{z}_{t,i}$. Covariate vector $\textbf{z}_{t,i}$ is a function of time
$t$ and is usually used to model the explanatory variables such as trend,
seasonality etc. In this model, the counts $Y_{i}(t)$, given the latent
process $X_{i}(t)$ and the covariates $\textbf{z}_{t,i}$, follow Poisson
distribution with conditional mean
$\mu_{i}(t)=\text{exp}(\textbf{z}_{t,i}^{T}\bm{\beta}_{i}+X_{i}(t))$, which is
denoted by

$\displaystyle
Y_{i}(t)|X_{i}(t),\textbf{z}_{t,i}\sim\text{Poisson}(\mu_{i}(t)).$ (1)

We assume the latent process $X_{i}(t)$ in the model (1) follows an AR(p)
process given by

$\displaystyle
X_{i}(t)=\sum_{j=1}^{n}a_{ij}(1)X_{j}(t-1)+\dots+\sum_{j=1}^{n}a_{ij}(p)X_{j}(t-p)+\epsilon_{i}(t),$
(2)

where $a_{ij}(.)$’s are autoregressive coefficients and $\epsilon_{i}(t)$ is a
Gaussian white noise having zero mean and variance $\sigma_{i}^{2}$. Let
$A_{k}$ = {$a_{ij}(k)$} be an $n\times n$ matrix of autoregressive
coefficients and $\mathbf{\epsilon}(t)$ :=
($\epsilon_{1}(t),\dots,\epsilon_{n}(t)$) be iid Gaussian random variables
with covariance matrix $\Sigma$ :=
diag($\sigma_{1}^{2},\dots,\sigma_{n}^{2}$). For $k=1,\dots,p$, the vector
notation of the AR(p) multivariate latent process is given by

$\displaystyle\mathbf{X}(t)=\sum_{k=1}^{p}A_{k}\mathbf{X}(t-k)+\mathbf{\epsilon}(t),\qquad\mathbf{\epsilon}(t)\sim
N(\mathbf{0},\Sigma).$ (3)

Let
$\bm{\beta}\coloneqq\begin{bmatrix}\bm{\beta}_{1}&\ldots&\bm{\beta}_{n}\end{bmatrix}\in\mathbb{R}^{q\times
n}$ and
$\bm{\sigma}\coloneqq\begin{bmatrix}\sigma_{1}&\ldots&\sigma_{n}\end{bmatrix}^{T}\in\mathbb{R}^{n\times
1}$. Also assume that the parameter space is denoted by $\Omega$. Define the
operator vec($\cdot$), which converts the matrix into a vector by setting
columns of the matrix as a vector. Then,
$\bm{\theta}\coloneqq\text{vec}\Big{(}\begin{bmatrix}\bm{\beta}^{T}&A_{1}&\ldots&A_{p}&\bm{\sigma}\end{bmatrix}\Big{)}\in\Omega\subseteq\mathbb{R}^{n(q+pn+1)\times
1}$ represents the complete set of unknown parameters described by models (1)
and (3).

###  2.2 Partial correlation graphs

Let $G=(\mathscr{V},\mathscr{E},\overrightarrow{\mathscr{E}})$ be a mixed
graph (graph with directed and undirected edges), where $\mathscr{V}$ is a set
of vertices, $\mathscr{E}$ is a set of undirected edges
$\mathscr{E}\subseteq\mathscr{V}\times\mathscr{V}$ and
$\overrightarrow{\mathscr{E}}$ is a set of directed edges
$\overrightarrow{\mathscr{E}}\subseteq\mathscr{V}\times\mathscr{V}$. The
undirected edge between vertices $i$ and $j$ is denoted as $(i,j)$. If
$(i,j)\in\mathscr{E}$ then also $(j,i)\in\mathscr{E}$. The directed edge
between vertices $i$ and $j$ with direction from $i$ to $j$ is denoted as
$(i,j)$. If $(i,j)\in\overrightarrow{\mathscr{E}}$ then
$(j,i)\notin\overrightarrow{\mathscr{E}}$.

Let $\mathbf{X}_{\mathscr{V}\backslash\\{i,j\\}}\in\mathbb{R}^{(n-2)}$ denote
$\mathbf{X}(t)$ except the $i$th and $j$th elements i.e., $X_{i}(t)$ and
$X_{j}(t)$ and
$\mathbf{Y}_{\mathscr{V}\backslash\\{i,j\\}}\in\mathbb{R}^{(n-2)}$ denote
$\mathbf{Y}(t)$ except $Y_{i}(t)$ and $Y_{j}(t)$. The partial correlation
between any two time series of multivariate second order stationary processes
is given in (Dahlhaus, 2000). Define

$\mathcal{E}_{i}(t)=X_{i}(t)-\mu_{i}^{*}-\sum_{u}\mathbf{a}^{*}_{i}(t-u)\mathbf{X}_{\mathscr{V}\backslash\\{i,j\\}},$
(4)

as the error process between $X_{i}(t)$ and best linear filter based on
$\mathbf{X}_{\mathscr{V}\backslash\\{i,j\\}}$ that minimizes
$E\big{(}\mathcal{E}_{i}(t)\big{)}^{2}$. Here $\mu_{i}^{*}$ and
$\\{\mathbf{a}^{*}_{i}(u)\\}$ are the optimal linear filters, exact
expressions for which are available in (Brillinger, 1981). Similarly, the
error process for $X_{j}(t)$ is $\mathcal{E}_{j}(t)$. If the cross-covariance
between $\mathcal{E}_{i}(t)$ and $\mathcal{E}_{j}(t)$ for all time lag
$h\in\mathbb{Z}$ is zero, i.e.,
$Cov(\mathcal{E}_{i}(t),\mathcal{E}_{j}(t-h))=0$, $\forall h\in\mathbb{Z}$,
then $X_{i}(t)$ and $X_{j}(t)$ are defined to be partially uncorrelated given
$\mathbf{X}_{\mathscr{V}\backslash\\{i,j\\}}$.

From (W. J. Granger, 1969), the process $X_{i}(t)$ is said to be causal for
the another process $X_{j}(t)$ if the prediction of $X_{j}(t)$ using all
available information except $X_{i}(t)$, can be improved by adding the
available information about $X_{i}(t)$.

We define the partial correlations and causalities between observed time
series of counts in terms of the partial correlations and causalities between
second order stationary latent processes as follows.

###### Definition 1

Two time series $Y_{i}(t)$ and $Y_{j}(t)$ are partially uncorrelated given
$\mathbf{Y}_{\mathscr{V}\backslash\\{i,j\\}}$ if $X_{i}(t)$ and $X_{j}(t)$ are
partially uncorrelated given $\mathbf{X}_{\mathscr{V}\backslash\\{i,j\\}}$ and
$Y_{i}(t)$ is not causal for $Y_{j}(t)$ if $X_{i}(t)$ is not causal for
$X_{j}(t)$.

Now, the partial correlation and causality graph for the multivariate time
series $\mathbf{Y}(t)$ is defined as,

######  Definition 2 (Partial correlation and causality graph of
$\mathbf{Y}(t)$)

The Partial correlation and causality graph of $\mathbf{Y}(t)$ modeled in (1)
and (3) is the mixed graph
$G=(\mathscr{V},\mathscr{E},\overrightarrow{\mathscr{E}})$ with vertex set
$\mathscr{V}=\\{Y_{1}(t),\dots,Y_{n}(t)\\}$, undirected edge set $\mathscr{E}$
and directed edge set $\overrightarrow{\mathscr{E}}$ such that there is no
undirected edge between the $i^{th}$ and the $j^{th}$ node if and only if the
corresponding latent variables $X_{i}(t)$ and $X_{j}(t)$ are partially
uncorrelated given $\mathbf{X}_{\mathscr{V}\backslash\\{i,j\\}}$, and there is
no directed edge from the $i^{th}$ to the $j^{th}$ node if and only if
$X_{i}(t)$ is not causal for $X_{j}(t)$.

Let $R_{XX}(h)=E[(\mathbf{X}(t))(\mathbf{X}(t-h))^{T}]$ be the autocovariance
function of $\mathbf{X}(t)$ with time lag $h\in\mathbb{Z}$ and $S_{XX}(h)$ be
the spectral density matrix of $\mathbf{X}(t)$ i.e.,

$\displaystyle
S_{XX}(\omega)=\sum_{h=-\infty}^{\infty}R_{XX}(h)e^{-jh\omega}.$ (5)

It is known that the partial correlation graph can be obtained from the
inverse spectral density matrix of multivariate second order stationary
process (Brillinger, 1981) which is given in the theorem below.

###### Theorem 3

(Brillinger, 1981) Consider the stationary time series
$\mathbf{X}(t)\in\mathbb{R}^{n}$ and assume that the spectral density matrix
$S_{XX}(\omega)$ of $\mathbf{X}(t)$ is invertible for all $\omega$. Then,
$X_{i}(t)$ and $X_{j}(t)$ are partially uncorrelated given
$\mathbf{X}_{\mathscr{V}\backslash\\{i,j\\}}$ if and only if
$((S_{XX}(\omega))^{-1})_{ij}=0,\forall\omega$.

Thus from (Songsiri and Vandenberghe, 2010) and Theorem 3, the partial
correlation relations for the stationary latent process $\mathbf{X}(t)$, which
follows an AR model (3), in terms of inverse spectral are given in the lemma
below.

###### Lemma 4

$(S_{XX}(\omega)^{-1})_{ij}=0\Leftrightarrow(W_{k})_{ij}=0~{}\text{and}~{}(W_{k})_{ji}=0~{}\text{for}~{}k=0,\dots,p,$

where

$W_{k}=\begin{cases}-\Sigma^{-1}+\sum_{l=1}^{p}A_{l}^{T}\Sigma^{-1}A_{l},&k=0\\\
-2\Sigma^{-1}A_{k}+2\sum_{l=1}^{p-k}A_{l}^{T}\Sigma^{-1}A_{l+k},&k=1,\dots,p\end{cases}$
(6)

with $A_{0}=I$.

The following result is evident from Definition 1, Theorem 3 and (4).

###### Proposition 5

$Y_{i}(t)$ and $Y_{j}(t)$ are partially uncorrelated given
$\mathbf{Y}_{\mathscr{V}\backslash\\{i,j\\}}$ if
$(\bm{W}_{k})_{ij}=0~{}~{}\text{and}~{}~{}(\bm{W}_{k})_{ji}=0,~{}\text{for}~{}k=0,\dots,p$.

The following theorem states that the causality graph of the second order
stationary process $X(t)$ can be found in terms of AR coefficients given in
(3) (TjØstheim (1981); Hsiao (1982)).

###### Theorem 6

The time series $X_{i}(t)$ does not cause $X_{j}(t)$ if and only if the
corresponding components $A_{k}(j,i)$ vanish for all $k$ i.e.,

$\displaystyle X_{i}(t)~{}\text{does not cause}~{}X_{j}(t)\Longleftrightarrow
A_{k}(j,i)=0,~{}\forall k\in\\{1,\dots,p\\}.$ (7)

The following definitions are useful for future development in later sections.
Let $IW_{i}$ be the total incoming edge weight of $i^{th}$ node from other
nodes in the causality graph, i.e.,

$\displaystyle IW_{i}=\sum_{k=1}^{p}\sum_{j=1,j\neq i}^{n}A_{k}(i,j)$ (8)

Similarly, let $OW_{i}$ be the total outgoing edge weight of $i^{th}$ node to
other nodes in the causality graph, i.e.,

$\displaystyle OW_{i}=\sum_{k=1}^{p}\sum_{j=1,j\neq i}^{n}A_{k}(j,i)$ (9)

###  2.3 Problem formulation

From Theorem 3 and Proposition 5, we need to estimate the inverse spectral
density matrix of multivariate latent process $\mathbf{X}(t)$ to estimate the
partial correlations between the observed multiple time series of counts
$\mathbf{Y}(t)$. Since the number of possible edges in the partial correlation
and causality graph can be large (up to $n^{2}$), overfitting can be a
potential issue in this estimation problem (Cawley and Talbot, 2010). LASSO
(Tibshirani, 1996) is one of the regularization methods to overcome this
difficulty by setting certain parameters to zero. This technique is used here
to introduce sparsity in the inverse spectral density matrix. In LASSO, the
maximum likelihood estimation problem is regularized using the $\ell_{1}$-norm
of the partial correlation constraints given in Proposition 5. The
$\ell_{1}$-norm of a vector
$\mathbf{x}=\begin{bmatrix}x_{1}&\dots&x_{n}\end{bmatrix}^{T}\in\mathbb{R}^{n}$
is defined as $\|\mathbf{x}\|_{1}=|x_{1}|+\dots+|x_{n}|$. Let
$\mathbf{Y}:=\begin{bmatrix}\mathbf{Y}(1)\dots\mathbf{Y}(N)\end{bmatrix}\in\mathbb{R}^{n\times
N}$ be the observed data matrix and
$\mathbf{X}:=\begin{bmatrix}\mathbf{X}(1)\dots\mathbf{X}(N)\end{bmatrix}\in\mathbb{R}^{n\times
N}$ be the latent process data matrix. From (Utazi et al., 2018), the joint
probability density function (PDF) of $\mathbf{Y}$ and $\mathbf{X}$ is

$\displaystyle
P(\mathbf{Y},\mathbf{X},\bm{\theta})=P(\mathbf{Y},\bm{\theta}|\mathbf{X})P(\mathbf{X},\bm{\theta})$
$\displaystyle=\Big{\\{}\prod_{i=1}^{n}\prod_{t=1}^{N}P\big{(}Y_{i}(t),\bm{\theta}|X_{i}(t)\big{)}\prod_{t=p+1}^{N}P\big{(}\mathbf{X}(t),\bm{\theta}|\mathbf{X}(t-1),$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\dots,\mathbf{X}(t-p)\big{)}\Big{\\}}P\big{(}\mathbf{X}(1),\dots,\mathbf{X}(p),\bm{\theta}\big{)}$
$\displaystyle=\prod_{i=1}^{n}\prod_{t=1}^{N}\text{exp}\big{(}X_{i}(t)Y_{i}(t)+\mathbf{z}_{t,i}^{T}\bm{\beta}_{i}Y_{i}(t)-e^{X_{i}(t)+\mathbf{z}_{t,i}^{T}\bm{\beta}_{i}}\big{)}$
$\displaystyle\qquad\qquad\qquad\qquad\prod_{t=p+1}^{N}\frac{1}{\sqrt{2\pi|\Sigma|}}\text{exp}\Big{\\{}-\frac{1}{2}\big{(}\epsilon(t)^{T}\Sigma^{-1}\epsilon(t)\big{)}\Big{\\}}$
$\displaystyle\qquad\qquad\qquad\qquad
P\big{(}\mathbf{X}(1),\dots,\mathbf{X}(p),\bm{\theta}\big{)}.$ (10)

where $\epsilon(t)=\mathbf{X}(t)-\sum_{k=1}^{p}A_{k}\mathbf{X}(t-k)$ and
$P\big{(}\mathbf{X}(1),\dots,\mathbf{X}(p),\bm{\theta}\big{)}$ is the joint
probability density function of initial values
$\mathbf{X}(1),\dots,\mathbf{X}(p)$. The calculation of probability density
function $P\big{(}\mathbf{X}(1),\dots,\mathbf{X}(p),\bm{\theta}\big{)}$ is
given in (Helmut, 2005), p.29. Let

$\displaystyle\mathscr{X}(t)$ $\displaystyle:=\begin{bmatrix}\mathbf{X}(t)\\\
\mathbf{X}(t-1)\\\ \vdots\\\
\mathbf{X}(t-p+1)\end{bmatrix}\in\mathbb{R}^{np\times
1},\quad\mathscr{W}(t):=\begin{bmatrix}\epsilon(t)\\\ 0\\\ \vdots\\\
0\end{bmatrix}\in\mathbb{R}^{np\times 1},$ $\displaystyle\mathbf{A}$
$\displaystyle:=\begin{bmatrix}A_{1}&A_{2}&\dots&A_{p-1}&A_{p}\\\
I_{n}&0&\dots&0&0\\\ 0&I_{n}&\dots&0&0\\\
\vdots&\vdots&\cdots&\vdots&\vdots\\\
0&0&\dots&I_{n}&0\end{bmatrix}\in\mathbb{R}^{np\times np}.$ (11)

Then, the model (3) becomes

$\displaystyle\mathscr{X}(t)=\mathbf{A}\mathscr{X}(t-1)+\mathscr{W}(t).$ (12)

Note that $\mathbf{X}(t)$ is a second-order stationary process, and
$\mathbf{X}(t)$ and $\epsilon(t)$ are uncorrelated. Thus, from (12) the auto-
covariance function of $\mathscr{X}(t)$ with zero time lag is

$\displaystyle
E[\mathscr{X}(t)\mathscr{X}(t)^{T}]=\mathbf{A}E[\mathscr{X}(t-1)\mathscr{X}(t-1)^{T}]\mathbf{A}^{T}+E[\mathscr{W}(t)\mathscr{W}(t)^{T}],$
$\displaystyle
R_{\mathscr{X}\mathscr{X}}(0)=\mathbf{A}R_{\mathscr{X}\mathscr{X}}(0)\mathbf{A}^{T}+R_{\mathscr{W}\mathscr{W}}(0),$
(13)

where

$\displaystyle R_{\mathscr{X}\mathscr{X}}(0)$
$\displaystyle=\begin{bmatrix}R_{XX}(0)~{}~{}~{}R_{XX}(1)~{}\dots~{}R_{XX}(p-1)\\\
R_{XX}(1)^{T}~{}~{}R_{XX}(0)~{}\dots~{}R_{XX}(p-2)\\\
\vdots~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\vdots~{}~{}~{}~{}~{}~{}\cdots~{}~{}~{}~{}~{}~{}\vdots\\\
R_{XX}(p-1)^{T}~{}R_{XX}(p-2)^{T}~{}\dots~{}R_{XX}(0)\end{bmatrix}\in\mathbb{R}^{np\times
np},$ $\displaystyle R_{\mathscr{W}\mathscr{W}}(0)$
$\displaystyle=\begin{bmatrix}\Sigma&0&\dots&0\\\ 0&0&\dots&0\\\
\vdots&\vdots&\cdots&\vdots\\\ 0&0&\dots&0\end{bmatrix}\in\mathbb{R}^{np\times
np}.$ (14)

The vectorized equation of (13) is

$\displaystyle\text{vec}(R_{\mathscr{X}\mathscr{X}}(0))$
$\displaystyle=(\mathbf{A}\otimes\mathbf{A})\text{vec}(R_{\mathscr{X}\mathscr{X}}(0))+\text{vec}(R_{\mathscr{W}\mathscr{W}}(0)),$
$\displaystyle\text{vec}(R_{\mathscr{X}\mathscr{X}}(0))$
$\displaystyle=(I_{np}-\mathbf{A}\otimes\mathbf{A})^{-1}\text{vec}(R_{\mathscr{W}\mathscr{W}}(0)),$
(15)

where ‘$\otimes$’ denotes the Kronecker product. Therefore, the initial values
$\mathbf{X}(1),\dots,\mathbf{X}(p)$ follows a Gaussian distribution (Helmut,
2005) with zero mean and variance is $R_{\mathscr{X}\mathscr{X}}(0)$ which can
be calculated from (2.3).

The exact log-likelihood function of $\mathbf{Y}$ and $\mathbf{X}$ is
$l(\mathbf{Y},\mathbf{X},\bm{\theta}):=log(P(\mathbf{Y},\mathbf{X},\bm{\theta}))$.
Let

$\displaystyle
L(\mathbf{Y},\mathbf{X},\bm{\theta}):=l(\mathbf{Y},\mathbf{X},\bm{\theta})-\gamma
h_{1}(W_{0},W_{1},\ldots,W_{p})$ (16)

where $\gamma\geq 0$ is a regularization parameter and the regularization
function
$h_{1}(W_{0},W_{1},\ldots,W_{p})=\sum_{j>i}\sum_{k=0}^{p}\big{\\{}|(W_{k})_{ij}|+|(W_{k})_{ji}|\big{\\}}$
is the $\ell_{1}$-norm of off-diagonal elements of matrices
$W_{0},W_{1},\ldots,W_{p}$ given in (6). The amount of sparsity of the
estimated inverse spectral density matrix is controlled by the regularization
parameter $\gamma$. As $\gamma$ varies, the sparsity pattern varies in the
estimated inverse spectral density matrix from dense ($\gamma$ small) to
diagonal ($\gamma$ large). We aim to maximize regularized log-likelihood
$L(\mathbf{Y},\mathbf{X},\bm{\theta})$ with respect to $\bm{\theta}$ over
$\Omega$ to estimate the optimal parameters. To maximize
$L(\mathbf{Y},\mathbf{X},\bm{\theta})$ directly, data matrix $\mathbf{X}$
should be known. However, in reality the latent data matrix $\mathbf{X}$ is
unknown. The marginal PDF of $\mathbf{Y}$ is
$P(\mathbf{Y},\bm{\theta}):=\int_{\mathbf{X}}P(\mathbf{Y},\mathbf{X},\bm{\theta})d\mathbf{X}$.
Let $l(\mathbf{Y},\bm{\theta}):=log(P(\mathbf{Y},\bm{\theta}))$ be the
marginal log-likelihood of $\mathbf{Y}$. The marginal log-likelihood of
$\mathbf{Y}$ with $\ell_{1}$-type regularization

$\displaystyle L(\mathbf{Y},\bm{\theta}):=l(\mathbf{Y},\bm{\theta})-\gamma
h_{1}(W_{0},W_{1},\ldots,W_{p}),$ (17)

is used to estimate $\bm{\theta}$. We next formulate the optimization problem
to estimate $\bm{\theta}$.

###### Problem 2.1

Find

$\displaystyle\bm{\theta}^{*}=\text{arg}\hskip
0.85358pt\max_{\bm{\theta}\in\Omega}~{}L(\mathbf{Y},\bm{\theta}).$ (18)

###  2.4 Challenges in solution of Problem 2.1

To solve Problem 2.1 directly, we need to calculate
$L(\mathbf{Y},\bm{\theta})$ explicitly, which can be difficult because of $nN$
multiple integrals. To overcome this difficulty and unknown latent data matrix
$\mathbf{X}$, we use the expectation and maximization (EM) algorithm (Dempster
et al., 1977) to solve Problem 2.1. The EM algorithm maximizes the function
$L(\mathbf{Y},\bm{\theta})$ by working with
$L(\mathbf{Y},\mathbf{X},\bm{\theta})$ and the conditional density function
$P(\mathbf{X},\bm{\theta}|\mathbf{Y})$ of the latent data matrix $\mathbf{X}$
given observed count data matrix $\mathbf{Y}$.

Note that the standard EM algorithm is modified with the regularization term
in the following Algorithm 2.1. The tolerance ($\delta$) is adjusted according
to the accuracy requirement in the estimated parameters.

Algorithm 2.1 Expectation and Maximization (EM) algorithm with $\ell_{1}$-type
regularization

Input: For $k=1$, initial condition $\bm{\theta}^{(1)}$.

Output:  The estimated parameter is $\bm{\theta}^{*}$.

1:while $tol>\delta$ do

2: E-step:  $\displaystyle Q($
$\displaystyle\bm{\theta};\bm{\theta}^{(k)})=E_{\bm{\theta}^{(k)}}\\{l(\mathbf{Y},\mathbf{X},\bm{\theta})|\mathbf{Y}\\}-\gamma
h_{1}(W_{0},W_{1},\ldots,W_{p}),$
$\displaystyle=\int_{X}l(\mathbf{Y},\mathbf{X},\bm{\theta})P(\mathbf{X},\bm{\theta}^{(k)}|\mathbf{Y})d\mathbf{X}-\gamma
h_{1}(W_{0},W_{1},\ldots,W_{p}).$

3: M-step: Maximize the function $Q(\bm{\theta};\bm{\theta}^{(k)})$ with
respect to $\bm{\theta}$ over space $\Omega$ and update the maximizer,
$\displaystyle\bm{\theta}^{(k+1)}=\text{arg}\hskip
0.85358pt\max_{\bm{\theta}\in\Omega}Q(\bm{\theta};\bm{\theta}^{(k)}).$ (19)

4:
$tol=\frac{\|\bm{\theta}^{(k+1)}-\bm{\theta}^{(k)}\|}{\|\bm{\theta}^{(k)}\|}$.

5: $k=k+1$.

6:end while

$\bm{\theta}^{*}=\bm{\theta}^{(k)}$.

The log-likelihood function $l(\mathbf{Y},\mathbf{X},\bm{\theta})$ is required
to be known to solve the Algorithm 2.1. However, the conditional density
function $P(\mathbf{X},\bm{\theta}^{(k)}|\mathbf{Y})$ of the latent data
matrix $\mathbf{X}$ given observed count data matrix $\mathbf{Y}$ is
impossible to calculate at the $k$th step in the Algorithm 2.1 because it is a
mixture of Gaussian and Poisson distributions. Here Algorithm 2.1 cannot be
implemented directly. To overcome this difficulty, Monte Carlo techniques are
used to draw samples from $P(\mathbf{X},\bm{\theta}^{(k)}|\mathbf{Y})$ and
thereby numerically compute the integral given in the E-step of Algorithm 2.1.
The Metropolis-Hastings algorithm (Casella and Berger, 2002) is used to
generate the samples from $P(\mathbf{X},\bm{\theta}^{(k)}|\mathbf{Y})$.

The generation of samples from $P(\mathbf{X},\bm{\theta}^{(k)}|\mathbf{Y})$
using Metropolis-Hastings algorithm for the parameter driven model with AR(1)
process is given in (Utazi et al., 2018). In (Utazi et al., 2018), the
proposal distributions are calculated in Metropolis-Hastings algorithm with
AR(1) process. Similarly, in this paper the proposal distributions for
Metropolis-Hastings algorithm with AR(p) process are given. Let

$\displaystyle\Lambda$ $\displaystyle:=R_{\mathscr{X}\mathscr{X}}(0)^{-1}$
$\displaystyle:=\begin{bmatrix}\Lambda_{11}&\Lambda_{12}&\dots&\Lambda_{1p}\\\
\Lambda_{12}^{T}&\Lambda_{22}&\dots&\Lambda_{2p}\\\
\vdots&\vdots&\cdots&\vdots\\\
\Lambda_{1p}^{T}&\Lambda_{1(p-1)}^{T}&\dots&\Lambda_{pp}\end{bmatrix}.$ (20)

For $t=1,\dots,p$, the proposal distribution of $\mathbf{X}(t)$ given
$\mathbf{X}(1),\dots,\mathbf{X}(t-1)$ from the past, and
$\mathbf{X}(t+1),\dots,\mathbf{X}(t+p)$ to the future follows a Gaussian
distribution with variance,

$\displaystyle\Sigma_{t}=\Big{(}\Lambda_{tt}^{-1}+\sum_{i=0}^{t-1}A_{p-i}^{T}\Sigma^{-1}A_{p-i}\Big{)}^{-1},$
(21)

and mean,

$\displaystyle\mu_{t}$
$\displaystyle=\Sigma_{t}\Big{(}\sum_{i=1}^{t}A_{p+1-i}^{T}\Sigma^{-1}\mathbf{X}(p+t+1-i)-\sum_{\begin{subarray}{c}i=1\\\
i\neq t\end{subarray}}^{p}\Lambda_{ti}\mathbf{X}(i)$
$\displaystyle\quad-\sum_{j=1}^{t}\sum_{\begin{subarray}{c}i=1\\\ i\neq
p+1-j\end{subarray}}^{p}A_{p+1-j}^{T}\Sigma^{-1}A_{i}\mathbf{X}(p+t+1-j-i)\Big{)}.$
(22)

For $t=p+1,\dots,N-p$, the proposal distribution of $\mathbf{X}(t)$ given
$\mathbf{X}(t-1),\dots,\mathbf{X}(t-p)$ from the past, and
$\mathbf{X}(t+1),\dots,\mathbf{X}(t+p)$ to the future follows a Gaussian
distribution with variance,

$\displaystyle\Sigma_{t}=\Big{(}\Sigma^{-1}+\sum_{i=1}^{p}A_{i}^{T}\Sigma^{-1}A_{i}\Big{)}^{-1},$
(23)

and mean,

$\displaystyle\mu_{t}$
$\displaystyle=\Sigma_{t}\Big{(}\sum_{i=1}^{p}\Sigma^{-1}A_{i}\mathbf{X}(t-i)+\sum_{i=1}^{p}A_{i}^{T}\Sigma^{-1}\mathbf{X}(t+i)$
$\displaystyle\quad-\sum_{j=1}^{p}\sum_{\begin{subarray}{c}i=1\\\ i\neq
j\end{subarray}}^{p}A_{j}^{T}\Sigma^{-1}A_{i}\mathbf{X}(t+j-i)\Big{)}.$ (24)

For $t=N-p+1,\dots,N$, the proposal distribution of $\mathbf{X}(t)$ given
$\mathbf{X}(t-1),\dots,\mathbf{X}(t-p)$ from the past, and
$\mathbf{X}(t+1),\dots,\mathbf{X}(N)$ to the future follows a Gaussian
distribution with variance,

$\displaystyle\Sigma_{t}=\Big{(}\Sigma^{-1}+\sum_{i=1}^{N-t}A_{i}^{T}\Sigma^{-1}A_{i}\Big{)}^{-1},$
(25)

and mean,

$\displaystyle\mu_{t}$
$\displaystyle=\Sigma_{t}\Big{(}\sum_{i=1}^{p}\Sigma^{-1}A_{i}\mathbf{X}(t-i)+\sum_{i=1}^{N-t}A_{i}^{T}\Sigma^{-1}\mathbf{X}(t+i)$
$\displaystyle\quad-\sum_{j=1}^{N-t}\sum_{\begin{subarray}{c}i=1\\\ i\neq
j\end{subarray}}^{p}A_{j}^{T}\Sigma^{-1}A_{i}\mathbf{X}(t+j-i)\Big{)}.$ (26)

Using these proposal distributions in Metropolis-Hastings algorithm given in
Algorithm 2.2, we generate $m$ samples $\mathbf{X}^{(1)}$,
$\mathbf{X}^{(2)}$,…, $\mathbf{X}^{(m)}\in\mathbb{R}^{n\times N}$ from
conditional density function $P(\mathbf{X},\bm{\theta}^{(k)}|\mathbf{Y})$.
Thus, we propose a Monte Carlo expectation and maximization (MCEM) algorithm
with $\ell_{1}$-type regularization next.

Algorithm 2.2 Metropolis-Hastings algorithm

Input: Initial values
$\mathbf{X}^{(1)}:=[\mathbf{X}^{(1)}(1)~{}\mathbf{X}^{(1)}(2)\cdots\mathbf{X}^{(1)}(N)]\in\mathbb{R}^{n\times
N}$ and number of samples $m$.

Output:  Generated samples
$\mathbf{X}^{(1)},\mathbf{X}^{(2)},\dots,\mathbf{X}^{(m)}\in\mathbb{R}^{n\times
N}$.

1:for $r=1,\dots,(m-1)$ do

2: for $t=1,\dots,N$ do

3: for $i=1,\dots,n$ do

4: Generate uniform random number $U_{i}$ and generate $X_{i}(t)$ from the
proposal distributions (21)-(2.4).

5: Calculate the acceptance probability,
$\displaystyle\rho_{i}=\text{min}\bigg{\\{}\frac{\text{exp}(X_{i}(t))Y_{i}(t)-\text{exp}(X_{i}(t)+\mathbf{z}_{t,i}^{T}\bm{\beta}_{i})}{\text{exp}(X_{i}^{(r)}(t))Y_{i}(t)-\text{exp}(X_{i}^{(r)}(t)+\mathbf{z}_{t,i}^{T}\bm{\beta}_{i})},1\bigg{\\}},$
where $X_{i}^{(r)}(t)$ is the $i$th element of $\mathbf{X}^{(r)}(t)$.

6: Update $X^{(r+1)}_{i}(t)$:  $\displaystyle
X^{(r+1)}_{i}(t)=\begin{cases}X_{i}(t),&\text{if}~{}U_{i}\leq\rho_{i}\\\
X_{i}^{(r)}(t),&\text{otherwise}\end{cases}.$

7: end for

8: end for

9:end for

##  3 Main results

In this section we discuss the Monte Carlo expectation and maximization (MCEM)
algorithm with $\ell_{1}$-type regularization. Further, the asymptotic
convergence of sequence of parameters generated by the MCEM algorithm with
$\ell_{1}$-type regularization is proved.

###  3.1 Monte Carlo Expectation and Maximization (MCEM) algorithm with
$\ell_{1}$-type regularization

In this section, we propose a MCEM algorithm with $\ell_{1}$-type
regularization to estimate parameters $\bm{\theta}$. As mentioned above, we
approximate the integral in the $E$-$step$ (step 2) of Algorithm 2.1. The MCEM
algorithm with $\ell_{1}$-type regularization is presented in Algorithm 3.1.
The tolerance ($\delta$) is adjusted according to the accuracy requirement in
the estimated parameters.

Algorithm 3.1 Monte Carlo Expectation and Maximization (MCEM) algorithm with
$\ell_{1}$-type regularization

Input: For $k=1$, initial condition $\bm{\theta}^{(1)}_{m}$.

Output:  The estimated parameter is $\bm{\theta}^{*}$.

1:while $tol>\delta$ do

2: Generate $m$ samples
$\mathbf{X}^{(1)},\mathbf{X}^{(2)},\dots,\mathbf{X}^{(m)}\in\mathbb{R}^{n\times
N}$ from the conditional PDF $P(\mathbf{X},\bm{\theta}^{(k)}_{m}|\mathbf{Y})$
using Metropolis-hastings algorithm 2.2.

3: E-step:
$Q_{m}(\bm{\theta};\bm{\theta}^{(k)}_{m})=\frac{1}{m}\sum_{i=1}^{m}l(\mathbf{Y},\mathbf{X}^{(i)},\bm{\theta})-\gamma
h_{1}(W_{0},W_{1},\ldots,W_{p}),$

4: M-step: Maximize the function (3) with respect to $\bm{\theta}$ over space
$\Omega$ and update the maximizer,
$\displaystyle\bm{\theta}^{(k+1)}_{m}=\text{arg}\hskip
0.85358pt\max_{\bm{\theta}\in\Omega}~{}Q_{m}(\bm{\theta};\bm{\theta}^{(k)}_{m}).$
(27)

5:
$tol=\frac{\|\bm{\theta}^{(k+1)}_{m}-\bm{\theta}^{(k)}_{m}\|}{\|\bm{\theta}^{(k)}_{m}\|}$.

6: $k=k+1$.

7:end while

$\bm{\theta}^{*}=\bm{\theta}^{(k)}_{m}$.

From this algorithm, the generated sequence is
{$\bm{\theta}^{(k)}_{m}:k=1,2,\dots$} and the estimated parameter is
$\bm{\theta}^{*}$. This estimated parameter $\bm{\theta}^{*}$ is the local
maximizer of the function $L(\mathbf{Y},\bm{\theta})$. We prove this claim in
the next section. Note that the generated sequence
{$\bm{\theta}^{(k)}_{m}:k=1,2,\dots$} is a sequence of random variables.

###  3.2 Asymptotic results

We obtain the sequence {$\bm{\theta}^{(k)}_{m}:k=1,2,\dots$} from the Monte
Carlo expectation and maximization (MCEM) algorithm with $\ell_{1}$-type
regularization given in Algorithm 3.1. We need the following lemmas to prove
the asymptotic convergence of this sequence. From step-3 of Algorithm 2.1,
define $M(\bm{\theta}^{(k)}):=\text{arg}\hskip
0.85358pt\max_{\bm{\theta}\in\Omega}Q(\bm{\theta};\bm{\theta}^{(k)})$ and
assume that $M$ is a continuous function on $\Omega$.

###### Lemma 7

$\bm{\theta}^{(k+1)}_{m}$ converges in probability to
$M(\bm{\theta}^{(k)}_{m})$ as $m\rightarrow\infty$.

Proof  The proof of this lemma follows from the convergence results of
Metropolis-Hastings algorithm (Mengersen and Tweedie, 1996; Roberts and Smith,
1994; Roberts and Tweedie, 1996).

###### Lemma 8

Let {$\bm{\theta}^{(k)}:k=1,2,\dots$} be a sequence generated by Algorithm
2.1. Then, the sequence $\\{L(\mathbf{Y},\bm{\theta}^{(k)})\\}$ is a non-
decreasing sequence.

Proof

The proof follows from page 78 of (McLachlan and Krishnan, 2008).

###### Lemma 9

Let {$\bm{\theta}^{(k)}_{m}:k=1,2,\dots$} be a sequence generated by Algorithm
3.1. Then,

$\displaystyle P(\\{L(\mathbf{Y},\bm{\theta}^{(k)}_{m})\leq
L(\mathbf{Y},\bm{\theta}^{(k+1)}_{m})\\})\rightarrow
1~{}~{}as~{}~{}m\rightarrow\infty.$ (28)

Proof  Consider the function given in E-step (step-3) of Algorithm 3.1,

$\displaystyle
Q_{m}(\bm{\theta};\bm{\theta}^{(k)}_{m})=\frac{1}{m}\sum_{i=1}^{m}l(\mathbf{Y},\mathbf{X}^{(i)},\bm{\theta})-\gamma
h_{1}(W_{0},W_{1},\ldots,W_{p}).$ (29)

The log-likelihood function
$l(\mathbf{Y},\mathbf{X}^{(i)},\bm{\theta})=l(\mathbf{Y},\bm{\theta})+l(\mathbf{X}^{(i)}|\mathbf{Y},\bm{\theta})$
where $l(\mathbf{Y},\bm{\theta}):=log(P(\mathbf{Y},\bm{\theta}))$ and
$l(\mathbf{X}^{(i)}|\mathbf{Y},\bm{\theta}):=log(P(\mathbf{X}^{(i)}|\mathbf{Y},\bm{\theta}))$.
Replace this joint log-likelihood in (29),

$\displaystyle Q_{m}(\bm{\theta};\bm{\theta}^{(k)}_{m})=\frac{1}{m}$
$\displaystyle\sum_{i=1}^{m}l(\mathbf{X}^{(i)}|\mathbf{Y},\bm{\theta})$
$\displaystyle+\big{(}l(\mathbf{Y},\bm{\theta})-\gamma
h_{1}(W_{0},W_{1},\ldots,W_{p})\big{)}.$ (30)

Let
$H_{m}(\bm{\theta};\bm{\theta}^{(k)}_{m}):=\frac{1}{m}\sum_{i=1}^{m}l(\mathbf{X}^{(i)}|\mathbf{Y},\bm{\theta})$.
Then, (3.2) becomes

$\displaystyle
Q_{m}(\bm{\theta};\bm{\theta}^{(k)}_{m})=H_{m}(\bm{\theta};\bm{\theta}^{(k)}_{m})+L(\mathbf{Y},\bm{\theta}).$
(31)

The function $L(\mathbf{Y},\bm{\theta})$ at
$\bm{\theta}=\bm{\theta}^{(k)}_{m}$ in (31) is

$\displaystyle
L(\mathbf{Y},\bm{\theta}^{(k)}_{m})=Q_{m}(\bm{\theta}^{(k)}_{m};\bm{\theta}^{(k)}_{m})-H_{m}(\bm{\theta}^{(k)}_{m};\bm{\theta}^{(k)}_{m}).$
(32)

The function $L(\mathbf{Y},\bm{\theta})$ at
$\bm{\theta}=\bm{\theta}^{(k+1)}_{m}$ in (31) is

$\displaystyle
L(\mathbf{Y},\bm{\theta}^{(k+1)}_{m})=Q_{m}(\bm{\theta}^{(k+1)}_{m};\bm{\theta}^{(k)}_{m})-H_{m}(\bm{\theta}^{(k+1)}_{m};\bm{\theta}^{(k)}_{m}).$
(33)

Subtract (32) from (33), then

$\displaystyle\\{L$
$\displaystyle(\mathbf{Y},\bm{\theta}^{(k+1)}_{m})-L(\mathbf{Y},\bm{\theta}^{(k)}_{m})\\}=\\{Q_{m}(\bm{\theta}^{(k+1)}_{m};\bm{\theta}^{(k)}_{m})$
$\displaystyle-
Q_{m}(\bm{\theta}^{(k)}_{m};\bm{\theta}^{(k)}_{m})\\}-\\{H_{m}(\bm{\theta}^{(k+1)}_{m};\bm{\theta}^{(k)}_{m})-H_{m}(\bm{\theta}^{(k)}_{m};\bm{\theta}^{(k)}_{m})\\}.$
(34)

From Algorithm 3.1, we know that

$\displaystyle Q_{m}(\bm{\theta}^{(k+1)}_{m};\bm{\theta}^{(k)}_{m})\geq
Q_{m}(\bm{\theta}^{(k)}_{m};\bm{\theta}^{(k)}_{m}).$ (35)

Consider the second part of right hand side of (3.2),

$\displaystyle H_{m}(\bm{\theta}^{(k+1)}_{m};\bm{\theta}^{(k)}_{m})$
$\displaystyle-H_{m}(\bm{\theta}^{(k)}_{m};\bm{\theta}^{(k)}_{m})$
$\displaystyle=\frac{1}{m}\sum_{i=1}^{m}\big{(}l(\mathbf{X}^{(i)}|\mathbf{Y},\bm{\theta}^{(k+1)}_{m})-l(\mathbf{X}^{(i)}|\mathbf{Y},\bm{\theta}^{(k)}_{m})\big{)}$
$\displaystyle=\frac{1}{m}\sum_{i=1}^{m}log\Big{\\{}\frac{P(\mathbf{X}^{(i)}|\mathbf{Y},\bm{\theta}^{(k+1)}_{m})}{P(\mathbf{X}^{(i)}|\mathbf{Y},\bm{\theta}^{(k)}_{m})}\Big{\\}}$
$\displaystyle\leq
log\Big{\\{}\frac{1}{m}\sum_{i=1}^{m}\frac{P(\mathbf{X}^{(i)}|\mathbf{Y},\bm{\theta}^{(k+1)}_{m})}{P(\mathbf{X}^{(i)}|\mathbf{Y},\bm{\theta}^{(k)}_{m})}\Big{\\}}.$
(36)

Consider right hand side of (3.2), from law of large numbers (Roberts and
Tweedie, 1996),

$\displaystyle\frac{1}{m}\sum_{i=1}^{m}\frac{P(\mathbf{X}^{(i)}|\mathbf{Y},\bm{\theta}^{(k+1)}_{m})}{P(\mathbf{X}^{(i)}|\mathbf{Y},\bm{\theta}^{(k)}_{m})}$
$\displaystyle\xrightarrow{p}\int_{X}\frac{P(\mathbf{X}|\mathbf{Y},\bm{\theta}^{(k+1)}_{m})}{P(\mathbf{X}|\mathbf{Y},\bm{\theta}^{(k)}_{m})}P(\mathbf{X}|\mathbf{Y},\bm{\theta}^{(k)}_{m})~{}d\mathbf{X}$
$\displaystyle=\int_{X}P(\mathbf{X}|\mathbf{Y},\bm{\theta}^{(k+1)}_{m})~{}d\mathbf{X}$
$\displaystyle=1.$ (37)

as $m\rightarrow\infty$. Then the inequality (3.2) holds in probability,

$\displaystyle
H_{m}(\bm{\theta}^{(k+1)}_{m};\bm{\theta}^{(k)}_{m})-H_{m}(\bm{\theta}^{(k)}_{m};\bm{\theta}^{(k)}_{m})\leq
0,$ (38)

as $m\rightarrow\infty$. Therefore, from (35) and (38), the equation (3.2)
becomes,

$\displaystyle L(\mathbf{Y},\bm{\theta}^{(k+1)}_{m})\geq
L(\mathbf{Y},\bm{\theta}^{(k)}_{m}).$ (39)

in probability as $m\rightarrow\infty$.

Following (Chan and Ledolter, 1995), Lemma 7, Lemma 8 and Lemma 9, the
asymptotic convergence results of sequence
{$\bm{\theta}^{(k)}_{m}:k=1,2,\dots$} follows from (Chan and Ledolter, 1995).

###### Theorem 10

Let {$\bm{\theta}^{(k)}_{m}:k=1,2,\dots$} be a sequence generated by Algorithm
3.1 based on sample size $m$. Suppose $\bm{\theta}^{*}$ is an isolated local
maximizer of function $L(\mathbf{Y},\bm{\theta})$. For any $\epsilon>0$, there
exist $K_{0}<\infty$ and $\delta>0$ such that for any starting value
$\bm{\theta}^{(1)}_{m}\in\mathscr{N}:=\\{\bm{\theta}:\|\bm{\theta}-\bm{\theta}^{*}\|\leq\delta\\}$,

$\displaystyle
P(\\{\|\bm{\theta}^{(k)}_{m}-\bm{\theta}^{*}\|<\epsilon~{}\text{for~{}some}~{}k\leq
K_{0}\\})\rightarrow 1~{}~{}as~{}~{}m\rightarrow\infty.$ (40)

##  4 Inference on Simulated data

The performance of the MCEM algorithm with $\ell_{1}$-type regularization to
estimate the partial correlation and causality graphs, is tested using
randomly generated data in this section.

We consider a parameter driven model (1) to generate counts randomly. The
increasing trend and yearly seasonality terms are considered in the true model
with number of time series $n=10$. The delay in the AR process (3) is taken to
be $p=2$. The conditional mean in the true parameter driven model (1) is

$\displaystyle\mu_{i}(t)=\text{exp}(\mathbf{z}_{t,i}^{T}\bm{\beta}_{i}+X_{i}(t))~{}\text{for}~{}i=1,2,\dots,10,$
(41)

where $\mathbf{z}_{t,i}=\begin{bmatrix}1&t&\text{cos}(2\pi
t/12)&\text{sin}(2\pi t/12)\end{bmatrix}$ and
$\bm{\beta}_{i}=\begin{bmatrix}0.5&0.005&0.5&0.5\end{bmatrix}$ for all
$i=1,2,\dots,10$. The AR(2) multivariate latent process
$\mathbf{X}(t)\in\mathbb{R}^{10}$ is

$\displaystyle\mathbf{X}(t)=\sum_{k=1}^{2}A_{k}\mathbf{X}(t-k)+\mathbf{\epsilon}(t),\qquad\mathbf{\epsilon}(t)\sim
N(\mathbf{0},\Sigma)$ (42)

where $10\times 10$ matrices $A_{1}$ and $A_{2}$ are randomly chosen with
elements $\pm 0.3$ and zeros. The covariance matrix of Gaussian noise,
$\Sigma$ is the diagonal matrix with diagonal elements equal to 0.01. We
generate the Poisson counts for $N=200$ from this true model. Let
$\mathbf{Y}:=\begin{bmatrix}\mathbf{Y}(1)&\ldots&\mathbf{Y}(200)\end{bmatrix}$
be an observed data matrix. From the observed multivariate data $\mathbf{Y}$,
the partial correlation and causality graph is estimated next. The partial
correlations and causalities between multiple count time series for this true
model (41) are shown in separate graphs given in Fig. 4(a) and Fig. 5(a).

###  4.1 Choice of regularization parameter $\gamma$

The sparsity in the inverse spectral density matrix is controlled by the
regularization parameter $\gamma$. As $\gamma$ varies, the sparsity in the
inverse spectral density matrix changes. There are several methods to estimate
the regularization parameter $\gamma$. Cross-validation (J.Z. Huang and Liu,
2006) is one of the methods to estimate $\gamma$. This method is not accurate
and it requires significant computations to get an accurate $\gamma$ if the
length of observed data is less. A method to select the ‘best’ regularization
parameter was given in (Songsiri, 2010) based on information scores. In this
method, the regularization parameter is estimated using the trade-off curve
(Songsiri, 2010) between the log-likelihood
($l(\mathbf{Y},\mathbf{X},\bm{\theta})$) and the $\ell_{1}$-type
regularization function ($h_{1}(W_{0},W_{1},\dots,W_{p})$). We collect several
$\gamma$ values from this trade-off curve. We use a further thresholding on
the elements of the inverse spectral density matrix for each value of $\gamma$
along the trade-off curve. This method is as follows.

Partial coherence spectrum : The inverse spectral density matrix normalized
with $diag(S_{XX}(\omega)^{-1})$ where $diag(S_{XX}(\omega)^{-1})$ is the
diagonal of $S_{XX}(\omega)^{-1}$, is called as the partial coherence spectrum
$R(\omega)$.

$R(\omega)=diag(S_{XX}(\omega)^{-1})^{-1/2}S_{XX}(\omega)^{-1}diag(S_{XX}(\omega)^{-1})^{-1/2}$
(43)

Thresholding: Let $\rho_{ij}$ be the $\ell_{\infty}$-norm of the entries for
the partial coherence spectrum $R(\omega)$ with respect to $\omega$ i.e.,

$\rho_{ij}=\sup_{\omega}|R(\omega)_{ij}|.$ (44)

The value of $\rho_{ij}$ signifies the partial correlation between $Y_{i}(t)$
and $Y_{j}(t)$ given $\mathbf{Y}_{\mathscr{V}\backslash\\{i,j\\}}$ and it
ranges from 0 to 1. The value of $\rho_{ij}$ being low signifies that the
partial correlation between $Y_{i}(t)$ and $Y_{j}(t)$ given
$\mathbf{Y}_{\mathscr{V}\backslash\\{i,j\\}}$ is negligible. A thresholding
approach to remove such negligible partial correlation has been discussed in
(Lounici, 2008). Let $\rho^{*}$ be a threshold value chosen by this method. If
$\rho_{ij}\leq\rho^{*}$, then remove the edge $(i,j)$ from the graph. Collect
all the partial coherence spectrum $(R(\omega))$ corresponding $\gamma$ values
along the trade-off curve by applying the threshold. Then assign the ranks
using Bayes information criteria (BIC) scores and select the partial
correlation graph which has the lowest score.

###  4.2 Topology selection from observed count data

For the observed multivariate count data $\mathbf{Y}(t)$ for $t=1,\dots,200$,
we applied Algorithm 3.1. We found $\gamma$ values 0, 0.0698, 0.2911, 0.5857,
0.6872, 0.9963, 1.8527, 2.6891 and 3 from the trade-off curve which is given
in Fig. 1. For varies values of regularization parameter on the trade off
curve, the partial coherence spectrum $R(\omega)$ from (43) is calculated by
setting entries with $\rho_{ij}\leq 0.1$ to zero. Thus, all the partial
coherence spectrum after thresholding corresponding to these $\gamma$ values
are given in Fig.2. We observed that sparsity has increased as we move from
$\gamma=0$ to $\gamma=3$. Next we rank all these partial coherence spectrum
with BIC score. The BIC scores corresponding to the $\gamma$ values are given
in Fig. 3. From the Fig. 3, it is apparent that the best regularization
parameter is selected to be $\gamma^{*}=0.2911$. The estimated partial
correlations and the causalities between multiple count time series
corresponding to $\gamma^{*}=0.2911$ are shown in separate graphs given in
Fig. 4(b) and Fig. 5(b). From the true and estimated partial correlation
graphs in Fig. 4, it is observed that the proposed algorithm has misclassified
only three edges as zeros and two edges as non-zeros. Similarly, from the true
and estimated causality graphs Fig .5, it is seen that only four directed
edges are incorrectly estimated.

Figure 1: Trade-off curve between the conditional log-likelihood
$l(\mathbf{Y},\mathbf{X},\bm{\theta})$ and $h_{1}(W_{0},W_{1},\dots,W_{p})$
Figure 2: The partial correlation graphs along the trade of curve in Fig. 1
Figure 3: The BIC scores along the trade of curve in Fig. 1

$Y_{1}(t)$$Y_{2}(t)$$Y_{3}(t)$$Y_{4}(t)$$Y_{5}(t)$$Y_{6}(t)$$Y_{7}(t)$$Y_{8}(t)$$Y_{9}(t)$$Y_{10}(t)$
(a) True partial correlation graph

$Y_{1}(t)$$Y_{2}(t)$$Y_{3}(t)$$Y_{4}(t)$$Y_{5}(t)$$Y_{6}(t)$$Y_{7}(t)$$Y_{8}(t)$$Y_{9}(t)$$Y_{10}(t)$
(b) Estimated partial correlation graph ($\gamma^{*}=0.2911$)

Figure 4: True and estimated partial correlation graphs

$Y_{1}(t)$$Y_{2}(t)$$Y_{3}(t)$$Y_{4}(t)$$Y_{5}(t)$$Y_{6}(t)$$Y_{7}(t)$$Y_{8}(t)$$Y_{9}(t)$$Y_{10}(t)$
(a) True causality graph

$Y_{1}(t)$$Y_{2}(t)$$Y_{3}(t)$$Y_{4}(t)$$Y_{5}(t)$$Y_{6}(t)$$Y_{7}(t)$$Y_{8}(t)$$Y_{9}(t)$$Y_{10}(t)$
(b) Estimated causality graph ($\gamma^{*}=0.2911$)

Figure 5: True and estimated causality graphs

##  5 Dengue spread in Greater Mumbai

Dengue is a viral infectious disease transmitted to humans through the bites
of infected Aedes aegypti mosquitoes. It is one of the most severe health
problems being faced globally. The Mumbai Metropolitan Area (previously
Greater Mumbai Metropolitan Area) is divided into twenty-four administrative
divisions known as wards. The dengue counts are observed weekly from January
2010 to December 2015 from the wards of Greater Mumbai by the public health
department of the Municipal Corporation of Greater Mumbai (MCGM). The length
of each time series is 318.

###  5.1 Model Estimation

In this section, we estimate the partial correlation and causality graph using
the proposed Algorithm 3.1 from the observed 24-dimensional dengue data set.
The time series plot of normalized dengue counts of some of the wards (A, B,
and C) are shown in Fig. 6.

###### Remark 11

In the plots below, the absolute counts are not shown due to data privacy
rules of MCGM. However the proposed algorithm and the final results are based
on the actual counts.

From these time series plots, the increasing trend and yearly seasonality with
high dengue counts being reported during monsoon months are observed. The non-
stationary behaviour of the dengue counts is also apparent from these plots.
Thus, a linear trend and yearly seasonality represented by a pair of sine and
cosine terms, are considered in the covariate vector $\textbf{z}_{t,i}$ for
all wards $i=1,\dots,24$ of the model (1). Then model (1) becomes

$\displaystyle
Y_{i}(t)|X_{i}(t),\textbf{z}_{t,i}\sim\text{Poisson}(\mu_{i}(t)).$ (45)

where $\textbf{z}_{t,i}=\begin{bmatrix}1&t&cos(2\pi t/52)&sin(2\pi
t/52)\end{bmatrix}^{T}$.

Figure 6: Time series plot of normalized dengue counts of A, B and C wards
where count* indicates the true count normalized with respect to the maximum
weekly count observed over the observed duration (318 weeks)

We estimate the unknown parameters of the model (45) and (3) with order
ranging from $p=0$ to $p=3$ using Algorithm 3.1. All these models are ranked
with second-order Akaike (AIC${}_{\text{c}}$) and Bayes information criteria
(BIC) which are given in Table. 1. The model of order $p=1$ which corresponds
to the lowest AIC${}_{\text{c}}$ and BIC score is selected.

| AIC${}_{\text{c}}$ | BIC
---|---|---
$p=0$ | 8924.5 | 13195.8
$p=1$ | 7176.2 | 12851.1
$p=2$ | 7377.6 | 16530.3
$p=3$ | 7405.1 | 19844.9
Table 1: AIC${}_{\text{c}}$ and BIC scores for different values of order $p$

Thus for the model with order $p=1$, we find the trade-off curve between the
conditional log-likelihood ($l(\mathbf{Y},\mathbf{X},\bm{\theta})$) and the
$\ell_{1}$-type regularization function ($h_{1}(W_{0},W_{1},\dots,W_{p})$)
which is shown in Fig. 7. After thresholding (see section 4.1) with threshold
$\rho^{*}=0.1$, the normalized inverse spectral density matrices from (43)
along the trade-off curve (Fig. 7) are given in Fig. 8. We observe that
sparsity pattern varies in the estimated inverse spectral density matrix from
dense ($\gamma$ small) to diagonal ($\gamma$ large). We rank all of the
partial correlation graphs along the trade-off curve with BIC scores which is
shown in Fig. 9. The regularization parameter ($\gamma$) values along the
trade-off curve and the corresponding BIC scores are given in Table. 2. The
estimated partial correlations and causalities between multiple count time
series corresponding to $\gamma^{*}=0.682$ are shown in separate graphs given
in Fig.10 and Fig. 11.

Figure 7: Trade-off curve between log-likelihood function ($l(\mathbf{Y},\mathbf{X},\bm{\theta})$) and regularization function ($h_{1}(W_{0},W_{1},\ldots,W_{p})$) Figure 8: Partial correlation graphs along the trade-off curve in Fig. 7 7 Figure 9: BIC scores along the trade-off curve in Fig. 7 $\gamma$ | BIC | $\gamma$ | BIC | $\gamma$ | BIC
---|---|---|---|---|---
0 | 12111.43 | 0.7969 | 10227.05 | 2.0186 | 10973.20
0.1706 | 10846.09 | 0.9795 | 10313.73 | 2.2213 | 11275.20
0.5084 | 10583.76 | 1.1266 | 10512.76 | 2.6271 | 11510.35
0.6829 | 10197.97 | 1.5951 | 10734.25 | 3.35 | 11811.79
Table 2: Regularization parameter ($\gamma$) values along the trade-off curve
in Fig. 7 and corresponding BIC scores
ABCDEF/NF/SG/NG/SH/EH/WK/EK/WP/NP/SR/CR/NR/SLM/EM/WNST Figure 10: Estimated
partial correlation graph with regularization ($\gamma^{*}=0.682$)
ABCDEF/NF/SG/NG/SH/EH/WK/EK/WP/NP/SR/CR/NR/SLM/EM/WNST Figure 11: Estimated
causality graph

###  5.2 Inference

Based on the number of undirected edges between the wards in the estimated
partial correlation graph in Fig. 10, we construct a colour coded map (see
Fig. 12) showing the number of incident edges for each ward. Similarly, based
on the number of outward edges of each ward in the estimated causality graph
of Fig. 11, we construct a colour coded map (see Fig. 13) showing the number
of outgoing edges for each ward. All estimated model parameters such as
undirected and directed edges for each ward along with normalized disease
count (see Remark 11), are presented in the Table 3 .

Figure 12: Each ward of Greater Mumbai is colour coded with the number of
undirected edges in the estimated partial correlation graph Fig. 10 Figure 13:
Each ward of Greater Mumbai is colour coded with the number of outgoing edges
in the estimated causality graph Fig. 11

The plots of the normalized dengue counts, number of undirected edges from the
estimated partial correlation graph Fig. 10 and number of incoming and
outgoing edges from the estimated causality graph Fig. 11 of each ward, are
shown in Fig. 14. In Fig .14, all counts are normalized with their
corresponding maximum.

Ward | Count | UE | IE | OE | $IW$ | $OW$
---|---|---|---|---|---|---
A | 0.252 | 0.6 | 0.411 | 0.823 | 0.508 | 0.631
B | 0.217 | 0.466 | 0.529 | 0.647 | 0.369 | 0.486
C | 0.201 | 0.066 | 0.176 | 0.470 | 0.202 | 0.295
D | 0.278 | 0.333 | 0.411 | 0.294 | 0.343 | 0.222
E | 0.984 | 0.933 | 0.764 | 0.705 | 0.799 | 0.429
F/N | 0.668 | 0.8 | 0.941 | 0.529 | 0.885 | 0.307
F/S | 0.481 | 0.6 | 0.647 | 0.588 | 0.624 | 0.430
G/N | 0.586 | 0.866 | 0.882 | 0.705 | 0.937 | 0.567
G/S | 0.964 | 0.666 | 0.647 | 0.588 | 0.727 | 0.317
H/E | 1 | 0.666 | 1 | 0.470 | 1 | 0.430
H/W | 0.370 | 0.533 | 0.411 | 0.411 | 0.421 | 0.354
K/E | 0.365 | 0.466 | 0.529 | 0.647 | 0.584 | 0.414
K/W | 0.372 | 0.466 | 0.235 | 0.764 | 0.358 | 0.477
P/N | 0.607 | 0.4 | 0.647 | 0.294 | 0.718 | 0.287
P/S | 0.431 | 0.266 | 0.470 | 0.411 | 0.510 | 0.276
R/C | 0.352 | 0.666 | 0.705 | 0.470 | 0.695 | 0.347
R/N | 0.638 | 0.533 | 0.705 | 0.705 | 0.753 | 0.425
R/S | 0.521 | 0.666 | 0.529 | 0.823 | 0.639 | 0.706
L | 0.522 | 1 | 0.647 | 1 | 0.766 | 1
M/E | 0.261 | 0.466 | 0.294 | 0.705 | 0.318 | 0.515
M/W | 0.116 | 0 | 0.058 | 0.352 | 0.143 | 0.419
N | 0.379 | 0.4 | 0.588 | 0.470 | 0.602 | 0.308
S | 0.608 | 0.6 | 0.882 | 0.470 | 0.862 | 0.362
T | 0.272 | 0.6 | 0.705 | 0.470 | 0.573 | 0.242
Table 3: All estimated model parameters such as undirected edges (UE),
incoming edges (IE), outgoing edges (OE), incoming weight (IW) and outgoing
weight (OW) for each ward along with disease counts. All values are normalized
with their corresponding maximum

In Fig. 15, the normalized dengue counts, number of undirected edges of each
ward, total incoming weights and total outgoing weights given in (8) and (9)
of each ward are shown.

Figure 14: Normalized number of dengue counts, number of undirected edges from
Fig. 10, number of incoming and outgoing edges from Fig. 11 of each ward
Figure 15: Normalized number of dengue counts, number of edges from Fig. 10,
total incoming weights and total outgoing weights from Fig. 11 of each ward

Several expected and some unexpected inferences are evident from Fig. 14 and
15:

  1. 1.

Dengue counts in some wards (A, E, F/N, G/N, L) are highly correlated with
counts in other wards. This is evident from the high number of undirected
edges of the estimated partial correlation graph incident upon these wards.
While this observation provides some insight into the spread pattern,
information about the direction of the spread is not available.

  2. 2.

Some of the wards (A, R/S, L) act as sources of dengue spread. This is evident
from the high number of outgoing edges of the estimated causality graph
emanating from these wards.

  3. 3.

Actual dengue count can be low for these source wards (A, R/S, L). This is
quite counter-intuitive and one of the key findings of this analysis. A very
high daily commuter flow might explain this phenomenon.

  4. 4.

From Fig. 14 and 15, the number of incoming edges seems to correlate well to
the disease count. This seems to verify that the dengue is spread primarily by
human movement.

##  6 Conclusion

In this paper, we have investigated graphical interaction models of
multivariate time series of counts using a parametric approach. The partial
correlations and causalities between observed multivariate count data are
defined in terms of the partial correlations and causalities between latent
processes. Further a joint MLE with $\ell_{1}$-type regularization is used to
estimate the inverse spectral density matrix. To overcome the computational
difficulties with the resulting mixture distributions, an MCEM algorithm with
$\ell_{1}$-type regularization is proposed. Asymptotic convergence results for
the sequence generated by Algorithm 3.1 are presented and the results are
verified with simulations. Finally, the partial correlation and causality
graphs are estimated for dengue count data observed weekly from each ward of
Greater Mumbai city over six years. Surprisingly some wards seem to act as
epicenters of disease spread even though their absolute disease counts are
relatively low. Such an inference may help correct public health intervention
policies in the future.

Acknowledgments

This work was partially supported by the Science and Research Engineering
Board (Department of Science and Technology, Government of India). The authors
thank Municipal Corporation of Greater Mumbai (MCGM) for providing the dengue
dataset.

## References

  * Alpago et al. (2018) D. Alpago, M. Zorzi, and A. Ferrante.  Identification of sparse reciprocal graphical models.  _IEEE Control Systems Letters_ , 2(4):659–664, Oct 2018.  ISSN 2475-1456.
  * Avventi et al. (2013) E. Avventi, A. G. Lindquist, and B. Wahlberg.  Arma identification of graphical models.  _IEEE Transactions on Automatic Control_ , 58(5):1167–1178, May 2013.  ISSN 0018-9286.
  * Bach and Jordan (2004) Francis R. Bach and Michael I. Jordan.  Learning graphical models for stationary time series.  _IEEE Trans. Signal Processing_ , 52(8):2189–2199, 2004.
  * Brillinger (1981) David R. Brillinger.  _Time Series Data Analysis and Theory_.  Society for Industrial and Applied Mathematics Philadelphia, 1981.
  * Brillinger (1996) David R. Brillinger.  Remarks concerning graphical models for time series and point processes.  _Revista de Econometria_ , 16:1–23, 1996.
  * Brockwell and Davis (2002) Peter J. Brockwell and Richard A. Davis.  _Introduction to Time Series and Forecasting_.  Springer, 2nd edition, mar 2002.  ISBN 0387953515.
  * Casella and Berger (2002) G. Casella and R.L. Berger.  _Statistical Inference_.  Duxbury advanced series in statistics and decision sciences. Thomson Learning, 2002.  ISBN 9780534243128.
  * Cawley and Talbot (2010) Gavin C. Cawley and Nicola L. C. Talbot.  On over-fitting in model selection and subsequent selection bias in performance evaluation.  _Journal of Machine Learning Research_ , 11:2079–2107, 2010\.
  * Chan and Ledolter (1995) K. S. Chan and Johannes Ledolter.  Monte carlo em estimation for time series models involving counts.  _Journal of the American Statistical Association_ , 90(429):242–252, 1995.  ISSN 01621459.
  * Ciccone et al. (2018) V. Ciccone, A. Ferrante, and M. Zorzi.  Robust identification of “sparse plus low-rank” graphical models: An optimization approach.  In _2018 IEEE Conference on Decision and Control (CDC)_ , pages 2241–2246, Dec 2018.
  * Dahlhaus (2000) R. Dahlhaus.  Graphical interaction models for multivariate time series.  _Metrika_ , (51):157–172, 2000.
  * Dahlhaus and Eichler (2003) Rainer Dahlhaus and Michael Eichler.  Causality and graphical models in time series analysis.  _Oxford Stat. Sci. Ser_ , 27, 01 2003.
  * David A. Bessler (2003) Jian Yang David A. Bessler.  The structure of interdependence in international stock markets, journal of international money and finance.  _Journal of International Money and Finance_ , 22(2):261–287, 2003.  ISSN 0261-5606.
  * Dempster (1972) A. P. Dempster.  Covariance selection.  _Biometrics_ , 28(1):157–175, 1972.  ISSN 0006341X, 15410420.
  * Dempster et al. (1977) Arthur Dempster, Natalie Laird, and Donald B. Rubin.  Maximum likelihood from incomplete data via the em algorithm.  _Journal of the Royal Statistical Society. Series B (Methodological)_ , 39:1–38, 01 1977.
  * Durbin and Koopman. (2000) J. Durbin and S. J. Koopman.  Time series analysis of non-gaussian observations based on state-space models from both classical and Bayesian perspectives.  _Journal of the Royal Statistical Society, Ser. B_ , 62:3–56, 2000.
  * Eichler (2006) Michael Eichler.  Fitting graphical interaction models to multivariate time series.  _Proceedings of the 22nd Conference on Uncertainty in Artificial Intelegence_ , 2006.
  * Eichler (2007) Michael Eichler.  Granger causality and path diagrams for multivariate time series.  _Journal of Econometrics_ , 137(2):334 – 353, 2007.  ISSN 0304-4076.
  * Eichler (2012) Michael Eichler.  Graphical modelling of multivariate time series.  _Probability Theory and Related Fields_ , 153(1):233–268, Jun 2012.  ISSN 1432-2064.
  * Fahrmeir and Tutz (1994) L. Fahrmeir and G. Tutz.  _Multivariate Statistical Modelling Based on Generalized Linear Models_.  Springer series in statistics. Springer-Verlag, 1994.  ISBN 9780387942339.  URL https://books.google.co.in/books?id=OionAQAAIAAJ.
  * Fahrmeir and Tutz (2001) L. Fahrmeir and G. Tutz.  _Multivariate Statistical Modeling Based on Generalized Linear Models_.  Springer-Verlag, New York, 2001.
  * Fahrmeir and Wagenpfeil (1997) L. Fahrmeir and S. Wagenpfeil.  Penalized likelihood estimation and iterative Kalman filtering for non-gaussian dynamic regression models.  _Computational Statistics and Data Analysis_ , 24:295–320, 1997.
  * Helmut (2005) Lütkepohl Helmut.  _New introduction to multiple time series analysis_.  Springer Berlin Heidelberg, 2005.
  * Hsiao (1982) Cheng Hsiao.  Autoregressive modeling and causal ordering of economic variables.  _Journal of Economic Dynamics and Control_ , 4:243 – 259, 1982.
  * Hue and Chiogna (2021) Nguyen Thi Kim Hue and Monica Chiogna.  Structure learning of undirected graphical models for count data.  _Journal of Machine Learning Research_ , 22:50–1, 2021\.
  * Jung et al. (2015) A. Jung, G. Hannak, and N. Goertz.  Graphical lasso based model selection for time series.  _IEEE Signal Processing Letters_ , 22(10):1781–1785, Oct 2015.  ISSN 1070-9908.
  * J.Z. Huang and Liu (2006) M. Pourahmadi J.Z. Huang, N. Liu and L. Liu.  Covariance matrix selection and estimation via penalised normal likelihood.  _Biometrika_ , 1(93):85–98, 2006.
  * Li (1994) W. K. Li.  Time series models based on generalized linear models: Some further results.  _Biometrics_ , 50:506–511, 1994.
  * Lounici (2008) Karim Lounici.  Sup-norm convergence rate and sign concentration property of lasso and dantzig estimators.  _Electron. J. Statist._ , 2:90–102, 2008.
  * McLachlan and Krishnan (2008) Geoffrey J. McLachlan and Thriyambakam Krishnan.  _The EM algorithm and extensions_.  Wiley series in probability and statistics. Wiley, Hoboken, NJ, 2. ed edition, 2008.
  * Mehdi Jalalpour and Levin (2015) Yulia Gel Mehdi Jalalpour and Scott Levin.  Forecasting demand for health services: development of a publicly available toolbox.  _Operations Research for Health Care_ , 5:1–9, 2015.
  * Mengersen and Tweedie (1996) K. L. Mengersen and R. L. Tweedie.  Rates of convergence of the hastings and metropolis algorithms.  _The Annals of Statistics_ , 24(1):101–121, 1996\.
  * Park and Park (2019) Gunwoong Park and Sion Park.  High-dimensional poisson structural equation model learning via $\backslash$ell_1-regularized regression.  _Journal of Machine Learning Research_ , 20:95–1, 2019\.
  * Park and Raskutti (2015) Gunwoong Park and Garvesh Raskutti.  Learning large-scale poisson dag models based on overdispersion scoring.  _Advances in neural information processing systems_ , 28, 2015.
  * Park and Raskutti (2017) Gunwoong Park and Garvesh Raskutti.  Learning quadratic variance function (qvf) dag models via overdispersion scoring (ods).  _Journal of Machine Learning Research_ , 18:224–1, 2017\.
  * Pierce and Haugh (1977) David A Pierce and Larry D Haugh.  Causality in temporal systems: Characterization and a survey.  _Journal of econometrics_ , 5(3):265–293, 1977\.
  * Qian et al. (2015) Long Qian, Yi Zhang, Li Zheng, Yuqing Shang, Jia-Hong Gao, and Yijun Liu.  Frequency dependent topological patterns of resting-state brain networks.  _PLOS ONE_ , 10(4):1–19, 04 2015.
  * R.K.Freeland and McCabe (2004) R.K.Freeland and B.P.M. McCabe.  Forecasting discrete valued low count time series.  _International Journal of Forecasting_ , 20:427–434, 2004\.
  * Roberts and Tweedie (1996) G. O. Roberts and R. L. Tweedie.  Geometric convergence and central limit theorems for multidimensional hastings and metropolis algorithms.  _Biometrika_ , 83(1):95–110, 1996.
  * Roberts and Smith (1994) G.O. Roberts and A.F.M. Smith.  Simple conditions for the convergence of the gibbs sampler and metropolis-hastings algorithms.  _Stochastic Processes and their Applications_ , 49(2):207 – 216, 1994.  ISSN 0304-4149.
  * Roy and Dunson (2020) Arkaprava Roy and David B Dunson.  Nonparametric graphical model for counts.  _Journal of Machine Learning Research_ , 21(229):1–21, 2020.
  * Sathish et al. (2019) Vurukonda Sathish, Debraj Chakraborty, and Siuli Mukhopadhyay.  Topology selection using monte carlo expectation and maximization algorithm with l1-type regularization for count data.  pages 6977–6982, 2019.
  * Sims (1972) Christopher A Sims.  Money, income, and causality.  _The American economic review_ , 62(4):540–552, 1972.
  * Songsiri (2010) Jitkomut Songsiri.  _Graphical Models of Time Series: Parameter Estimation and Topology Selection_.  2010\.
  * Songsiri (2013) Jitkomut Songsiri.  Sparse autoregressive model estimation for learning granger causality in time series.  In _2013 IEEE International Conference on Acoustics, Speech and Signal Processing_ , pages 3198–3202. IEEE, 2013.
  * Songsiri (2015) Jitkomut Songsiri.  Learning multiple granger graphical models via group fused lasso.  In _2015 10th Asian Control Conference (ASCC)_ , pages 1–6. IEEE, 2015.
  * Songsiri and Vandenberghe (2010) Jitkomut Songsiri and Lieven Vandenberghe.  Topology selection in graphical models of autoregressive processes.  _J. Mach. Learn. Res._ , 11:2671–2705, December 2010.  ISSN 1532-4435.
  * Songsiri et al. (2009) Jitkomut Songsiri, Joachim Dahl, and Lieven Vandenberghe.  Graphical models of autoregressive processes.  _Convex Optimization in Signal Processing and Communications_ , 01 2009.  doi: 10.1017/CBO9780511804458.004.
  * Tibshirani (1996) R. Tibshirani.  Regression shrinkage and selection via the lasso.  _Journal of the Royal Statistical Society (Series B)_ , 58:267–288, 1996.
  * TjØstheim (1981) Dag TjØstheim.  Granger-causality in multiple time series.  _Journal of Econometrics_ , 17(2):157 – 176, 1981\.
  * Utazi et al. (2018) C. Edson Utazi, Emmanuel O. Afuecheta, and C. Christopher Nnanatu.  A bayesian latent process spatiotemporal regression model for areal count data.  _Spatial and Spatio-temporal Epidemiology_ , 25:25 – 37, 2018.  ISSN 1877-5845.
  * W. J. Granger (1969) Clive W. J. Granger.  Investigating causal relations by econometric models and cross-spectral methods.  _Econometrica_ , 37:424–38, 02 1969.  doi: 10.2307/1912791.
  * Wei and Tanner (1990) Greg C. G. Wei and Martin A. Tanner.  A monte carlo implementation of the em algorithm and the poor man’s data augmentation algorithms.  _Journal of the American Statistical Association_ , 85(411):699–704, 1990.
  * Wright (1921) Sewall Wright.  Correlation and causation.  _J. agric. Res._ , 20:557–580, 1921.
  * Wright (1934) Sewall Wright.  The method of path coefficients.  _The annals of mathematical statistics_ , 5(3):161–215, 1934.
  * Zeger and Qaqish (1988) Scott L. Zeger and Bahjat Qaqish.  Markov regression models for time series: A quasi-likelihood approach.  _Biometrics_ , 44(4):1019–1031, 1988.
  * Zorzi and Sepulchre (2016) M. Zorzi and R. Sepulchre.  Ar identification of latent-variable graphical models.  _IEEE Transactions on Automatic Control_ , 61(9):2327–2340, Sep. 2016.  ISSN 0018-9286.