Adjacent-category models for ordinal time series and their application to climate-dependent spruce budworm defoliation dynamics

We denote $Y_{t}$ as the response variable. In our example in Section 3, $Y_{t}$ corresponds to the defoliation level for each area classified on a scale of 0 to 3, corresponding to no defoliation (approx. $0\%$) or a low (approx. $1\%-35\%$), moderate ($36\%-70\%$) or high ($71\%-100\%$) fraction of the year’s foliage defoliated by SBW. Our model is given by

	$\displaystyle Y_{k,t}=\mathds{1}\left\{\sum_{k=0}^{j-1}\pi_{k,t}\leq U_{t}<\sum_{k=0}^{j}\pi_{k,t}\right\},\quad$	$\displaystyle 0\leq j\leq K,$		(2.1)
	$\displaystyle\log\frac{\pi_{j,t}}{\pi_{j-1,t}}=:\eta_{j,t}=\omega_{j}+\gamma^{\top}X_{t-1}+\alpha^{\top}\overline{Y}_{t-1}+\beta_{j}\eta_{j,t-1},$	$\displaystyle\quad 1\leq j\leq K,t\in\mathbb{Z},$		(2.2)

where $\omega_{j}\in\mathbb{R},\beta_{j}\in\mathbb{R}^{*},j=1,\ldots,K,\gamma\in\mathbb{R}^{P}\text{ and }\alpha\in\mathbb{R}^{K}$. The vector $\overline{Y}_{t}=(Y_{1,t},\ldots,Y_{K,t})\in\{0,1\}^{K}$, where for $1\leq j\leq K,\leavevmode\nobreak\ Y_{j,t}=1$ when $Y_{t}=j.$ The vector $(Y_{0,t},\ldots,Y_{K,t})$ is then the one-hot encoding of the categorical random variable $Y_{t}$ valued in $\{0,\ldots,K\}.$ It is well known that the parameters $\beta_{j},j=1,\ldots,K$ act as feedback parameters; $\beta_{j},j=1,\ldots,K$ allow us to extend the memory of Model (2.1)-(2.2) beyond lag 1 when the effect of lagged values decreases quickly when $j=1,\ldots,K,|\beta_{j}|<1.$

From Equation (2.2),

$$\eta_{t}=\omega+\Gamma X_{t-1}+A\overline{Y}_{t-1}+B\eta_{t-1},$$

where $\omega=(\omega_{j})_{1\leq j\leq K},\Gamma^{\top}=(\gamma|\cdots|\gamma),A^{\top}=(\alpha|\cdots|\alpha)$ and $B=\mathrm{diag}(\beta_{j},1\leq j\leq K)$. The vector of the parameters of the model is denoted by $\theta=(\omega^{\top},\mathrm{vec}(\Gamma)^{\top},\alpha^{\top},\beta_{1},\ldots,\beta_{K})^{\top}$, and we assume that $\Theta(\ni\theta)$ is a compact set. The true value of the parameter is $\theta_{0}$. For what follows, we write $\eta_{t}(\theta)$ (resp. $\eta_{j,t}(\theta),\pi_{j,t}(\theta),1\leq j\leq K)$ to make the latent process depending on the parameter $\theta.$

For $k=1,\ldots,K$,

$$\pi_{k,t}(\theta)=\frac{\mathrm{e}^{\sum_{j=1}^{k}\eta_{j,t}(\theta)}}{1+\sum_{i=1}^{K}\mathrm{e}^{\sum_{j=1}^{i}\eta_{j,t}(\theta)}},\quad\theta\in\Theta,$$

and the parameter $\theta$ can by estimated by

$$\hat{\theta}_{n}=\underset{\theta\in\theta}{\mathrm{argmin}}-\sum_{t=1}^{n}\left(\sum_{k=1}^{K}Y_{k,t}\left(\sum_{j=1}^{k}\eta_{j,t}(\theta)\right)-\log\left(1+\sum_{k=1}^{K}\mathrm{e}^{\sum_{j=1}^{k}\eta_{j,t}(\theta)}\right)\right)$$

(2.6)

with initial values of $\eta_{0}$ (resp. $\overline{Y}_{0}$). This estimator corresponds to the conditional maximum likelihood estimator of the Model (2.1)-(2.2). Let us set

$\displaystyle s_{t}(\theta)=-\sum_{k=1}^{K}\left(\sum_{j\geq k}Y_{j,t}-\frac{\sum_{j\geq k}\mathrm{e}^{\sum_{i=1}^{j}\eta_{i,t}(\theta)}}{1+\sum_{k=1}^{K}\mathrm{e}^{\sum_{j=1}^{k}\eta_{j,t}(\theta)}}\right)\nabla_{\theta}\eta_{k,t}(\theta)$

and

	$\displaystyle h_{t}(\theta)$	$\displaystyle=$	$\displaystyle\sum_{k=1}^{K}\left\{\frac{1}{\left(1+\sum_{k=1}^{K}\mathrm{e}^{\sum_{j=1}^{k}\eta_{j,t}(\theta)}\right)^{2}}\left[\left(\sum_{j\geq k}\mathrm{e}^{\sum_{i=1}^{j}\eta_{i,t}(\theta)}\right)\sum_{\ell=1}^{k-1}\left(1+\sum_{j=1}^{\ell-1}\mathrm{e}^{\sum_{i=1}^{j}\eta_{i,t}(\theta)}\right)\nabla_{\theta}\eta_{\ell,t}(\theta)+\right.\right.$
			$\displaystyle\left.\left.\left(1+\sum_{j=1}^{k-1}\mathrm{e}^{\sum_{i=1}^{j}\eta_{i,t}(\theta)}\right)\sum_{\ell=k}^{K}\left(\sum_{j\geq\ell}\mathrm{e}^{\sum_{i=1}^{j}\eta_{i,t}(\theta)}\right)\nabla_{\theta}\eta_{\ell,t}(\theta)\right]\right\}\nabla_{\theta}^{\top}\eta_{k,t}(\theta)$

with $\nabla_{\theta}^{\top}\eta_{k,t}(\theta)=\left(\frac{\partial{\eta_{k,t}}}{\partial{\theta_{i}}}\right)_{1\leq i\leq 3K+P}$ (see Appendix for the computations).

We do not prove the previous result. Interested readers will find some elements of the proof in Moysiadis and Fokianos, (2014), proof of Theorem 2. The same lines of proof can be found in Straumann et al., (2006) for the GARCH model or Debaly and Truquet, (2023) for a multivariate model for time series of mixed data. A condition for the identifiability assumption, $\eta_{0}(\theta)=\eta_{0}(\theta_{0})\leavevmode\nobreak\ \mathbb{P}_{\theta_{0}}-a.s\Rightarrow\theta=\theta_{0},$ based on Debaly and Truquet, (2023), is given in the Appendix. We asume the invertibility of the matrix $J$. Finally, in the estimation procedure for an observation-driven model, one often requires the stationary condition only for the true value of parameter $\theta_{0}$ and the condition for contraction for the map that defines the latent process (see Chap. 7 in Francq and Zakoian, (2019) or Straumann et al., (2006)). In our case, these coincide.

Here, we test the adequation of the Model (2.1)-(2.2) and the set of assumptions of Proposition 2. The null hypothesis is then that Model (2.1)-(2.2) holds with $\theta=\theta_{0}.$ One can look at the residuals:

$$e_{t}(\theta)=\left(\sum_{j\geq k}Y_{j,t}-\frac{\sum_{j\geq k}\mathrm{e}^{\sum_{i=1}^{j}\eta_{i,t}(\theta)}}{1+\sum_{k=1}^{K}\mathrm{e}^{\sum_{j=1}^{k}\eta_{j,t}(\theta)}}\right)_{1\leq k\leq K},t\in\mathbb{Z},\theta\in\Theta$$

(2.7)

and the autocorrelations:

$$\rho_{h}(\theta)=\frac{1}{n}\sum_{t=1}^{n}e_{t}(\theta)\odot e_{t-h}(\theta),h=1,\ldots,n,\theta\in\Theta.$$

(2.8)

The vector of $q$ successive autocorrelations is $\rho_{1:q}(\theta)=\mathrm{vec}((\rho_{1}(\theta)|\cdots|\rho_{q}(\theta))^{\top})$, and its evaluation at $\hat{\theta}_{n}$ is $\rho_{1:q}(\hat{\theta}_{n})=\mathrm{vec}((\rho_{1}(\hat{\theta}_{n})|\cdots|\rho_{q}(\hat{\theta}_{n}))^{\top}).$ The vector $\rho_{1:q}(\theta)$ evaluated at $\theta_{0}$ is denoted by $\rho_{1:q}.$ Note that $\rho_{1:q}=n^{-1}\sum_{t=1}^{n}m_{t}$ , and under the null hypothesis $\mathbb{E}(m_{t}|\mathcal{F}_{t-1})=0,$ i.e., $(m_{t})_{t\in\mathbb{Z}}$, is a martingale difference sequence. Obviously, there are many choices for residuals that make $(e_{t}(\theta_{0})\odot e_{t-h}(\theta_{0}))_{t\in\mathbb{Z}}$ a martingale difference sequence, such as

$$\left(Y_{k,t}-\frac{\mathrm{e}^{\sum_{j=1}^{k}\eta_{j,t}(\theta)}}{1+\sum_{k=1}^{K}\mathrm{e}^{\sum_{j=1}^{k}\eta_{j,t}(\theta)}}\right)_{1\leq k\leq K},t\in\mathbb{Z},\theta\in\Theta.$$

(2.9)

The residuals of interest here is relevant because

$$s_{t}(\theta)=-\sum_{k=1}^{K}e_{k,t}(\theta)\nabla_{\theta}\eta_{k,t}(\theta)$$

with $e_{k,t}(\theta)$ as the $k$-th element of $e_{t}(\theta)$. This allows us to derive the asymptotic distribution of $\sqrt{n}\rho_{1:q}.$

The Jacobian matrix $\xi_{t}$ of $e_{t}$ with respect to $\theta$ is

	$\displaystyle\xi_{t}(\theta)(k,\cdot)^{\top}$	$\displaystyle=$	$\displaystyle-\frac{1}{\left(1+\sum_{k=1}^{K}\mathrm{e}^{\sum_{j=1}^{k}\eta_{j,t}(\theta)}\right)^{2}}\left[\left(\sum_{j\geq k}\mathrm{e}^{\sum_{i=1}^{j}\eta_{i,t}(\theta)}\right)\sum_{\ell=1}^{k-1}\left(1+\sum_{j=1}^{\ell-1}\mathrm{e}^{\sum_{i=1}^{j}\eta_{i,t}(\theta)}\right)\nabla_{\theta}\eta_{\ell,t}(\theta)+\right.$
			$\displaystyle\left.\left(1+\sum_{j=1}^{k-1}\mathrm{e}^{\sum_{i=1}^{j}\eta_{i,t}(\theta)}\right)\sum_{\ell=k}^{K}\left(\sum_{j\geq\ell}\mathrm{e}^{\sum_{i=1}^{j}\eta_{i,t}(\theta)}\right)\nabla_{\theta}\eta_{\ell,t}(\theta)\right].$

For the sake of readability, we write $\hat{\rho}_{1:q}$ (resp. $e_{t}$ and $\xi_{t}$) for $\rho_{1:q}(\hat{\theta}_{n})$ (resp. $e_{t}(\theta_{0})$ and $\xi_{t}(\theta_{0})$).

Let us set $C_{q}=(c_{1}|\cdots|c_{q})^{\top}$ with $c_{k}=(c_{k,1}|\cdots|c_{k,q})$. For $h=1,\ldots,q,c_{k,h}^{\top}=\mathbb{E}e_{k,-h}\xi_{0}(k,\cdot),k=1,\ldots,K$, and $\hat{C}_{K}$ is its empirical counterpart. Let us set

$$d_{i,j}({u,v})=\mathbb{E}\mathrm{e}_{i+1,0}\mathrm{e}_{i+1,-u}\mathrm{e}_{j+1,0}\mathrm{e}_{j+1,-v},i,j=0,\dots,K-1,u,v=1,\dots.q,$$

$G=(g_{1}|\cdots|g_{K})^{\top}$, $g_{k}=(g_{k,1}^{\top}|\cdots|g_{k,q}^{\top})$ with

$$g_{k,h}=\mathbb{E}e_{k,-h}e_{k,0}\sum_{\ell=1}^{K}e_{\ell,0}\nabla_{\theta_{0}}^{\top}\eta_{\ell,0}(\theta){J^{-1}}^{\top},1\leq h\leq q,$$

and $D$ is the $Kq\times Kq$ matrix with elements given by

$$D[(iq+1):(i+1)q,(jq+1):(j+1)q]=[d_{i,j}(u,v)]_{1\leq u,v\leq q},0\leq i,j\leq K-1,$$

where $\hat{D}$, $\hat{G}$ and $\hat{S}$ are the empirical counterparts of $D$, $G$ and $S.$ Let us set $W=D+C_{q}J^{-1}L{J^{-1}}^{\top}C_{q}^{\top}+GC_{q}^{\top}+C_{q}G^{\top}$ and $\hat{W}$ as its empirical counterpart.

We require the following additional assumption, which represents the invertibility of $W.$ For $k=1,\ldots,K$, let us set

$$\iota_{k,t}=\left(0_{q}^{\top},\ldots,0_{q}^{\top},\underbrace{e_{k,t-1}(\theta_{0}),\ldots,e_{k,t-q}(\theta_{0})}_{\text{k-th block of q elements}},0_{q}^{\top},\ldots,0_{q}^{\top}\right)^{\top}.$$

INV-1

For $\lambda\in\mathbb{R}^{Kq},$ if

$$\lambda^{\top}\sum_{k=1}^{K}e_{k,0}\left(\iota_{k,0}+C_{q}J^{-1}\nabla_{\theta_{0}}\eta_{k,0}(\theta)\right)=0\text{ a.s.},$$

then $\lambda=0.$

It is well known that Portemanteau tests are omnibus tests (Francq and Zakoian,, 2019). They are less powerful than specific alternative tests that correspond to the nullity of the coefficient of the extra-lag model. Our study focuses only on a model of lag $1.$

When we are interested in two or more sites for the same process, it may happen that we have the same models for certain sites with coefficient values that are a priori different. We propose a two-by-two model comparison test. For example, when we are interested in two ecological sites, we assume that the model on one site is driven by $\theta_{0}^{(1)}$ and that of the second site is driven by $\theta_{0}^{(2)}$. We want to test whether the $p-th$ covariate has the same effect on both sites, that is $\theta_{p,0}^{(1)}=\theta_{p,0}^{(2)}$ , and also globally $\theta_{0}^{(1)}=\theta_{0}^{(2)}$. For $j=1,2$, let us denote $(Y_{t}^{(j)},X_{t}^{(j)})_{t\in\mathbb{Z}}$ as the processes $(Y_{t},X_{t})$ coinciding with the site $j$. The issue here is testing the equality of the distributions of $(Y_{t}^{(1)},X_{t}^{(1)})_{t\in\mathbb{Z}}$ and $(Y_{t}^{(2)},X_{t}^{(2)})_{t\in\mathbb{Z}}$. Tests of the laws of two strictly stationary processes have been investigated previously. They consist of tests on: 1. marginal distributions (Doukhan et al.,, 2020, 2015); 2. jointly finite distributions (Dhar et al.,, 2014); and 3. all possible $d-$dimensional joint distributions of both processes (Pommeret et al.,, 2022). With regards to the parametric form of Equation (2.1)-(2.2), we can perform the test of equality of laws of entire processes without requiring their finite distributions.

For $\theta\in\Theta$, $(s_{t}^{(j)}(\theta))_{t\in\mathbb{Z}}$ and $(h_{t}^{(j)}(\theta))_{t\in\mathbb{Z}}$ are respectively the score and Hessian processes evaluated on $\theta$ for the process $(Y_{t}^{(j)},X_{t}^{(j)})_{t\in\mathbb{Z}}.$ Under the conditions of Proposition 2 for $(Y_{t}^{(j)},X_{t}^{(j)})_{t\in\mathbb{Z}},$

$$\lim_{n\to\infty}\frac{1}{\sqrt{n}}\sum_{t=1}^{n}s_{t}^{(j)}(\theta_{0}^{(j)})=\mathcal{N}(0,L_{j})$$

and a.s.

$$\lim_{n\to\infty}\frac{1}{n}\sum_{t=1}^{n}h_{t}^{(j)}(\theta_{0}^{(j)})=J_{j},j=1,2.$$

Let us set $S=\mathbb{E}s_{0}^{(2)}(\theta_{0}^{(2)})s_{0}^{(1)}(\theta_{0}^{(1)})^{\top}$, $V=J_{1}^{-1}L_{1}{J_{1}^{-1}}^{\top}+J_{2}^{-1}L_{2}{J_{2}^{-1}}^{\top}+J_{2}^{-1}S{J_{1}^{-1}}^{\top}+{J_{1}^{-1}}S^{\top}{J_{2}^{-1}}^{\top}$ and $\hat{V}$ its empirical counterpart. We denote $\hat{\theta}_{n}^{(1)}$ (resp. $\hat{\theta}_{n}^{(2)}$) the estimator of $\theta_{0}^{(1)}$ (resp. $\theta_{0}^{(2)}$), given by Equation (2.6) , and $\hat{\theta}_{p,n}^{(1)}$ (resp. $\hat{\theta}_{p,n}^{(2)}$) as the $p-th$ element of $\theta_{0}^{(1)}$ (resp $\theta_{0}^{(2)}$).

We require the following additional assumption for the invertibility of $V$.

INV-2

For $\lambda\in\mathbb{R}^{3K+P},$ if

$$\lambda^{\top}\left(J_{1}^{-1}s_{t}^{(1)}(\theta_{0}^{(1)})+J_{1}^{-2}s_{t}^{(2)}(\theta_{0}^{(2)})\right)=0\text{ a.s.},$$

then $\lambda=0.$

The numerical implementation of Model (2.1)-(2.2) and the statistical inference methods are implemented using R (R Core Team,, 2022). We choose $\beta_{j,0}=0.5,j=1,\ldots,K$ for the initial values of the latent processes, whereas the first observation of the paths acts to initialize the sequence $(Y_{t})_{t\in\mathbb{Z}}$. The initial value of the score sequence is chosen to be equal to $0_{3K+P}.$ The loss function (Equation (2.6)) is optimized for $20$ different random initializations using the L-BFGS-B method (Byrd et al.,, 1995), a derivative-free optimization routine that allows box constraints, run in the R function optim. The estimated $\theta$ values are those corresponding to the minimum loss function value over all initializations. Let us set $\theta_{j}$ as the $j-th$ element of the parameters vector $\theta$, and the set $\Theta\ni\theta$ is

$$\Theta=:\Theta_{\varepsilon}=\left\{\theta\in\mathbb{R}^{3K+P}:\theta_{3K+P-j}\in[-1+\varepsilon,1-\varepsilon],j=1,2,3\text{ and }\theta_{j}\in[-1/\varepsilon,1/\varepsilon],j=1,\ldots,2K+P\right\}$$

for a fixed and known $\varepsilon>0$. Here we fix $\varepsilon=1\mathrm{e}^{-6}$. All methods are made available using two top level functions poly_autoregression and comparison_test_dependent. The first takes as an argument the observed paths of sequences $(\overline{Y}_{t})_{t\in\mathbb{Z}}$ and $(X_{t})_{t\in\mathbb{Z}}$ and returns, among others, the estimated values of parameters $\theta$, its standard values and the result of the Portemanteau test in Proposition 3. The second takes as arguments the returned values of two models fitted with the poly_autoregression function and gives the statistics and the p-values of the test in Proposition 4.

Recall that the vector of the parameters is denoted by $\theta=(\omega^{\top},\mathrm{vec}(\Gamma)^{\top},\alpha^{\top},\beta_{1},\ldots,\beta_{K})^{\top}$. For numerical experimentation, we use a simulated adjacent categorical response variable with four levels, as for our real data example, and five simulated covariates. We consider the sets of parameters of Table 1, chosen to meet the condition $|\beta_{j}|<1,j=1,\ldots K$ in Proposition 1.