21.1 FORECASTING COUNT DATA

21.1 FORECASTING COUNT DATA

Often we are interested in forecasting outcomes for variables restricted to take on integer values. For example, ${\gamma }_t$ could have a support $\{0,1,2,\dots ,M\}$ , where M could be infinite. In such situations it seems appropriate to take this support into account both for the elicitation of the loss function as well as in the construction of the forecasting model. It would seem as though nearly all forecasting problems fit this situation; for example, forecasting the number of jobs added to the economy fits this situation since the number of jobs is a natural number. This also holds when forecasting the level of GDP since this is often rounded to millions of dollars and in any case cannot be fractions of cents. However, for large values of $M,$ the “count data” nature of the problem can appropriately be ignored because taking it into consideration considerably complicates the forecasting process with no hope of any gains to either predictive accuracy or our understanding of the problem. Hence, applications of count data really refer to situations where M is quite small, typically no more than 5 or 10. Examples include forecasting the size of a household, or the number of cars or computers owned by the household, or the number of individual stocks held by private investors.

There appears to be very little work done on understanding loss functions for count data problems. Nearly all point forecasting methods revolve around the provision of estimates of the conditional mean, implicitly assuming a squared error loss function.

21.1.1 Parametric Models for Count Data

Parametric approaches to modeling count data choose distributions for which the support includes only nonnegative integers. The parametric model of choice in this literature is the Poisson model. The Poisson(λ), $\lambda \ge 0$ model has support $y=$ $\{0,1,2,\dots \}$ , so M is infinite and assigns probabilities $P\left[Y=y\right]=e^{-\lambda }{\lambda }^y\overline{/(y!)}$ . Both the mean and variance of the Poisson distribution equal λ and all moments exist for this distribution.

More generally, in a regression setting we could allow the Poisson parameter λ to be a function of conditioning variables, $\lambda =\lambda (z)$ , once again imposing that $\lambda (z)$ is positive. The standard approach to modeling λ(z) is to use the exponent of a linear index,

$$\begin{array}{r}\lambda (z)=\exp \left({\beta }^{\prime }z\right).\end{array}\text{\hspace{1em}(21.1)}$$

Assuming that the parametric model is correctly specified, maximum likelihood procedures are available. Results are also available for other estimators. The literature on estimation of these models is specialized and lengthy; Cameron and Trivedi (1998) provide an excellent overview.

Models of the type in (21.1) can include lags of the counts ${\gamma }_t$ in $z_t$ . In the Poisson model this would also make the variance of the outcome depend on lags of the counts. Alternatively, a wide set of dynamic models are aimed at generalizing autoregressive models. McKenzie (2003) provides a wide ranging review of such models.

A popular generalization of autoregressive specifications for use in modeling and forecasting count data is the integer autoregressive model of order $P,$ denoted $\mathrm{INAR}(p)$ , introduced by Al-Osh and Alzaid (1987). The most common version of this model takes the form

$$y_t=a_1\circ y_{t-1}+a_2\circ y_{t-2}+\cdot \cdot \cdot +a_{\mathcal{P}}\circ y_{t-k}+{\epsilon }_t,\text{\hspace{1em}(21.2)}$$

where $\begin{array}{r}{a}_i\circ y_{t-i}={\sum }_{j=1}^{y_{t-i}}{B}_{j,i}\end{array}$ , where each of the $B_{j,i}$ terms are Bernoulli(ai ) 0–1 random variables with $0<a_i<1$ . The motivation for this approach comes from $a_i\circ y_{t-i}$ always being integer valued—as opposed to, say, including a lagged dependent variable with a fixed coefficient. Typically, ${\epsilon }_t$ is drawn from a parametric distribution that is nonnegative and independent of $B_{j,i}$ and past values of $\gamma .$ . This ensures that ${\epsilon }_t$ is independent of the lagged component of the model, which is integer valued, and so requires ${\epsilon }_t$ to also be integer valued. Again, the most common specification for ${\epsilon }_t$ is the Poisson(λ) distribution, although other distributions such as the binomial are also used. The distribution for ${\epsilon }_t$ is nonnegative and so has a nonzero mean, which explains the absence of a constant term in the model (21.2). Using that $y_{t-i}$ and $B_{j,i}$ are independent, we have $E\left[a_i\circ y_{t-i}|y_{t-i}\right]=a_i{y}_{t-i}$ . As $a_i<1$ , past values for the count have a damped effect on the current count.1 For stationary models we require that ${\sum }_{i=1}^p{a}_i<\bar{1}$ . Numerous variations on this model result in different properties; see McKenzie (2003).

Estimation of the INAR(1) model with $T$ observations on ${\gamma }_t$ is relatively straightforward. Assuming ${\epsilon }_t$ to be Poisson(λ), the parameters $(a,\lambda )$ can be estimated using either maximum likelihood or simply by OLS regression of ${\gamma }_t$ on a constant and ${\gamma }_{t-1}$ with aˆ being the OLS coefficient on the lagged dependent variable, followed by setting $\begin{array}{r}\hat{\lambda }=(T-1{)}^{-1}{\sum }_{t=2}^T(y_t-\hat{a}{y}_{t-1})\end{array}$ . Using these plug-in estimators, a prediction of the mean of ${\gamma }_{T+1}$ can be computed as $\hat{\alpha }{y}_T+\hat{\lambda }$ , where the extra term arises through the mean of ${\epsilon }_T$ . Alternatively, since the model is fully parametric, we can construct the forecast distribution directly (Freeland and McCabe, 2004). This approach does not generally result in an integer-valued prediction, although the forecast can be rounded to the nearest integer value.

Using recursive substitution and taking expectations results in the h-step-ahead forecast,2

$$E\left[{\gamma }_{t+h}|{\gamma }_t\right]=a^h{y}_t+\frac{1-a^h}{1-a}\lambda ,\text{\hspace{1em}(21.3)}$$

where the first term is familiar from the AR(1) model and the second term arises due to the nonzero mean of ${\epsilon }_t$

To construct the forecast distribution directly, consider the INAR(1) model with Poisson(λ) residuals:

$$\begin{array}{l}P\left(Y_{T+h}=y|y_T\right)\\ =\sum _{s=0}^{\min (y,y_T)}\left(\frac{y_T}{s}\right)a^{hs}(1-a^h{)}^{y_T-s}\frac{1}{(y-s)!}\exp \left(-\lambda \frac{1-a^h}{1-a}\right){\left(\lambda \frac{1-a^h}{1-a}\right)}^{y-s}\end{array}\text{\hspace{1em}(21.4)}$$

for each $\gamma =0,1,2,\dots$ . To derive this expression, note that there are a number of combinations of the autoregressive component and the residual that sum to $y.$ For each of these combinations, the probability is simply the product of the binomial and Poisson terms that give rise to that sum. This follows from the independence of these terms. If the binomial part equals some value $\tilde{y}\le y$ , then the Poisson part must be $(\mathbf{y}-\tilde{\mathbf{y}})$ for the total sum to be $\gamma .$ Plug-in estimates of $(\hat{a},\hat{\lambda })$ can be used to generate an estimated h-step-ahead forecast distribution and, provided these parameters are consistently estimated, the probabilities constructed in this manner are also consistent. Bayesian prediction methods have also been developed for this class of models; see McCabe and Martin (2005).

Less parametric approaches have also been considered. McCabe, Martin, and Harris (2011) study the $\mathrm{INAR}(p)$ model (21.2) but leave the distribution of ${\epsilon }_t$ unrestricted apart from imposing that its support is the set of nonnegative integers. Instead they suggest a nonparametric MLE approach for estimating the distribution of the residuals; for details of this method, see section 2.3 of their paper.