Modelling Volatility of Spatio-temporal Integer-valued Data
with Network Structure and Asymmetry

Integer-valued time series can be observed in a wide range of scientific fields, such as the yearly trading volume of houses on real estate market (De Wit et al., , 2013), number of transactions of stocks (Jones et al., , 1994), daily mortality from Covid-19 (Pham, , 2020) and daily new cases of Covid-19 (Jiang et al., , 2022). A first idea to model integer-valued time series is using a simple first-order autoregressive model (AR)

$$X_{t}=\alpha X_{t-1}+\varepsilon_{t},$$

(1.1)

where $0\leq\alpha<1$ is a parameter. However in (1.1) $X_{t}$ is not necessarily an integer given integer-valued $X_{t-1}$ and $\varepsilon_{t}$, due to the multiplication structure $\alpha X_{t-1}$. To circumvent such problem, McKenzie, (1985) and Al-Osh and Alzaid, (1987) proposed an integer-valued counterpart of AR model (INAR). They replaced the ordinary multiplication $\alpha X_{t-1}$ by the binomial thinning operation $\alpha\circ X_{t-1}$, where $\alpha\circ X|X\sim Bin(X,\alpha)$. This was ground-breaking and led to various extensions of thinning-based linear models including integer-valued moving average model (INMA) (Al-Osh and Alzaid, , 1988) and INARMA model (McKenzie, , 1988). An alternative approach to the multiplication problem is to consider the regression of the conditional mean $\lambda_{t}:=\mathbb{E}(X_{t}|\mathcal{H}_{t-1})$, where $H_{t-1}$ is the $\sigma$-algebra generated by historical information up to $t-1$. Based on this idea, integer-valued GARCH-type models (INGARCH) were proposed by Heinen, (2003), Ferland et al., (2006) and Fokianos et al., (2009) with conditional Poisson distribution of $X_{t}$, i.e.

		$\displaystyle X_{t}|\mathcal{H}_{t-1}\sim Poisson(\lambda_{t}),$		(1.2)
		$\displaystyle\lambda_{t}=\omega+\sum_{i=1}^{p}\alpha_{i}X_{t-i}+\sum_{j=1}^{q}\beta_{j}\lambda_{t-j},$

where parameters satisfy $\omega>0,\alpha_{i}\geq 0,i=1,\cdots,p,\beta_{j}\geq 0,j=1,\cdots,q$. Other variations of INGARCH models with different specifications of conditional distribution include negative binomial INGARCH (Zhu, , 2010; Xu et al., , 2012) and generalized Poisson INGARCH (Zhu, , 2012).

Integer-valued models mentioned above are all limited to one-dimensional time series. The development of multivariate integer-valued GARCH-type models is still at its early stage. For example, the bivariate INGARCH models (Lee et al., , 2018; Cui and Zhu, , 2018; Cui et al., , 2020), multivariate INGARCH models (Fokianos et al., , 2020; Lee et al., , 2023) and some counterparts of the continuous-valued network GARCH model (Zhou et al., , 2020) are on fixed low-dimensional time series of counts. For high-dimensional integer-valued time series with increasing dimension, there is the Poisson network autoregressive model (PNAR) by Armillotta and Fokianos, (2024). The PNAR allows for integer-valued time series with increasing network dimension. However, it adopted a ARCH-type structure without considering the autoregressive term on the conditional mean/variance to describe persistence in volatility, and it can not capture asymmetric characteristics of volatility. Tao et al., (2024) proposed a grouped PNAR model which has a GARCH structure, but its network dimension is fixed and not suitable to spatio-temporal integer-valued data. In this paper, we propose a Poisson threshold network GARCH model (PTNGARCH) and discuss its asymptotic inferences. The specific contributions are as follows.

1. 1.

   A threshold structure is designed so that it is able to capture asymmetric property of high-dimensional volatility for integer-valued data. The threshold effect can also be tested under such a framework.

2. 2.

   Our PTNGARCH includes an autoregressive term on the conditional mean/variance so that it provides a parsimonious description of dynamic volatility persistence of high-dimensional integer-valued time series.

3. 3.

   Asymptotic theory, when both sample size ($T$) and network dimension ($N$) are large, of maximum likelihood estimation for the proposed model is established by the limit theorems for arrays of weakly dependent random fields in Pan and Pan, (2024). This enables the proposed model to be suitable for fitting dynamic spatio-temporal integer-valued data (or the so-called ultra-high dimensional integer-valued time series).

The remainder of this paper is arranged as follows. The PTNGARCH model is introduced and its stationarity over time is discussed under fixed network dimension in section 2. Section 3 presents the MLE for parameters, including the threshold, and the consistency and asymptotic normality for estimates of coefficients under large sample size ($T$) and large network dimension ($N$). Section 4 gives a Wald test which can be applied to testing of the existence of threshold effect (i.e. asymmetric property), GARCH effect and network effect. Section 5 conducts simulation studies to verify the asymptotic properties of the MLE, and applies the proposed model to daily numbers of car accidents that occurred in 41 neighbourhoods in New York City, with interpretation of the results of analysis. All proofs of theoretical results are postponed to the appendix.

Consider an non-directed and weightless network with $N$ nodes. Define adjacency matrix $A=(a_{ij})_{1\leq i,j\leq N}$, where $a_{ij}=1$ if there is a connection between node $i$ and $j$, otherwise $a_{ij}=0$. In addition, self-connection is not allowed for any node $i$ by letting $a_{ii}=0$. As an interpretation of the network structure, $A$ is symmetric since $a_{ij}=a_{ji}$, hence for any node $i$, the out-degree $d_{i}^{(out)}=\sum_{j=1}^{N}a_{ij}$ is equal to the in-degree $d_{i}^{(in)}=\sum_{j=1}^{N}a_{ji}$, and we use $d_{i}$ to denote both for convenience. To embed a network into statistical models, it is often convenient to use the row normalized adjacency matrix $W$ with its $(i,j)$ element $w_{ij}=\frac{a_{ij}}{d_{i}}$.

For any node $i\in\{1,\cdots.N\}$ in this network, let $y_{it}$ be an non-negative integer-valued observation at time $t$, and $\mathcal{H}_{t-1}$ denotes the $\sigma$-algebra consisting of all available information up to $t-1$. In our Poisson threshold network GARCH model, for each $i=1,2,...,N$ and $t\in\mathbb{Z}$, $y_{it}$ is assumed to follow a conditional (on $\mathcal{H}_{t-1}$) Poisson distribution with $(i,t)$-varying variance (mean) $\lambda_{it}$. A PTNGARCH(1,1) model has the following form:

		$\displaystyle y_{it}|\mathcal{H}_{t-1}\sim Poisson(\lambda_{it}),$		(2.1)
		$\displaystyle\lambda_{it}=\omega+\left(\alpha^{(1)}1_{\{y_{i,t-1}\geq r\}}+\alpha^{(2)}1_{\{y_{i,t-1}<r\}}\right)y_{i,t-1}+\xi\sum_{j=1}^{N}w_{ij}y_{j,t-1}+\beta\lambda_{i,t-1},$
		$\displaystyle i=1,2,\cdots,N.$

The threshold parameter $r$ is a positive integer, and $1_{\{\cdot\}}$ denotes an indicator function. To assure the positiveness of conditional variance, we need to assume positiveness of the base parameter $\omega$, and non-negativeness of all the coefficients $\alpha^{(1)}$, $\alpha^{(2)}$, $\xi$, $\beta$.

Let $\{N_{it}:i=1,2,...,N,t\in\mathbb{Z}\}$ be independent Poisson processes with unit intensities. Depending on $\lambda_{it}$, $y_{it}$ can be interpreted as a Poisson distributed random variable $N_{it}(\lambda_{it})$, which is the number of occurrences during the time interval $(0,\lambda_{it}]$, i.e. $\mathbb{P}(y_{it}=n|\lambda_{it}=\lambda)=\frac{\lambda^{n}}{n!}e^{-\lambda}.$ We could rewrite (2.1) in a vectorized form as follows:

$$\begin{cases}&\mathbb{Y}_{t}=(N_{1t}(\lambda_{1t}),N_{2t}(\lambda_{2t}),...,N_{Nt}(\lambda_{Nt}))^{\prime},\\ &\Lambda_{t}=\omega\mathbf{1}_{N}+A(\mathbb{Y}_{t-1})\mathbb{Y}_{t-1}+\beta\Lambda_{t-1},\end{cases}$$

(2.2)

where

	$\displaystyle\Lambda_{t}=(\lambda_{1t},\lambda_{2t},...,\lambda_{Nt})^{\prime}\in\mathbb{R}^{N},$
	$\displaystyle\mathbf{1}_{N}=(1,1,...,1)^{\prime}\in\mathbb{R}^{N},$
	$\displaystyle A(\mathbb{Y}_{t-1})=\alpha^{(1)}S(\mathbb{Y}_{t-1})+\alpha^{(2)}(I_{N}-S(\mathbb{Y}_{t-1}))+\xi W,$
	$\displaystyle S(\mathbb{Y}_{t-1})=\operatorname*{diag}\left\{1_{\{y_{1,t-1}\geq r\}},1_{\{y_{2,t-1}\geq r\}},...,1_{\{y_{N,t-1}\geq r\}}\right\}.$

Note that $\mathbb{Y}_{t}\in\mathbb{N}^{N}$ with dimension $N$ and $\mathbb{N}=\{0,1,2,\cdots\}$. For a fixed dimension $N$, the stationarity condition of time series $\mathbb{Y}_{t}$ can be obtained. The proof of the following theorem is given in the appendix.

Denote parameter vector of model (2.1) by $\nu=(\theta^{\prime},r)^{\prime}$ with $\theta=(\omega,\alpha^{(1)},\alpha^{(2)},\xi,\beta)^{\prime}$. Let $\Theta=\{\theta:\omega>0,\alpha^{(1)}\geq 0,\alpha^{(2)}\geq 0,\xi\geq 0,\beta\geq 0\}$ and $\mathbb{Z}_{+}=\{0,1,2,\cdots\}$. For a given threshold $r$, the reasonable parameter space for $\theta$ is a sufficiently large compact subset of $\Theta$, such as $\Theta_{\delta}=\{\theta:\omega\geq\delta,\alpha^{(1)}\geq 0,\alpha^{(2)}\geq 0,\xi\geq 0,\beta\geq 0\}$ with sufficient small positive real number $\delta$.

Let $D_{NT}=\{(i,t):i=1,\cdots,N;t=1,\cdots,T\}$ be the index set. Suppose that the samples $\{y_{it}:(i,t)\in D_{NT},NT\geq 1\}$ are generated by model (2.1) with respect to true parameters $\nu_{0}=(\omega_{0},\alpha^{(1)}_{0},\alpha^{(2)}_{0},\xi_{0},\beta_{0},r_{0})^{\prime}$.

The log-likelihood function (ignoring constants) is

$$\left\{\begin{array}[]{ll}&L_{NT}(\nu)=\frac{1}{NT}\sum_{(i,t)\in D_{NT}}l_{it}(\nu),\\ &l_{it}(\nu)=y_{it}\log\lambda_{it}(\nu)-\lambda_{it}(\nu)\end{array}\right.$$

(3.1)

where $\lambda_{it}(\nu)$ is generated from model (2.1) as

	$\displaystyle\lambda_{it}(\nu)=$	$\displaystyle\omega+\alpha^{(1)}1_{\{y_{i,t-1}\geq r\}}y_{i,t-1}+\alpha^{(2)}1_{\{y_{i,t-1}<r\}}y_{i,t-1}$		(3.2)
		$\displaystyle+\xi\sum_{j=1}^{N}w_{ij}y_{j,t-1}+\beta\lambda_{i,t-1}(\nu).$

In practice, (3.1) can not be evaluated without knowing the true values of $\lambda_{i0}$ for $i=1,2,...,N$. We need to approximate (3.1) by (3.3) below, using specified initial values $\lambda_{i0}=\tilde{\lambda}_{i0},i=1,2,...,N$:

$$\left\{\begin{array}[]{ll}&\tilde{L}_{NT}(\nu)=\frac{1}{NT}\sum_{(i,t)\in D_{NT}}\tilde{l}_{it}(\nu),\\ &\tilde{l}_{it}(\nu)=y_{it}\log\tilde{\lambda}_{it}(\nu)-\tilde{\lambda}_{it}(\nu).\end{array}\right.$$

(3.3)

Therefore, the maximum likelihood estimate (MLE) of parameter $\nu$ is defined as

$$\hat{\nu}_{NT}=\operatorname*{argmax}_{\nu\in\Theta\times\mathbb{Z}_{+}}\tilde{L}_{NT}(\nu).$$

(3.4)

However, the solution that maximizes the target function $\tilde{L}_{NT}(\nu)$ can not be directly obtained by solving $\frac{\partial\tilde{L}_{NT}(\nu)}{\partial\nu}=0,$ since $r\in\mathbb{Z}_{+}$ is discrete-valued, therefore the partial derivative of $\tilde{L}_{NT}(\nu)$ w.r.t. $r$ is invalid. According to Wang et al., (2014), such optimization problem with integer-valued parameter $r$ could be broken up into two steps as follows:

1. 1.

   Find

   $$\hat{\theta}_{NT}^{(r)}=\operatorname*{argmax}_{\theta\in\Theta}\tilde{L}_{NT}(\theta,r)$$

   for each $r$ in a predetermined range $[\underline{r},\bar{r}]\subset\mathbb{Z}_{+}$;

2. 2.

   Find

   $$\hat{r}_{NT}=\operatorname*{argmax}_{r\in[\underline{r},\bar{r}]}\tilde{L}_{NT}(\hat{\theta}_{NT}^{(r)},r).$$

Then $\hat{\nu}_{NT}=\left(\hat{\theta}_{NT}^{(\hat{r}_{NT})^{\prime}},\hat{r}_{NT}\right)^{\prime}$ would be the maximizer of $\tilde{L}_{NT}(\nu)$.

Assumption 3.1 below is a regular condition on the parameter space. Assumptions 3.2 and 3.3 are necessary for obtaining $\eta$-weak dependence of $\{l_{it}(\nu):(i,t)\in D_{NT},NT\geq 1\}$ and their derivatives. For the definition of $\eta$-weak dependence for random fields, see Pan and Pan, (2024). Then the consistency of MLE in Theorem 2 can be proved based on the law of large numbers (LLN) of $\eta$-weakly dependent arrays of random fields in Pan and Pan, (2024).

Since $\hat{r}_{NT}$ is an integer-valued consistent estimate of $r_{0}$, $\hat{r}_{NT}$ will eventually be equal to $r_{0}$ when the sample size $NT$ becomes sufficiently large. Therefore, $\hat{\nu}_{NT}=\left(\hat{\theta}_{NT}^{(\hat{r}_{NT})^{\prime}},\hat{r}_{NT}\right)^{\prime}$ is asymptotically equal to $\left(\hat{\theta}_{NT}^{(r_{0})^{\prime}},r_{0}\right)^{\prime}$. In this way, the problem of investigating the asymptotic distribution of $\hat{\nu}_{NT}$ degenerates to investigating the asymptotic distribution of $\hat{\theta}_{NT}^{(r_{0})}$.

Based on Theorem 2 and Theorem 3, for sufficiently large sample region such that $\hat{r}_{NT}=r_{0}$, we are able to design a Wald test for the null hypothesis

$$H_{0}:\Gamma\theta_{0}=\eta,$$

(4.1)

where $\Gamma$ is an $s\times 5$ matrix with rank $s$ ( $1\leq s\leq 5$) and $\eta$ is an $s$-dimensional vector. For example, to test the existence of threshold effect, simply let $\Gamma=(0,1,-1,0,0)$ and $\eta=0$, and the null hypothesis (4.1) becomes

$$H_{0}:\alpha^{(1)}_{0}=\alpha^{(2)}_{0};$$

To see if the autoregressive term is necessary (i.e. GARCH effect), take $\Gamma=(0,0,0,0,1)$ and $\eta=0$, and the hypothesis becomes

$$H_{0}:\beta_{0}=0\quad v.s.\quad H_{1}:\beta_{0}>0;$$

To test the network effect, just let $\Gamma=(0,0,0,1,0)$ and $\eta=0$, and the question becomes testing the hypothesis

$$H_{0}:\xi_{0}=0\quad v.s.\quad H_{1}:\xi_{0}>0.$$

Corresponding to the asymptotic normality of $\hat{\theta}_{NT}^{(r_{0})}$ in Theorem 3, we define a Wald test statistic as follows

$$W_{NT}:=(\Gamma\hat{\theta}_{NT}^{(r_{0})}-\eta)^{\prime}\left\{\frac{1}{NT}\Gamma\widehat{\Sigma}_{NT}^{-1}\Gamma^{\prime}\right\}^{-1}(\Gamma\hat{\theta}_{NT}^{(r_{0})}-\eta),$$

(4.2)

where $\widehat{\Sigma}_{NT}$ is defined as (3.6). The following Theorem 4 shows that $W_{NT}$ has an asymptotic $\chi^{2}$-distribution with degree of freedom $s$.

We have proposed a model to describe spatio-temporal integer-valued data which are observed at nodes of a network and have asymmetric property. Although the proposed model is called GARCH-type for volatility (conditional variance), it also can be used to model the conditional mean, because we assume that the conditional distribution is Poisson distribution. Asymptotic inferences, parameter estimation and hypothesis testing of the proposed model, are discussed by applying limit theorems for nonstationary arrays of random fields. Simulations for different network structures and application to a real dataset show that our methodology is useful for modelling dynamic spatio-temporal integer-valued data.