For modeling integer-valued time series (or any non-Gaussian time series), I would opt for a State-Space model that allows the latent dynamic process to evolve independently of the observations. This is very useful when dealing with observations that have restrictions, such as non-negative integers or proportional values, because we can let the latent dynamic process be real-valued and use the convenience of link functions (like we do in Generalized Linear Models) to translate to the observation scale. My package {mvgam} was designed specifically for this kind of situation because I frequently have to analyse and forecast multivariate sets of count-valued time series with missing values, many zeros and overdispersion, and none of the more commonly used methods (such as the INAR) are capable of dealing with all of these features. I also wanted the ability to include nonlinear smooth functions of covariates (using Generalized Additive Models) in both the latent process model and in the observation model, because again many real-world time series have observation error that needs to be captured.
The general formula for Dynamic Generalized Additive Models is:
$$for~i~in~1:N_{series}~...$$
$$for~t~in~1:N_{timepoints}~...$$
$$g^{-1}(\tilde{\boldsymbol{y}}_{i,t})=\alpha_{i}+\sum\limits_{j=1}^J\boldsymbol{s}_{i,j,t}\boldsymbol{x}_{j,t}+Z\boldsymbol{z}_{i,t}\,,$$
Here $\alpha$ are the unknown intercepts, the $\boldsymbol{s}$'s are unknown smooth functions of covariates $(\boldsymbol{x})$'s, which can potentially vary among the response series, and $\boldsymbol{z}$ are dynamic latent processes. Each smooth function $\boldsymbol{s_j}$ is composed of basis expansions whose coefficients, which must be estimated, control the functional relationship between $\boldsymbol{x}_{j}$ and $g^{-1}(\tilde{\boldsymbol{y}})$. But note that we can also include smooth functions of covariates in the process models $\boldsymbol{z}$. The size of the basis expansion limits the smooth’s potential complexity, where a larger set of basis functions allows greater flexibility.
To capture time series dynamics, we can impose a wide variety of temporal dynamic structures (such as Random Walk, AR processes, Continuous Time AR processes, or even Vector Autoregressions). The $Z$
matrix affords even more flexibility by letting some series share the same latent process model (i.e. perhaps two observation series are tracking the same hidden process, but with different observation errors). For more information on GAMs and how they can smooth through data, see this blogpost on how to interpret nonlinear effects from Generalized Additive Models.
To see how this can work for integer-valued time series, you can see a recent webinar with a short worked example.
