Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
102,036 characters · 26 sections · 72 citation commands
Latent Gaussian dynamic factor modeling and forecasting for multivariate count time series
\affil[1]{Cornell University} \affil[2]{Friedrich-Alexander-Universit\"at Erlangen-N\"urnberg} \affil[3]{The Pennsylvania State University} \affil[4]{University of North Carolina at Chapel Hill}
\footnotetext{Corresponding author. Email: [email removed]}
\keywordsmsc{Primary 62M10, 62H12; secondary 62H20.}
This work develops theory, estimation, and forecasting methods for dynamic factor modeling of discrete-valued multivariate, possibly high-dimensional time series. Count time series are widespread in the natural, health, social and other sciences, for example, monthly counts of earthquakes or the amount of rainfall above a certain magnitude, daily counts of virus infections over spatial locations, item responses on surveys, and number of followers in social network over time. Mathematically, we consider a $d-$vector time series $\{X_{t}\}_{t\in\mathbb{Z}}=\{(X_{i,t})_{i=1,\ldots,d}\}_{t\in\mathbb{Z}}$, where $X_{i,t}\in\mathbb{N}_0=\{0,1,2,\ldots\}$ and $t\in\mathbb{Z}$ represents time. While the model can handle the range $\mathbb{Z}$, discrete-valued models are commonly applied to counts and the range $\mathbb{N}_0$. Any set of finite discrete values can be represented as a subset of $\mathbb{N}_0$. The primary focus will be on stationary models, though we also discuss the inclusion of covariates and differencing.
In general, modeling time series with discrete (count) values is delicate. In the continuous case, the class of autoregressive moving average (ARMA) models parsimoniously spans all non-deterministic stationary series (by the classical Wold decomposition). In the count setting, the landscape is much less established, with no single class of models dominating in popularity. In fact, researchers have developed numerous methods for constructing stationary count time series. The majority of work on count time series has been devoted to the univariate case. Popular approaches include those based on thinning operators mckenzie1985some,alzaid1993some and the generalized state-space models davis2016handbook, including hidden Markov models (HMMs) macdonald1997hidden, Bayesian dynamic models gamerman2015dynamic. Integer-valued autoregressive conditional heteroskedasticity modeling ferland2006integer,fokianos2009poisson,zhu2011negative is another popular observation-driven approach. Recent reviews of this research area are given in weiss2018introduction,davis2021count and fokianos2021multivariate.
Multivariate, and potentially high-dimensional, count time series have received considerably less attention. A recent popular approach uses generalized linear model (GLM) constructions in high-dimensional settings with component series means, conditionally on the past, depending on their past values or those of the counts themselves akin to vector autoregression (VAR) models; see chen2017multivariate,hall2018learning,mark2018network,mark2019estimating, and fokianos2020multivariate. A variation of the approach is to use dynamic factor model (DFM) constructions instead; see jung2011dynamic,cui2014generalized,wang2018modeling, and brauning2020dynamic. This approach posits conditional distributions of the counts and is convenient for likelihood estimation procedures. In a collection of articles on count time series, karlis2016models surveys relatively recent works on multivariate discrete-valued time series models.
In a recent paper, jia2023latent proposed a new count time series model driven by latent Gaussian time series. The model offers a flexible and most general correlation structure that can accommodate any count marginal distribution. The marginal distribution can exhibit over- or under-dispersion, or have zero-inflation. The autocovariance function (ACVF) of the model is as general as possible for a given marginal distribution, in particular, capable of achieving the most negative pairwise correlations. kong2023seasonal further extended the model to incorporate periodic and seasonal features by replacing the vanilla ARMA with periodic autoregressive moving average (PARMA) and seasonal autoregressive moving average (SARMA) models. The model bins the latent Gaussian series into discrete values and is particularly suitable for data that can be thought of as a discretization of an underlying continuous-valued signal (e.g., when a discrete scale is used for response over a continuous scale).
In this work, we propose a multivariate, possibly high-dimensional extension of jia2023latent where the latent Gaussian series follows a DFM. A recent study by duker2024high considered a similar extension but with the latent Gaussian series following a high-dimensional, sparse VAR model. jia2023latent considered several estimation methods, including the efficient maximum likelihood estimation based on particle approximations of likelihood function via sequential Monte Carlo. The likelihood approximation becomes computationally infeasible with available tools in higher dimensions. Both here and duker2024high, a computationally efficient scheme is employed based on a relationship between second-order properties of the observed and latent Gaussian series. While duker2024high consider a sparse VAR for the latent Gaussian series, we assume the latter to follow a DFM. DFMs are models of choice when the cross-sectional dependence across the variables is strong, for example, manifested through correlations as in our application (see Section (ref) and Figure (ref)). Sparse VARs typically do not have this property and, in fact, there is substantial literature combining the two types of models lin2020regularized,chen2023community. Our approach is also shown to have theoretical guarantees, in the high-dimensional regime. Here, we rely on the concentration bounds for the covariances proved in duker2024high but substantial technical work is still needed to exploit them for the considered DFM and also to show that our factor-based model satisfies certain assumptions in duker2024high. Furthermore, we suggest a novel cross-validation scheme related to model selection, namely, the number of factor series. Though we do not employ likelihood estimation based on particle approximations, we use the latter for forecasting once the model is fitted to data. This still carries an expensive computational cost but it is manageable compared to likelihood inference. Kalman filtering techniques alleviate some of the computational cost. Our modeling and forecasting approaches are illustrated in simulations and on real psychometrics data.
We also note that our model is closely related to another strand of the literature, namely, that on copula models and particularly Gaussian copula models. Indeed, as noted in jia2023latent, for any finite collection of times, our model can be expressed in terms of a Gaussian copula function, depending on our latent model parameters. Related work includes Gaussian copula regression masarotto2012gaussian. Factor-type copula models, not necessarily Gaussian or specifically discrete-valued, were considered by murray2013bayesian,nikoloulopoulos2015factor,kadhem2021factor, and others. These references are far from being exhaustive. What sets this work apart are the time series setting, a new estimation method, the possibility of high-dimensional regime, and theoretical guarantees.
In summary, our main contributions and highlights of the paper are as follows:
The rest of the paper is organized as follows. Section (ref) introduces the latent Gaussian dynamic factor model and establishes relationships between the second-order dependence structures of count and underlying Gaussian models. The estimation procedure is described in Section (ref), followed by Section (ref), providing its theoretical guarantees. Section (ref) concerns forecasting. Numerical experiments can be found in Section (ref), followed by an illustrative application in Section (ref). We close our paper with comments for future work in Section (ref). Finally, Appendix (ref) contains the proofs of our theoretical results and Appendix (ref) contains some details for our forecasting approach.
For a $d-$vector time series $\{X_{t}\}_{t\in\mathbb{Z}}=\{(X_{i,t})_{i=1,\ldots,d}\}_{t\in\mathbb{Z}}$, a latent Gaussian dynamic factor model is defined as follows. For $i=1,\ldots,d$, each component series $X_{i,t}$ at time $t$ is given by
where $F_{i}$ is a cumulative distribution function (CDF), $F_{i}^{-1}(u)=\inf\{v:F_{i}(v)\geq{u}\}$ is its (generalized) inverse, $\Phi(z)$ is the CDF of $\mathcal{N}(0,1)$ distribution, and $Z_{i,t}$ is a zero mean, unit variance, Gaussian stationary series defined below. The CDFs $F_{i}$ are thought to come from parametric families, parameterized by a (possibly vector) parameter $\theta_{i}$. Note that by construction (ref), $\{X_{i,t}\}$ is stationary and its marginal distribution is $F_{i}$. We focus here on discrete distributions $F_{i}$ taking nonnegative integer values $\mathbb{N}_0$. For example, if $F_{i}$ is the CDF of a Bernoulli distribution with parameter $p_{i}$, Bern$(p_{i})$, the model (ref) becomes $X_{i,t}=\boldsymbol{1}_{\{\Phi(Z_{i,t})>1-p_{i}\}}=\boldsymbol{1}_{\{Z_{i,t}>\Phi^{-1}(1-p_{i})\}}$, where $\boldsymbol{1}_{\{\cdot\}}$ is the indicator function. More generally, if $F_i$ is a CDF whose support lies in $\mathbb{N}_0$, for example a Poisson distribution with parameter $\theta_{i}$, then $X_{i,t}$ is represented through $Z_{i,t}$ by
In view of (ref), the model discretizes the continuous-valued series $\{Z_{i,t}\}$ and is particularly natural to use in the context where $\{X_{i,t}\}$ can be thought as resulting from such discretization. In fact, (ref) defines a count random variable $X_{i,t}$ represented through $Z_{i,t}$ that follows any marginal distribution $F_{i}$. However, we shall use the examples of Bernoulli, Poisson, negative binomial, and categorical marginal distributions for illustration throughout the paper. For example, for categorical marginal counts, the infinite sum in (ref) reduces to a finite sum by the number of categories, say $N$, where $C_{i,N} = 1$. Whereas the Poisson distribution with a large parameter is close to Gaussian and the latent $\{Z_{i,t}\}$'s are effectively observed, note that the Bernoulli case lies at the other extreme and is expected to be most difficult to deal with in our tasks.
We are interested in the scenario where the underlying Gaussian series $\{Z_{i,t}\}$ obeys a DFM. More specifically, we suppose that the $d-$vector time series $\{Z_{t}\}_{t\in\mathbb{Z}}=\{(Z_{i,t})_{i=1,\ldots,d}\}_{t\in\mathbb{Z}}$ satisfies
where $\Lambda$ is a $d\times{r}$ loadings matrix, and $\varepsilon_{t}$ are $\mbox{i.i.d. }\mathcal{N}(0,\Sigma_{\varepsilon})$ random $d-$vectors (independent of $Y_{t}$'s) and $r-$vector factor series $\{Y_{t}\}_{t\in\mathbb{Z}}=\{(Y_{k,t})_{k=1,\ldots,r}\}_{t\in\mathbb{Z}}$ follows a causal (stable) stationary VAR model of order $p$, VAR($p$), given by
where $\Psi_{1},\ldots,\Psi_{p}$ are $r\times{r}$ matrices and $\eta_{t}$ are $\mbox{i.i.d. }\mathcal{N}(0,\Sigma_{\eta})$ random $r-$vectors. The VAR model (ref) is flexible to capture temporal dependence from a practical standpoint. Note that the DFM (ref) is in the so-called static form. The generalized DFMs where (ref) includes lags of $Y_t$ forni2000generalized,forni2005generalized go beyond the scope of this work. Gaussianity is assumed for the various components of (ref) and (ref) given $Z_{i,t}$ is Gaussian in (ref). The factor structure in (ref) imposes dependence of $Z_{i,t}$ and hence also $X_{i,t}$ across $i=1,\ldots,d$.
Note that the unit variance of $Z_{i,t}$ is assumed in (ref). For general $Z_{i,t}$, one can standardize it to have unit variance. More generally, the ACVF $\Sigma_{Z}(h) = {\mathbb E}[ Z_{t+h}Z_{t}']$ of $\{Z_t\}$ at lag $h$ can similarly become ACF $R_Z(h)$ as
We use $R_{Z}(h)$ for the rest of the analysis so the unit variance assumption for $Z_{i,t}$ is made throughout.
Our estimation procedure is based on the following property of the model (ref). It is known pipiras2017long that, for any $i,j=1,\ldots,d$,
or, in short, and entry-wise,
where $L_{ij}:[-1,1]\mapsto[-1,1]$ are functions to be referred to as link functions (and $L$ as a link function). Furthermore, $L_{ij}$ depends only on the CDFs $F_{i}$ and $F_{j}$ and can be expressed as described below. As $F_i$ depends on parameter $\theta_i$, we shall sometimes write $L_{\theta_i,\theta_j}$ instead of $L_{ij}$ to indicate this dependence.
For $k=0,1,\ldots$, let $H_{k}(z)=(-1)^{k}e^{z^{2}/2}(d^{k}e^{-z^{2}/2}/dz^{k})$ be the Hermite polynomial of order $k$ and
be the corresponding Hermite coefficient of the function $G_{i}(z)$ in (ref), so that $G_{i}(z)=\sum_{k=0}^{\infty}g_{i,k}H_{k}(z)$. For $G_{i}(z)$ associated with the CDF $F_i$ on nonnegative integers, jia2023latent showed that
where $Q_{i,n}=\Phi^{-1}(C_{i,n})$ and $C_{i,n}={\mathbb P}(X_{i,t}\leq{n})=F_{i}(n)$. When $Q_{i,n}=\pm\infty$ for $C_{i,n}=0$ or $1$, the summand $e^{-Q_{i,n}^{2}/2}H_{k-1}(Q_{i,n})$ is interpreted as zero. For example, for $F_{i}=$ Bern$(p_{i})$, $g_{i,k}=e^{-\Phi^{-1}(1-p_{i})^{2}/2}H_{k-1}(\Phi^{-1}(1-p_{i}))/(k!\sqrt{2\pi})$. Similarly, $g_{i,k}$ can be computed when $F_i = \textrm{Pois}(\lambda_i)$ but (ref) will have infinitely many terms. However, the number of terms will be finite and small in practice. This is because the Poisson distribution is light-tailed so that $Q_{i,n}$ is indistinguishable from 1 numerically even at moderate $n$.
The link functions $L_{ij}$ can now be expressed as
Under mild assumptions, they can be shown to be monotonically increasing on the interval $(-1,1)$ jia2023latent with values in $(L_{ij}(-1),L_{ij}(1))$. Note that $L_{ij}(0)=0$ regardless of the marginal distribution. The quantities $\rho_{+,ij} := L_{ij}(1)$ and $\rho_{-,ij} := L_{ij}(-1)$ are given by
for $Z=\mathcal{N}(0,1)$. When $i=j$, $\rho_{+,ij}=1$ but usually $\rho_{-,ij}>-1$. As noted in jia2023latent, $\rho_{+,ij}$ and $\rho_{-,ij}$ are the largest and smallest correlations that two dependent count variables with marginals $F_i$ and $F_j$ can achieve; by (ref), they are achieved with construction $G(Z)$ and hence within our considered model. For example, when $F_{i}=$ Bern$(p_{i})$, it can be shown by using (ref) that
Since the link functions $L_{ij}$ are monotonically increasing, the inverse link functions can be defined as $L_{ij}^{-1}:[\rho_{-,ij},\rho_{+,ij}]\mapsto[-1,1]$. We will discuss the numerical calculation of the inverse $L_{ij}^{-1}$ in Section (ref) below. Thus, (ref) and (ref) imply that
or, in short and entrywise,
We will exploit (ref) in estimation in the next section.
The relation (ref) suggests a natural estimation procedure for the model parameters of the latent series $\{Z_t\}$. Indeed, recall that the component function $L_{ij}$ of $L$ is defined through the marginal CDFs $F_i$ and $F_j$ which depend on the marginal parameters $\theta_i$ and $\theta_j$, respectively. As expanded upon in Section (ref) of Supplemental Material, the marginal parameters $\theta_i$ can be estimated marginally from the stationary data $X_{i,1},\ldots,X_{i,T}$. This leads to estimates $\widehat{\theta}_i$ and $\widehat{\theta}_j$, and hence $\widehat{L}_{ij}$ and $\widehat{L}$. Computation of the inverse link functions $L_{ij}^{-1}$ (or $\widehat{L}_{ij}^{-1}$) is discussed in Section (ref) below. On the other hand, $R_{X}(h)$ in (ref) can be estimated through sample ACF $\widehat{R}_X(h)$ of the data. Substituting the sample quantities $\widehat{L}$ and $\widehat{R}_X(h)$ into the right-hand side of (ref) leads to estimate $\widehat{R}_Z(h)$, which characterizes the second-order properties of the latent process $\{Z_t\}$. For the latent VAR-DFM considered in this work, the model parameters can be estimated from the model second-order properties and hence $\widehat{R}_Z(h)$ as described in Section (ref). Taken together, the approach allows estimating the marginal and latent series parameters of our model.
We describe here estimation of the parameters $\Lambda,\Psi_{1},\ldots,\Psi_{p},\Sigma_{\varepsilon}$ and $\Sigma_{\eta}$ of the latent series. We assume here that the marginal CDFs $F_{i}$ and hence the link functions $L_{ij}$ are known. In practice, the link functions can be estimated from marginal stationary data as noted above. We also assume here that $r$ and $p$ are known. Their estimation is descibed in Section (ref) below and Section (ref) of Supplemental Material. Note that the relation (ref) allows one to estimate $R_{Z}(h)=\Sigma_{Z}(h)$, $h=0,\ldots,p$, as
where $\widehat{R}_{X}(h)$ is the sample matrix ACF of the data $X_{1},\ldots,X_{T}$. $\widehat{R}_{X}(h)$, $h\neq0$, are not necessarily symmetric. Similarly, $\widehat{R}_Z(0)$ is symmetric but not necessarily nonnegative definite. The estimated covariance $\widehat{R}_Z(0)$ of the latent Gaussian process will be used in forecasting described in the next section. If needed, we employ a small positive shift of the eigenvalues to make $\widehat{R}_Z(0)$ nonnegative definite. However, we do not shift eigenvalues through the estimation procedure.
Since $\{Z_t\}$ follows the dynamic factor model (ref) and (ref), the loadings matrix $\Lambda$, $\Sigma_{Y}(0)$, and $\Sigma_{\varepsilon}$ can be estimated through principal component analysis (PCA). More specifically, since the dynamic factor model (ref) implies
it is natural to estimate $\Lambda\Sigma_{Y}(0)\Lambda'$ as an $r-$rank approximation of $R_{Z}(0)$, and take $\Sigma_{\varepsilon}$ as the approximation error. We thus proceed as follows. Consider the covariance $\widehat{R}_Z(0)$ estimated by (ref). Let $\widehat{R}_{Z}(0)=\widehat{U}\widehat{E}\widehat{U}'$ be the eigendecomposition with $\widehat{E}=\operatorname{diag}(\widehat{e}_{1},\ldots,\widehat{e}_{d})$ consisting of ordered eigenvalues $\widehat{e}_{1}\geq\ldots\geq\widehat{e}_{d}$, and $\widehat{U}=(\widehat{u}_{1},\ldots,\widehat{u}_{d})$ being the orthonormal eigenvector matrix. Setting $\widehat{U}_{r}=(\widehat{u}_{1},\ldots,\widehat{u}_{r})$ and $\widehat{E}_{r}=\operatorname{diag}(\widehat{e}_{1},\ldots,\widehat{e}_{r})$ as $d\times r$ and $r \times r$ matrices, respectively, a rank-$r$ approximation of $\Sigma_{Z}(0)$ can be taken as
The relations (ref) and (ref) suggest estimating the $d \times r$ loadings matrix $\Lambda$ and the $r \times r$ covariance matrix $\Sigma_{Y}(0)$ of the factor series as
The choice (ref) identifies $\Lambda,\Sigma_{Y}(0)$ up to a non-singular $r \times r$ transformation, assuming that
The identifiability condition (ref) is common in factor models doz2011two,bai2013principal. Another identifiability condition used in DFMs proposed by bai2015identification is to make the first $r \times r$ block of the loadings matrix be identity, that is,
where $\Lambda_{2}$ is $(d-r)\times r$. Note that with the convention (ref) above,
The relations (ref) and (ref) also suggest setting
to define both $\widehat{\Sigma}_{Y}(0)^{1/2}$ and $\widehat{\Lambda}_{2}$. Either (ref) or (ref) lead to estimators $\widehat{\Lambda}$ and $\widehat{\Sigma}_{Y}(0)$. The estimator $\widehat{\Sigma}_{\varepsilon}$ can now be defined as
We consider the first identifiability condition (ref) for the various simulation settings in Section (ref), while the second condition (ref) is used in Section (ref) with the application.
Note also that DFM (ref) yields
Setting $\widehat{R}_{Z}(h)=L^{-1}(\widehat{R}_{X}(h))$ as in (ref) naturally suggests the estimators
Alternatively, the estimators (ref) also solve
where $\| A \|_{F}^2 = \sum_{i,j=1}^d |A_{ij}|^2$ denotes the Frobenius norm for a matrix $A =(A_{ij})_{i,j=1,\dots,d} \in \mathbb R^{d \times d}$. Note that the estimator $\widehat{\Sigma}_{Y}(h)$ can be obtained regardless of the identifiability condition, either (ref) or (ref). Having these estimators, the Yule-Walker equations can now be used to obtain the rest of the required estimators $\widehat{\Psi}_{1},\ldots,\widehat{\Psi}_{p}$ and $\widehat{\Sigma}_{\eta}$. That is, $\widehat{\Psi}_{1},\ldots,\widehat{\Psi}_{p}$ solve the system of matrix linear equations
and $\widehat{\Sigma}_{\eta}$ is
In estimation, the link functions $L_{ij}$ entering $L$ in (ref) are replaced by their sample counterparts $\widehat{L}_{ij}$. We expand here on how $\widehat{L}_{ij}$ and their inverse $\widehat{L}_{ij}^{-1}$ are computed.
Link functions $L_{ij}$ in (ref) are defined in terms of the Hermite coefficients $g_{i,k},g_{j,k}$ which are obtained from marginal distributions $F_{i},F_{j}$ through (ref). If $F_{i},F_{j}$ are characterized by model parameters $\theta_i,\theta_j$, then $g_{i,k}=g_{i,k}(\theta_i),g_{j,k}=g_{j,k}(\theta_j)$ are functions of these parameters. The parameters $\theta_i,\theta_j$ can be estimated through observations for marginal $i,j$ as discussed in Section (ref), leading to estimated Hermite coefficients $\widehat{g}_{i,k}=g_{i,k}(\widehat{\theta}_i),\widehat{g}_{j,k}=g_{j,k}(\widehat{\theta}_j)$ and link function coefficients $\widehat{\ell}_{ij,k}$. Link function estimates are then $\widehat{L}_{ij}(u)=\sum_{k=1}^K \widehat{\ell}_{ij,k}u^k$ for large $K$, say $K=100$ or more. To simplify the notation, we write $L_{ij}$ for $\widehat{L}_{ij}$, and consider the computation of $L_{ij}^{-1}$ next.
The idea to calculate $L_{ij}^{-1}$ is as follows. Partition the interval $[-1,1]$ into $u_0<u_1<\dots<u_M$ that satisfies $u_{0}=-1$, $u_{M}=1$, and set
so that one has the points $(v_{m},L_{ij}^{-1}(v_{m}))=(v_{m},u_{m})$ on the curve $L_{ij}^{-1}(v)=u$. The value of $L_{ij}^{-1}(v)$ for other points $v$ can then be obtained through piecewise linear interpolation. In addition, we use finer grids for $u$ near $\pm1$ while wider grids are used around $u=0$. From a numerical standpoint, the $M+1$ points $(v_m,u_m)$ satisfy $v_m = L_{ij}^{-1}(u_m)$, and the interpolation of $L_{ij}^{-1}$ is defined as
for $v\in(v_{m},v_{m+1})$, $m=0,1,\ldots,M$, with $M$ pieces of the linear spline polynomials
Figure (ref) depicts $L_{ij}(u),L_{ij}^{-1}(v)$ and its interpolation $\widetilde{L}_{ij}^{-1}(v)$ for four representative marginal count distributions: Bernoulli, categorical, Poisson, and negative binomial. For example, the Bernoulli case considers a pair of $F_i=\textrm{Bern}(p_i)$ and $F_j=\textrm{Bern}(p_j)$ with four different choices of combinations $(p_i,p_j)$. Several combinations of parameters for other types of distributions, including categorical, Poisson, and negative binomial distributions, are also considered. As seen from the plots, the inverse link function $L_{ij}^{-1}$ obtained by flipping the axes and numerical inverse $\widetilde{L}_{ij}^{-1}$ obtained through the linear interpolation are nearly indistinguishable. These combinations of marginal distributions are used in the simulation study in Section (ref).
In practice, the number of factor series $r$ is unknown and needs to be chosen for the model estimation. We consider here several methods for this task based on available approaches and also introduce cross-validation schemes tailored for our model.
One practical approach is to examine a scree plot of the eigenvalues of $\widehat{\Sigma}_{Z}(0) = \widehat{L}^{-1}(\widehat{\Sigma}_{X}(0))$ where the presence of a “knee” suggests the value of $r$. More formally, one can design an algorithm that determines the “knee.” For example, onatski2010determing suggested the Edge Distribution (ED) estimator, whereby
where $\widehat{e}_1 \geq \widehat{e}_2 \geq \ldots \geq \widehat{e}_{d}$ are the ordered eigenvalues of $\widehat{R}_Z(0)$ and $\delta$ is calibrated through the algorithm described in Section 4 of that paper.
Alternatively, one could rely on information criteria (IC). That is, from $q=1,\ldots,r_{\max}$, we choose $q$ as the estimate $\widehat{r}$ that minimizes
where $\widehat{\Sigma}_{\varepsilon}(q)$ is the estimator in (ref) from the rank-$q$ approximation of $\Sigma_{Z}(0)$ in (ref), and $g_{i}(d,T)\to 0$ and $\min(d,T)g_i(d,T)\to \infty$ as $d,T \to \infty$. The recommended choices of the penalty functions $g_{i}(d,T)$ are
which corresponds to $\textrm{IC}_{p1}(r)$--$\textrm{IC}_{p3}(r)$ studied by bai2002determining.
We now propose a block cross-validation (BCV)-based rank selection. The idea goes back at least to browne1989single, and continues being utilized in psychometrics haslbeck2022estimating. Differences from the setting considered here are that our observations are serially correlated and that our factor model is latent. To account for temporal dependence, we do not partition observations randomly but rather into equally sized consecutive blocks. The latent nature of the factor model will be dealt with by exploiting the idea in Remark (ref). We thus fold $\{X_t\}$ along the time into $B$ blocks. The superscript $(b)$ will refer to the $b$th block, to be used for test data. The superscript $(-b)$ will refer to the $b$th block being excluded, to be used for training data. Let $\widehat{R}_{Z}^{(b)}(0) = \widehat{L}^{-1}(\widehat{R}_{X}^{(b)}(0))$ be the sample matrix ACF of the latent Gaussian series at lag 0, computed from the sample matrix ACF of the observations from the $b$th block, substituted into the inverse link function. Similarly, one can compute $\widehat{R}_{Z}^{(-b)}(0) = \widehat{L}^{-1}(\widehat{R}_{X}^{(-b)}(0))$ that excludes the $b$th block. Then, for each candidate rank $q$ on a grid $q=1,\ldots,r_{\max}$, the mean square error (MSE) of BCV is defined as
where $\widehat{R}_{Z}^{(-b,q)}(0)$ is the rank-$q$ approximation of $\widehat{R}_{Z}^{(-b)}(0)$ plus a diagonal matrix from the estimated covariance matrix of the innovations of the factor model. The PCA-based estimation procedure described in Section (ref) can be used. The minimizer $q$ of the MSE is chosen as the estimate of the number of factor series $r$.
Instead of our PCA-based estimation, other estimation approaches can and have been used as well. For example, in the minimum residual factor analysis harman1966factor, $\widehat{R}_{Z}^{(-b,q)}(0)$ is the rank-$q$ approximation of $\widehat{R}_{Z}^{(-b)}(0)$ obtained by minimizing the difference between $\widehat{R}_{Z}^{(-b)}(0)$ and the sum of $\widehat{R}_{Z}^{(-b,q)}(0)$ and a diagonal matrix, the latter accounting for the variances of the error terms. This estimation approach is quite popular in factor analysis, with the estimation procedure implemented through R package psych by revelle2023psych. Other estimation approaches for factor analysis can be found in bertsimas2017certifiably.
The selection of lag order $p$ in (ref) can be similarly conducted by using cross-validation, which is explained in Supplemental Material. Note that determining the number of factors can be performed regardless of the lag order of the factor series. We thus recommend choosing $r$ first, followed by selecting $p$.
In this section, we investigate theoretical properties of our estimators for the latent factor model (ref). Our results are derived for the PCA estimators of the latent factor model in Section (ref). The proofs are based on duker2024high who provide a theoretical foundation to investigate the model (ref) for general latent Gaussian processes. We base our investigations on results in doz2011two who prove that the PCA estimators are consistent when the series following a factor model is observed and also on concentration results of bounded functions of Markov chains in fan2021hoeffding. We start with briefly recalling our estimation procedure and introduce some additional notation in Section (ref). We collect assumptions in Section (ref) and then state our main results in Section (ref). Section (ref) concludes with a particular case satisfying the assumptions made in Section (ref).
We first recall and supplement some of the quantities and their estimated counterparts given in Section (ref). Recall the function $L$ from (ref) and introduce its non-standardized counterpart
where $g_{i,k}$ are the Hermite coefficients defined in (ref). Then, using (ref), (ref) can also be written as
with $R_{Z}(h)$ as in (ref). We write $C_{n}(\theta_{i}):=C_{i,n}$ and $Q_{n}(\theta_{i}):=Q_{i,n}$ to emphasize dependence on $\theta_{i}$. An estimator of $\ell$ is written as $\widehat{\ell}$ and computed by replacing $C_{i,n}$ with $\widehat{C}_{i,n} = C_{n}( \widehat{\theta}_{i})$ in $Q_{i,n} = \Phi^{-1}(C_{i,n})$ in (ref). Estimation can then be conducted using (ref) and a regular autocovariance estimator $\widehat{\Sigma}_{X}(h)$ for $\Sigma_{X}(h)$ based on our observed count series $\{X_t\}$ such that
We also introduce $\bm{R}_{Z} = (R_Z(r-s))_{r,s=1,\dots,p}$ and $\bm{\Sigma}_{X} = (\Sigma_X(r-s))_{r,s=1,\dots,p}$ with $p$ denoting the lag order of the latent factors. With a slight abuse of notation, we write
We now move on to estimation of the latent model parameters in (ref) and (ref). Based on the decomposition (ref), we get
On the population level, we write $E_r$ for the diagonal matrix whose diagonal entries are the eigenvalues of $\Lambda' \Lambda$ and $Q_r$ for the matrix of a set of unitary eigenvectors associated with $E_r$. Set further
and recall from (ref) the relations
In order to estimate the transition matrices in the VAR($p$) model (ref), we write (ref) as a $p r$-dimensional VAR($1$) model, that is,
On the population level, we introduce $\Sigma_{Y}(h) = {\mathbb E}[ Y_{t+h} Y_{t}']$, $\bm{\Sigma}_Y^{(p)} := (\Sigma_Y(r-s))_{r,s=1,\dots,p}$ such that $\bm{\Sigma}_{QY}^{(p)} := (Q\Sigma_Y(r-s)Q')_{r,s=1,\dots,p} = \bm{Q} \bm{\Sigma}_{Y}^{(p)} \bm{Q}'$ with $\bm{Q} = I_p \otimes Q_r$. Consistency results are derived below up to transformation with the orthogonal matrix $Q_r$. Since the factors are defined up to a pre-multiplication by an invertible matrix, we choose $\Sigma_Y(0) = I_r$ and maintain this assumption throughout the paper.
Using the introduced notation, the equation in (ref) can be written as
with $S_1 =
$, $S_2 =
$, $0_{p}$ being a $p-$dimensional column vector with all entries set to zero. Furthermore, we can infer from (ref) that
is the solution of (ref) and we get $\bm{\Psi} = S_2 \bm{\Sigma}_{QY}^{(p+1)} S'_1 (\bm{\Sigma}_{QY}^{(p)})^{-1}$ on the population level. We further have $ \widehat{\Sigma}_{Y}(h) = \widehat{E}_r^{-1/2} \widehat{U}'_r \widehat{R}_Z(h) \widehat{U}_r \widehat{E}_r^{-1/2} $ and
We work with two sets of assumptions. The first assumption set applies to Proposition 3.1 in duker2024high, which relates the probability of how much $\bm{\widehat{R}}_Z$ deviates from $\bm{R}_{Z}$ to the analogous probabilities for the respective quantities of the observed series $\{ X_t \}$.
The constant $\bm{c}_{Z}$ in Assumption (ref) controls the dependence of the series $\{Z_t\}$ and appears in the non-asymptotic bound derived in duker2024high which are used in our work. The assumption is expected to hold for some $\bm{c}_{Z}$ for the processes considered both here and in duker2024high. Assumption (ref) is shown to hold for several common count distributions under Assumption (ref); see Appendix E in duker2024high. Note that we require our moment conditions to hold uniformly in a neighborhood around $\theta_{i}$. This allows us to infer finiteness on a compact subset of the parameter space of $\theta_{i}$. The second set of assumptions is needed to derive consistency results for estimators of the latent factor model. The asymptotics described below are understood as $d,T\to\infty$.
Assumptions (ref) and (ref) are crucial to ensure consistent estimation of the transition matrices of the latent factors. We point out that the rate $\sqrt{ \frac{\log(dT)}{T}}$ might not be optimal. Results in doz2011two, where estimation for the latent factors is done assuming that the series following a factor model is observed, suggest that a rate of $\frac{1}{\sqrt{T}}$ might be possible. See also Remark (ref) below. We prove that Assumptions (ref) and (ref) are satisfied under mild assumptions on the marginal distributions in Section (ref).
We remark also that Assumptions (ref) and (ref) are the analogues of Assumptions C.1 and C.2 in duker2024high. Assumptions C.1 and C.2 in duker2024high ensure consistent estimation of $R_Z$ when $\{Z_t\}_{t\in\mathbb{Z}}$ follows a $d-$dimensional VAR($p$) model. Those assumptions are concentration inequalities, suggesting non-asymptotic bounds for the deviation between estimated and true quantities. We use here a representation in terms of stochastic boundedness which is a weaker assumption but sufficient for our purposes since we aim to recover the results in doz2011two who also state their results in terms of stochastic boundedness.
The third set of assumptions is on the latent factor model. The assumptions below and their interpretation are essentially the same as in doz2011two. The notation $\| \cdot \|$ stands for the spectral norm of a matrix, and $\| \cdot \|_{\max}$ for the maximum of its elements (in absolute value). For a symmetric matrix $M$, $\lambda_{\min}(M)$ and $\lambda_{\max}(M)$ will refer to its minimum and maximum eigenvalues.
Our main results state that certain quantities of the latent model can be consistently estimated using the relationship (ref). This will yield consistency of the loadings matrix and the transition matrices in the VAR$(p)$ representation of the latent factors. The proofs of the stated results are in Appendix (ref).
The following results are analogous to Lemma 2(i), and Proposition 3 in doz2011two. While doz2011two assume that the observed series follows a factor model, we recover their results even though we observe a count series with latent factor model structure. Note that our rate differs from that derived by doz2011two by a $\sqrt{\log(dT)}$ term; see also Remark (ref) below.
Under some additional assumptions on the marginal distribution of the observed series in (ref), one can verify Assumptions (ref) and (ref).
Define the function $G: \mathbb R^d \to \mathbb R^d$ with $z \mapsto G(z)$ such that our multivariate count model can be written as
Assumptions (ref) and (ref) cover important cases in the modeling of count time series including Bernoulli marginals. We expect that the boundedness of $G$ in Assumption (ref) can be relaxed but are constrained by concentration results currently available in the literature.
Our proposed verification of Assumptions (ref) and (ref) employs concentration results for Markov chains. In Appendix (ref), we define a reversible Markov chain $\{ \mathcal{V}_t \}$ with Gaussian transition kernel $K(x,\cdot)$ with mean $\bm{\Psi}_{\mathcal{V}} x$ and covariance matrix $\Sigma_{\xi}$, where, with $\widetilde{q} = d(p-1)$,
with
We assign more meaning to the above quantities in the proofs; see Appendix (ref).
The used concentration results for Markov chains are expressed in terms of $\lambda_{\operatorname{r}}:=\lambda_{\operatorname{r}}(P)$, the rightmost value of the spectrum $[-\lambda, \lambda]$ of the Markov operator $P$ induced by the transition kernel $K$. We refer to $1-\lambda_{\operatorname{r}}$ as the right spectral gap of the Markov chain. For the reversible Markov chain $\{ \mathcal{V}_t \}$ with Markov operator $P$ and stationary distribution $\pi$, $\lambda_{\operatorname{r}}$ is defined as
where $\langle h_1,h_2 \rangle_{\pi} = \pi(h_1 h_2)$ and $\| h \|_{\pi}$ is its induced norm.
Assumption (ref) is expected to be satisfied under quite general assumptions. We refer to Section (ref) for a discussion on Assumption (ref).
The following lemma formalizes that the described setting is sufficient for Assumptions (ref) and (ref) to be satisfied.
A fitted latent Gaussian dynamic factor model in (ref)--(ref) can naturally be used to forecast the series $X_{t}$. We will show how this can be carried out for fixed model parameters. For example, the model parameters can be obtained through estimation (in which case our forecast will not reflect any uncertainty from estimation error). More specifically, for a given $t$ (typically, $t=T$, the sample length), we are interested in the distribution of
where $h=1,2,\ldots$, the vertical bar indicates the conditioning and $x_{1},\ldots,x_{t}$ are the observed values of $X_{1},\ldots,X_{t}$. This is equivalent to finding
for arbitrary function $V$, where the subscript $x_{1:t}$ in $\mathbb{E}_{x_{1:t}}$ refers to the conditioning on $\{x_{1},\ldots,x_{t}\}$ as in (ref). For example, suppose we have a $d-$dimensional vector $x$. With $V(x) = V((x_{i})_{i=1,\ldots,d})=\boldsymbol{1}_{\{x_{i}=n_{i},\ i=1,\ldots,d\}}$ and $d-$dimensional integer $n$, the quantity (ref) becomes $\mathbb{P}_{x_{1:t}}(\widehat{X}_{t+h|t}=n)$, which is of primary interest in forecasting $X_{t+h}$, when the components of $X_t$ are integer-valued.
Let $\widehat{Z}_{t+h|t}=\widehat{Z}_{t+h}(Z_{1:t})=H_{t1}^{(h)}Z_{t}+\ldots+H_{tt}^{(h)}Z_{1}$ be the $h$-step-ahead linear prediction of $Z_{t+h}$ from $Z_{1:t}$. Define $\widehat{R}_{t+h|t} = \mathbb{E}[(Z_{t+h}-\widehat{Z}_{t+h|t})(Z_{t+h}-\widehat{Z}_{t+h|t})']$ as the corresponding covariance matrix of prediction error of $Z_{t+h}$. One can compute the latent prediction $\widehat{Z}_{t+h|t}$ by using Kalman recursions as recalled in Appendix (ref). This exploits the state-space structure of the model and is more efficient computationally than a direct application of e.g. Durbin-Levinson algorithm. Then, the quantity (ref) can be expressed as
where
jia2023latent. The right-hand side of (ref) will be approximated through a Monte Carlo scheme below; direct numerical calculation of underlying integrals is too cumbersome.
Note that conditioning on $x_{1:t}$ does not determine the exact path of $z_{1:t}$. Indeed, recall from (ref) that for $i=1,\ldots,d$,
where $C_{i,n}=\mathbb{P}(X_{i,t}\leq n)$ and $A_{i,x_{i}}=\left(\Phi^{-1}(C_{i,x_{i}-1}),\Phi^{-1}(C_{i,x_{i}}) \right]$. That is, each entry of the realization of $X_{t}$ is determined by the range of the corresponding entry of $Z_t$ at each time $t$. For $d-$dimensional observations, one has
where $A_{x_t}=A_{1,x_{1,t}}\times \ldots \times A_{d,x_{d,t}}$. The notation $A_{i,x_{i}}$ and $A_x$ will be used below. In the Bernoulli case, for example,
The Bernoulli marginal distribution thus has the largest ranges for $Z_t$. In a Monte Carlo approximation of (ref), one will be generating $Z_{i,t}\in A_{i,x_{i,t}}$ and producing their forecast $\widehat{Z}_{i,t+h}$.
The following presentation extends that of jia2023latent to the multivariate setting relevant to the Monte Carlo approximation problem. The quantity (ref) is known to be well approximated through Sequential Monte Carlo (SMC) by generating “particles” over time $t$ doucet2001introduction,doucet2009tutorial. The main difference from jia2023latent is that the model here has two latent processes $\{Z_t\}_{t\in\mathbb{Z}}$ and $\{Y_t\}_{t\in\mathbb{Z}}$. To deal with the two latent processes, we use Kalman recursions to forecast and update the latent process $\{Y_t\}$ and approximate the distribution of $\{Z_t\}$ conditioning on $\{X_t\}$. This approach is called Rao-Blackwellization and the method is adapted from the partially observed Gaussian state-space models andrieu2002particle,briers2010smoothing.
The following is the SMC algorithm for particle filtering to generate particles $\{\widetilde{Z}_t^{(k)}\}_{t=1,\ldots,T}$, $k=1,\ldots,N$, over time, whose weighted average then approximates (ref). The particle $\{\widetilde{Z}_{t}^{(k)}\}_{t=1,\ldots,T}$ can be regarded as the realization of the underlying latent process. Additionally, we add a resampling step which is often used in sequential Monte Carlo algorithms. We will explain the necessity of resampling below.
Sequential Importance Sampling and Resampling (SIS/R): Set the initial importance weights $w_{0}^{(k)}=1$ for all $k$, initialize $\widetilde{Y}_{0|0}^{(k)} \sim \mathcal{N}(0,\widetilde{Q}_{0|0})$, where $\widetilde{Q}_{0|0}=\textrm{Var}(Y_0)$. For $p=1$, $\widetilde{Q}_{0|0}$ is approximated by $\sum_{m=0}^M \Psi^m \Sigma_{\eta} (\Psi')^m$ for large $M$. Then, recursively over $t=1,\ldots,T$, do the following steps: For each $k=1,\ldots,N$:
Finally, the SMC approximation of (ref) becomes
where $\widehat{Z}_{t+h|t}$, $h\geq1$, are computed through forecasting step in the Kalman recursions. See equation (25) in jia2023latent for the justification of an analogous approximation.
The SMC is known to suffer from the so-called weight degeneracy of particles snyder2008obstracles which occurs when the variance of normalized weights becomes inflated. The latter happens and becomes worse as the sample size increases. To overcome this, it is suggested to remove the particles with small weights. By following doucet2009tutorial, we resample only when the effective sample size (ESS) exceeds $N/2$, as a rule of thumb, for the criteria of resampling, where the ESS is defined as
More specifically, we resample particles by following systematic resampling. That is, sample $U_1 \sim \mathcal{U}(0,1/N)$ and set $U_k = U_1 + \frac{k-1}{N}$, $k=2,\ldots,N$. Then, compute
with $\sum_{j=1}^0=0$ as a convention. This is used in Step 3 of the SIS/R algorithm above. Alternatively, one can use categorical resampling, which is resampling particles by regarding $\{\widetilde{w}_t^{(k)}\}$ as a categorical probability distribution. Many other resampling methods exist (see douc2005comparison for more information).
The computation burden for sequential Monte Carlo sampling is substantial. The majority of the cost is due to sampling (doubly) truncated multivariate Gaussian random variables $\{\xi_{t}^{(k)}\}$ in (ref), and the fact that the algorithm runs for $t=1,\ldots,T$. Currently, we implement sampling through the R package TruncatedNormal developed by botev2017normal. But a significant improvement in the computation speed for generating truncated multivariate normal random variables is not expected.
To reduce the computational cost, we note that the covariance matrices $\widehat{R}_{t|t-1}$ of prediction error typically converge within a few steps. This is due to a similarly quick convergence of the covariance matrix $\widehat{Q}_{t|t-1}$ of prediction error of the factor series in (ref), and the covariance matrix $\widetilde{Q}_{t|t}$ and Kalman gain $K_t$ described in (ref).
From the pair of covariance matrices in (ref) and (ref), one has the recursive equation
with the given initial condition $\widehat{Q}_{1|0}=\Psi\widetilde{Q}_{0|0}$. The covariance matrices of the prediction error converge to a positive definite matrix $Q$ satisfying the discrete algebraic Riccati equation (DARE),
One has a similar equation for the covariance matrices $\widetilde{Q}_{t|t}$. It is the convergence to these equations that happens within a few time steps substantially shorter than the length of observations $T$. Since the purpose of the SIS/R algorithm is to obtain the importance weights $\{\widetilde{w}_{T}^{(k)}\}$ along the particles $\{Z_{1:T}^{(k)}\}$ and we presume the stable factor series, it is reasonable to run the forecasting algorithm with only a few last observations. In simulation and application below, we employ the forecasting algorithm with the last 5 observations.
In this section, we briefly discuss what to expect from the forecasting method when the forecasting horizon $h$ becomes longer. Recall that the latent factor series is stationary and follows a stable VAR model. The latent process $\{Z_t\}_{t\in\mathbb{Z}}$ is also stationary and its long-term prediction converges to its mean, which is a zero vector. From (ref) and (ref), the predicted particles are therefore expected to converge eventually to zero vectors as well. Thus, for longer horizon $h$, we expect
As in Remark (ref), consider $V(x) = V(x_i) = \boldsymbol{1}_{\{x_i=n_i\}}$, $i=1,\ldots,d$. For $\widehat{Z}_{t+h|t}=z=0$, (ref) becomes
where $(\widehat{R}_{t+h|t})_{ii}$ is the $i$th diagonal entries of covariance matrix $\widehat{R}_{t+h|t}$. For longer horizon $h$, $(R_{t+h|t})_{ii}\approx \textrm{Var}(Z_{i,t})=1$. We thus expect from (ref) and (ref) that for longer horizon $h$,
where $Z \sim \mathcal{N}(0,1)$ and $X_i$ has a CDF $F_i$ with the parameter $\theta_i$. Hence, for longer horizon $h$, our forecasting method can be thought as choosing $n_i$ among possible values that maximizes the most likely count value according to the distribution $F_i$ with the parameter $\theta_i$. Two observations are worth making in this regard.
First, in some instances, e.g. Bernoulli and categorical distributions, $\theta_i$ represents proportion of count values and is estimated as the corresponding sample proportions. In these cases, for long horizon $h$, we therefore expect our forecast to yield the most likely observed count (modulo the issue of ties). This is the case with the application considered in Section (ref). On the other hand, for many other distributions, this observation may not necessarily hold. For example, for the Poisson distribution, the parameter is taken as the sample mean and the most likely count according to this Poisson distribution does not need to be the most likely observed value, let alone be in the sample.
Second, the relation (ref) might be confusing from the following point of view. As argued above, for longer horizon $h$, we expect $\widehat{Z}_{t+h|t}\approx0$. In fact, we see this clearly in the application of Section (ref). The relation (ref) might be read as saying that 0 belongs to the bin $A_{i,n_i}$ with the highest standard normal probability. This is, however, not necessarily the case. It will be the case when $n_i$ is the median of the distribution $F_i$ (i.e. $F_{i}(n_i-1)<1/2$ and $F_{i}(n_i)\geq1/2$), a quite likely scenario in practice especially for “bell-shaped” distribution with the most likely value $n_i$ being at the center, but the statement will not hold in general.
The goodness of fit of a time series model for univariate count data can be assessed through a probability integral transformation (PIT) histogram. See czado2009predictive,kolassa2016evaluating,jia2023latent. PIT involves computing predictive distributions (see also below). We are not aware of PIT extensions to the multivariate setting, though natural possibilities can certainly be proposed. As with forecasting discussed in Sections (ref)--(ref), the bigger issue for our model and potentially high-dimensional setting is that obtaining a multivariate predictive distribution would be computationally infeasible. As in those sections, we shall compromise by computing marginal predictive distributions but conditioning on the past of all variables. More specifically, we proceed as described next.
Suppose that our model parameters $\theta,\Lambda,\Psi_{1},\ldots,\Psi_{p},\Sigma_{\varepsilon}$ and $\Sigma_{\eta}$ are estimated as in Sections (ref), (ref), and (ref). Then, the marginal predictive distribution of $i$th variable is defined as
It can be approximated similarly to the SMC approximation (ref) by
where $w_{t}^{(k)}=w_{t-1}^{(k)}w_{t}(\widehat{Z}_{t|t-1}^{(k)})$ as in (ref) and $\Omega_{N,t}=\sum_{k=1}^Nw_{t}^{(k)}$, and by using (ref),
where $\widehat{Z}_{t|t-1}^{(k)}$ is computed as in Section (ref) using $x_1,\ldots,x_{t-1}$, and $\widehat{R}_{t|t-1}$ is derived by Kalman recursions in Appendix (ref).
The (nonrandom) sample PIT for the $i$th component is defined as
where
In practice, the sample PIT of the $i$th component is a histogram with the height $\bar{F}_{i}(m/M) - \bar{F}_{i}((m-1)/M)$, $m=1,\ldots,M$, where $M$ is the number of bins (often $M=10$). If the model fits the data well, the sample PIT histogram should be close to a uniform distribution.
In this section, we assess the performance of estimation and forecasting procedures introduced in Sections (ref) and (ref). We focus on Bernoulli, categorical, Poisson, and negative binomial marginal distributions with several parameter values. Our goal is to check how the performance of our proposed methods compares to that of the traditional PCA methods in the factor modeling literature lam2011estimation,doz2011two,doz2012quasi.
The simulation settings are as follows. First, we fix the order $p=1$ for factor series $\{Y_t\}_{t\in\mathbb{Z}}$ in (ref) and take the number of factor series as $r=2$ or $r=5$. These values are assumed to be given in this section. Then, by following the standard DFM literature doz2011two, we generate the latent Gaussian DFM (ref)--(ref) with parameters:
These model parameters do not ensure that $\mathbb{E}Z_{i,t}^2=1$ as needed for (ref), so further standardization is applied. We take $d=15,30,60,90$ and $T=100,200$. Finally, to define the marginal count distributions, let $I_1,I_2,I_3$ be the sets of indices that partition $\{1,2,\ldots,d\}$ such that $I_1=\{1,\ldots,d/3\}$, $I_{2}=\{d/3+1,\ldots,2d/3\}$, and $I_3=\{2d/3+1,\ldots,d\}$. Then,
Monte Carlo simulations are based on 100 replications for each setting. For the 100 replications, we report the mean and the standard deviations of $\ell_2$ losses $\|\widehat{a}-a\|_2/\sqrt{b}$ of the estimators, where $\widehat{a}$ is the (vectorized) estimator of the parameter $a$ and $b$ is a scalar; $b=d$ for $\theta$ in (ref) and $\Lambda,\Sigma_{\varepsilon}$ in (ref), and $b=r$ for $\Psi_1,\Sigma_{\eta}$ in (ref).
Table (ref) summarizes the estimation results for several $r,d,T$ values, and marginal distributions. The results align with what is expected from the standard factor models. For example, the means of the losses decrease with increasing dimension $d$ and sample size $T$. When the number of factors $r$ increases, both the averages of the losses of the estimators $(\widehat{\Lambda},\widehat{\Sigma}_{\varepsilon})$ for the factor model (ref) and those of $(\widehat{\Psi},\widehat{\Sigma}_{\eta})$ for the factor series (ref) increase. Another interesting point is that the means of the losses for the factor series are larger than those for the factor model. On the observation level, the means of the losses of estimators $\widehat{\theta}$ of marginal distributions are relatively small, regardless of $d,T$, and even $r$. Furthermore, the magnitudes of the means of the losses are not substantially different across different marginal count distributions.
We investigate the performance of the selection methods of the number of factors $r$, suggested in Section (ref). The same model parameters as in Section (ref) with $d=30,50$, and 100 are used with fixed $p=1$. We denote the scree plot method of finding the “knee” described in (ref) by ED. The IC methods as combinations of (ref) with the three different penalty functions (ref) are denoted by IC1--IC3, respectively. For the BCV-based approach, we employ two different estimation procedures. Our principal component-based estimation in Section (ref) is denoted by PC. On the other hand, Fac refers to MINRES estimation by harman1966factor, which is briefly described in Section (ref). For each simulation setting, 100 replications are performed.
Figure (ref) depicts the frequencies of estimated $r$ in 100 replications. From the figure, PC method outperforms all baselines (IC1--3, ED) in selecting true $r$, and performs similarly to Fac. Interestingly, the traditional information criterion-based or scree plot-based approaches fail for this model. The quality of estimation by the BCV-based approaches follows the pattern in the estimation of model parameters. That is, as the dimension $d$ and sample length $T$ increase, so do the percentages of correctly estimated $r$. Another observation is that the larger number of factor series deteriorates the performance of cross-validations. This can be seen for both cross-validation schemes, which is also an expected phenomenon with the standard factor models. In terms of marginal distributions, cross-validation schemes work best for negative binomial, followed by Poisson, categorical, and Bernoulli marginal distributions. This indicates that cross-validation schemes improve when marginal distributions tend to take a larger number of values so that they are more akin to continuous distributions.
In this section, we assess forecasting performance in the simulation settings considered in Section (ref), with sample lengths for the data generation extended by 12 observations. We then hold out the last 12 observations so that the targeted forecasting horizons are $H=1,2,3,6$, and $12$. We assume that the model parameters are given, allowing for the estimation of covariances $\widehat{Q}_{t|t-1}$ and $\widehat{Q}_{t|t}$, the Kalman gain $K_t$, and the importance weights ${\widetilde{w}_{T}^{(k)}}$ to be conducted within the last 5 observations from the sample, as discussed in Section (ref). With these estimators, we generate $h$-step-ahead predictions $V(\widehat{X}_{T+h|T})$ through (ref). The remainder of the prediction follows the Kalman recursions as discussed in Appendix (ref).
Assuming the data is generated from our model, we would like to examine our forecasting scheme, especially compared to alternative approaches. We refrain from comparing forecasting based on other multivariate discrete-valued time series models mentioned in Section (ref). The MATLAB codes for Bayesian DFMs by cui2014generalized and nonstationary DFMs by wang2018modeling are publicly available. However, those models were not considered for forecasting, necessitating additional methodological and implementation considerations. Furthermore, the assumed dimensions in those models are relatively low, typically five or fewer, compared to 15 or more in our study. As a result, we consider three naive baselines instead. First, as the most naive forecasting method, we use the last observation. We refer to this method as Last. Second, we consider each dimension $i=1,\ldots,d$ separately and define the predictions as the likeliest previous value for that dimension. For example, if the most frequent value $X_{1,t}$ for the first dimension $i=1$ within the observation window is 3 for the Poisson case, we take 3 as predictions for all 5 steps ahead. We refer to this approach as Marginal; see also Section (ref). Third, we also consider the case discussed at the end of Section (ref) where the forecast is taken as the discrete value associated with the bin containing the value $Z=0$. For example, if the $i$th variable follows Bernoulli distribution with $p_i=0.4$, the forecast is 0 since $0<\Phi^{-1}(1-0.4)$. On the other hand, if $p_i=0.7$, the forecast becomes 1 since $0>\Phi^{-1}(1-0.7)$. Forecasts for other distributions can be defined similarly by using (ref). We refer to this approach as Null.
The Monte Carlo simulations are based on 100 replications. We consider two types of forecasting measures: one for the latent series and one for the observed series. For the latent Gaussian series, we compute means and standard deviations of the root mean square error (RMSE) of $H$-step ahead forecast error averaged over $N$ particles defined by
For the observed series, we consider an out-of-sample loss, similar to (ref) by replacing $N=100$ and $\widehat{Y}_{T+H|T}^{(k)},Y_{T+H}$ or $\widehat{Z}_{T+H|T}^{(k)},Z_{T+H}$ with $N=1$ and $\widehat{X}_{T+H|T},X_{T+H}$, respectively. In addition, applying this performance measure to $X_t$ may be less adequate due to the discrete nature of values. So, forecasting accuracy for the observations using the proposed methods, including the two baselines, is also measured by the accuracy of $H$-step ahead forecasting,
Table (ref) reports the results of means and standard deviations of the RMSE of $H$-step ahead forecast error from the two latent processes $Z_t$ and $Y_t$, defined in (ref). The increasing trends for the RMSEs of those two forecasting errors along the forecasting horizon $H$ are what one would expect from forecasting results for typical time series models. The RMSEs of the forecasting errors of the factor series (ref) are comparably larger than those of the factor models (ref). Furthermore, note that the RMSEs of the forecasting errors depend on both the dimension $d$ and the number of factors $r$ but in different ways: while the increase of $d$ leads to smaller errors of both the factor series and the factor models, the increase of $r$ causes the increase and decrease of errors in the factor models and factor series, respectively.
Table (ref) reports the results from the observations $X_t$, in terms of the mean of RMSEs of forecasting errors and accuracy of forecasting including three benchmarks. Similarly to the two latent processes, the RMSEs of the forecasting errors increase as the forecasting horizon increases. Furthermore, the RMSEs of the forecasting errors for the observed series tend to be larger than those for the latent processes. It is interesting to note the larger forecasting errors for marginal distributions that take more values, such as negative binomial, followed by Poisson, categorical, and Bernoulli distributions. This was not the case for the RMSEs of the forecasting errors of latent processes. Similar tendencies as for the RMSEs of the forecasting errors are observed with the accuracy measure. Finally, compared to Last, Marginal, and Null approaches, consistently better forecasting performance is noted for the proposed method (the larger accuracy is associated with better predictors).
To demonstrate the utility of the proposed model, we consider individual-level time series consisting of daily self-reported measures of personality, collected by borkenau1998big. That study was designed to explore items describing personal emotions that represent the “Big Five” personality factors. The structure of the data is as follows. Time series for 30 emotion-related items were collected, where groups of 6 items are known to correspond to one of the five personality factors. These data were collected for 22 students over 90 days.
All 30 items are believed to correspond to at least one of the “Big Five” personality factors. These factors are known to follow a structured categorization, denoted as categories 1 through 5 below. Each evening, participants in the study were instructed to evaluate their daily behavior for each item on a scale from 0 to 6, with 6 indicating the strongest endorsement of the emotion that day. For illustrative purposes, we focus on the data from one student out of the 22 available. This choice of the student was largely motivated by the following consideration: for many students, responses to some items exhibited very little variability (i.e., remained mostly constant), making them less interesting from a time series modeling perspective. To mitigate the influence of extreme observations 0 and 6, we merged them with 1 and 5, respectively, resulting in a new scale ranging from 1 to 5. Practitioners intending to use our model (or any other) should first conduct similarly basic exploratory data analysis.
The data set and its time series plots are presented in Figure (ref). Among 90 consecutive observations, we use the first 85 observations for estimation and reserve the last 5 as a holdout sample to evaluate our forecasts. The corresponding 30 items are grouped by the identified categories C1 through C5, each containing 6 items. The time series of individual items exhibit substantial variability. In addition, the time plots reveal that the dynamics of the 6 items within each category share common patterns. This suggests the plausible existence of a latent factor structure underlying these dynamics.
The following remarks provide further evidence for the latent factor structure. Several estimated correlation matrices are depicted in Figure (ref). The top panel presents the sample autocorrelation matrices of the observed series $\{X_t\}_{t=1,\ldots,85}$, while the bottom panel shows the estimated autocorrelation matrices of the latent series $\{Z_t\}_{t=1,\ldots,85}$. In both panels, the left plot corresponds to lag 0, and the right plot corresponds to lag 1. One can see that both autocorrelation matrices for the same lag order are nearly indistinguishable. At lag 0, the plots in both panels show clear block patterns characteristic of the factor structure. Furthermore, the factor structure appears to be preserved through temporal dependence, as suggested by the plots of the sample ACFs at lag 1, though it is less discernible.
We work with the model assuming $r=5$ and $p=1$ for the sake of illustration. We also assume categorical distributions on $\{1,2,3,4,5\}$ as marginal distributions. For the loadings matrices, we use the second identifiability condition (ref). The left panel of Figure (ref) presents the estimated loadings matrix, while the right panel depicts the transition matrix of the factor series, estimated using the method described in Section (ref). This pattern is consistent with what we observe in Figure (ref). Note that the transition matrix has relatively large values not only on the diagonal but also for some off-diagonal entries. The large off-diagonal values indicate that there are cross-correlations across the factor series.
For the model diagnostics, Figure (ref) shows the PIT histograms for 30 individual items, following the procedure in Section (ref). With the exception of a few items, the heights of the PIT histograms across 10 bins tend to be uniformly equal, aligning with the 0.1 relative frequencies indicated by the black dashed lines.
With the estimated model, we forecast the next 5 steps as discussed in Remark (ref) and compare the values with the true ones held out of the sample. For comparison, we consider two simple forecasting approaches, Last and Marginal, as described in Section (ref).
Figure (ref) presents the generated particles $\{\widehat{Z}_t^{(k)}\}_{t=1,\ldots,T}$ for 30 items within the last 5 observations $t=81,\ldots,85$ and the predicted values of latent process $\{\widehat{Z}_{T+h|T}^{(k)}\}_{h=1,\ldots,5}$ for the next 5 forecasting steps, respectively. The items are grouped into five underlying categories. The three or four parallel lines in Figure (ref) represent the thresholds, with each pair of lines forming a bin. To distinguish the values, we use different line types for each threshold. As explained above, the particles are generated at each time point to belong to a specific bin that corresponds to the discrete observation for that dimension. This explains why the particles stay within bins at all time points in the left panel. For some of the items, a thresholding line extends beyond the given vertical scale due to the way the bins are defined. Since each bin is estimated through the observations, some values do not appear if they are not realized in the observation period. These values are also excluded from the candidate forecasting values. On the other hand, since no further observations are assumed to be given after the observation period, the particles are generated by forecasting the latent process. As a result, the particles in the right panel do not need to stay in the same bin. Note that all of the particles seem to converge to zero for the longer forecasting horizon. As discussed in Section (ref), this behavior is expected when a stable VAR is used for forecasting. In addition, the particles are rather close to zero even for the first few horizons. This is a consequence of the interplay between the levels of signal (factors) and noise (errors) in the estimated model. The forecasted value naturally approaches 0 when forecasting noise, which in turn downweights the factor forecast as the noise magnitude increases, given that the variance of our latent process is 1 in each dimension.
Finally, Figure (ref) shows the plots of the absolute differences between forecasts and the true values for each item. The items are arranged according to their categories, under the expectation that each factor primarily influences its corresponding category. Overall, the proposed forecasting approach slightly outperforms the reference methods, generally exhibiting small absolute differences. As discussed in Section (ref), the forecasting performance becomes identical to the Marginal for the longer horizon. But for shorter horizons, our approach performs better for 4 items, whereas the Marginal method performs better for 2 items. While this suggests a potential advantage, it is not sufficient to draw overarching conclusions based on this evidence.
In this work, we considered a multivariate discrete-valued times series model, wherein component count series are obtained by binning the continuous values of latent Gaussian dynamic factor processes. We introduced an estimation method based on second-order properties of the count and latent processes, and PCA. We also suggested additional model selection approaches for determining the number of factor series and their lag orders through cross-validation. We provided the theoretical guarantees of the estimators by applying available concentration results for the considered model with general latent Gaussian processes. Facilitated by the state-space formulation of our model, we employed a sequential Monte Carlo method with resampling, for forecasting. Our estimation and forecasting methods were examined on simulated data and an empirical example. The R code used for the illustrations in Section (ref), the simulation study of Sections (ref), and the data analysis of Section (ref) are available on GitHub at \href{https://github.com/yk748/latent_Gaussian_TS}{https://github.com/yk748/latent_Gaussian_TS}.
While our study advances a framework for latent Gaussian time series modeling of categorical observations collected over time, important questions remain. The question of how to employ time-varying covariates is discussed in Supplemental Material but left for future work. There are also potential improvements to make in terms of accuracy and computing time of our forecasting methods. Instead of employing standard particle filtering strategies, one could try other variants of sequential Monte Carlo sampling, for example, ensemble Kalman filtering in a high-dimensional regime katzfuss2020ensemble.
MD was supported partially by the FAU Emerging Talents Initiative. VP was supported partially by the grants NSF DMS-1712966, DMS-2113662, and DMS-2134107. The authors also thank two anonymous Reviewers for their many helpful comments.