FINANCE , STATISTICS, FORECASTING. ACCURACY TESTING AND ERROR MATRICES

Omega 40 (2012) 703–712

Contents lists available at SciVerse ScienceDirect

Omega

0305-04

doi:10.1

n Corr

E-m

john.bo

journal homepage: www.elsevier.com/locate/omega

Forecast horizon aggregation in integer autoregressive moving average (INARMA) models

Maryam Mohammadipour n, John E. Boylan

Buckinghamshire New University, High Wycombe Campus, Queen Alexandra Road, High Wycombe, Buckinghamshire, HP11 2JZ, UK

a r t i c l e i n f o

Available online 21 October 2011

Keywords:

Discrete time series

INARMA model

Temporal aggregation

Cross-sectional aggregation

Forecast horizon aggregation

Yule�Walker estimation

83/$ - see front matter & 2011 Elsevier Ltd. A

016/j.omega.2011.08.008

esponding author. Tel.: þ44 7758357575.

ail addresses: [email protected] (M. Mo

[email protected] (J.E. Boylan).

a b s t r a c t

This paper addresses aggregation in integer autoregressive moving average (INARMA) models.

Although aggregation in continuous-valued time series has been widely discussed, the same is not

true for integer-valued time series. Forecast horizon aggregation is addressed in this paper. It is shown

that the overlapping forecast horizon aggregation of an INARMA process results in an INARMA process.

The conditional expected value of the aggregated process is also derived for use in forecasting.

A simulation experiment is conducted to assess the accuracy of the forecasts produced by the

aggregation method and to compare it to the accuracy of cumulative h-step ahead forecasts over the

forecasting horizon. The results of an empirical analysis are also provided.

& 2011 Elsevier Ltd. All rights reserved.

1. Introduction

Time series aggregation is a widely discussed subject for contin- uous-valued time series. It goes back over 50 years [1] and since then many papers have considered different aspects of aggregation for continuous-valued time series (see for example [2–9]). Three types of aggregation have been identified in the literature which can be classified as cross-sectional aggregation, temporal aggregation and forecast horizon aggregation. The first class of aggregation produces forecasts based on aggregated data series and the other two produce forecasts based on aggregated periods.

Cross-sectional or contemporaneous aggregation is conducted across individual series rather than time. For example, in demand forecasting of many products with a short demand history, similar products are grouped in a product family and the demand forecast is built for the family rather than individuals, which may produce more reliable forecasts than the forecasts for individual items.

Temporal aggregation, also called flow scheme, refers to aggregation in which a low frequency time series (e.g. annual) is derived from a high frequency time series (e.g. quarterly or monthly). The low frequency variable is the sum of k consecutive periods of the high frequency variable. For example, the annual observations are the sum of the monthly observations every twelve periods (k¼12).

Finally, forecast horizon aggregation refers to the case in which there is requirement for a forecast of the total value over a number of time periods ahead. For example, in demand forecasting in a supply

ll rights reserved.

hammadipour),

chain, when there is a lead time between ordering by a manufacturer and receiving the order from a supplier, the demand over that lead time has to be forecasted in order to prevent shortage during the lead time period (see for example [10]).

With respect to temporal/forecast horizon aggregation, we must distinguish between overlapping and non-overlapping cases. In non- overlapping aggregation, the demand series are divided into con- secutive non-overlapping blocks of equal length. In overlapping aggregation, the blocks are of equal lengths but, at each period, the oldest observation is dropped and the newest is included. This paper focuses on the case of overlapping aggregation.

Although many papers examine continuous-valued time series, issues related to their application [11,12] and different types of aggregation in them [2–9], the same is not true for time series of counts. Brännäs et al. [13] first studied temporal and cross-sectional aggregation of an Integer Auto-Regressive process of order one, INAR(1). To our knowledge, forecast horizon aggregation in more general Integer Auto-Regressive Moving Average (INARMA) models has not been studied before. This is the motivation for this study to begin to address this issue and present some new results. The paper is structured as follows. The forecast horizon aggregation of INAR- MA(p,q) processes is discussed in detail in Section 2. A simulation experiment is designed and performed in Section 3 to assess the accuracy of the aggregated forecasts of Section 2. An empirical analysis, based on two datasets, is performed in Section 4. The conclusions are provided in the final section of the paper.

2. Forecast horizon aggregation and forecasting

In this section, aggregation and forecasting over a forecast horizon are discussed. This has applications in many areas. Some

M. Mohammadipour, J.E. Boylan / Omega 40 (2012) 703–712704

application areas require forecasts of the whole distribution [14,15]. However, other application areas need forecasts of the conditional mean [13,16]. This study concentrates on estimation of the conditional mean, while further research will focus on forecasting the whole distribution. This research examines whether there is any benefit to be leveraged from INARMA models in forecasting the conditional mean. It is shown that there is such a benefit in certain circumstances.

First, the conditional expected values of the aggregated INAR(1) and INMA(1) processes are presented. Then, it is shown that the forecast horizon aggregation of an INARMA process is an INARMA process. The conditional expected value of the aggre- gated INARMA process is also obtained.

The most common forecasting procedure discussed in the time series literature is using the conditional expectation. The main advantage of this method, apart from being simple, is that it produces forecasts with minimum mean square error (MMSE). This forecasting procedure is adopted in this paper.

2.1. INAR(1) models

A Poisson INAR(1) process, PoINAR(1) is defined by

Y t ¼a3Y t�1þZt ð1Þ

where aA(0,1] and {Zt} is a sequence of i.i.d. non-negative integer- valued Poisson distributed random variables, with mean and finite variance l. Zt and Yt�1 are assumed to be stochastically independent for all points in time. The thinning operation ‘‘ ’’ of Sueutel and van Harn [17] is defined by a3Y ¼

PY i ¼ 1 Xi where {Xi}

is a sequence of i.i.d. Bernoulli random variables with P(Xi¼1)¼a for i¼1,y,Y.

It follows from [14] that the conditional mean of the aggre- gated process is

E Xl i ¼ 1

Y tþi9Y t

! ¼ að1�alÞ

1�a Y tþ

l 1�a

l� Xl j ¼ 1

aj 0 @

1 A ð2Þ

As shown in Appendix A, the conditional variance of the aggregated process is as follows:

var Xl i ¼ 1

Y tþi9Y t

" # ¼Y t

Xl j ¼ 1

ajð1�ajÞþ l

1�a l� Xl j ¼ 1

aj 2 4

3 5

þ 2l

1�a Xl j ¼ 1

a2j�1 ðl�jÞ� að1�al�jÞ

1�a

� � ð3Þ

For an INAR(1) process of (1), the cumulative value over horizon l is given by

Xl j ¼ 1

Y tþj ¼Y tþ1þY tþ2þ���þY tþ l ¼ða3Y tþZtþ1Þ

þða23Y tþa3Ztþ1þZtþ2Þþ �� �þðal3Y tþal�13Ztþ1 þal�23Ztþ2þ ���þZtþ lÞ ð4Þ

Bearing in mind that a3Xþb3X aðaþbÞ3X (the LHS is the sum of two binomial random variables with the same number of trials and different success probabilities), the above equation can be written in the following form:

Xl j ¼ 1

Y tþj ¼ Xl j ¼ 1

Xn1j i ¼ 1

c1ij3Y tþ Xl j ¼ 1

Xn2j i ¼ 1

c2ij3Ztþkij ð5Þ

where n1 j

is the number of Yt terms in each of fY tþjg l j ¼ 1

in (4), c1ij is the corresponding coefficient for each Yt, n

2 j

is the number of

Ztþkij terms in each of fY tþjg l j ¼ 1

in (4) and c2ij is the corresponding

coefficient for each Ztþkij. Further details about the coefficients are given by [18].

It can be seen that, based on (5), the conditional expected value of the aggregated PoINAR(1) process is

E Xl j ¼ 1

Y tþ j9Y t

0 @

1 A¼ Xl

j ¼ 1

Xn1j i ¼ 1

c1ij

0 @

1 AY tþ Xl

j ¼ 1

Xn2j i ¼ 1

c2ij

0 @

1 Al ð6Þ

The above equation is the same as (2). At time T, when YT is observed, the aggregated forecast can be obtained from

E Xl j ¼ 1

Y T þj9Y T

0 @

1 A¼ að1�alÞ

1�a Y T þ

l 1�a

l� Xl j ¼ 1

aj 2 4

3 5 ð7Þ

2.2. INMA(1) models

For an INMA(1) process of Y t ¼b3Zt�1þZt , where bA(0,1] and {Zt} is as before, the cumulative value over horizon l is given by

Xl j ¼ 1

Y tþ j ¼Y tþ1þY tþ2þ���þY tþ l ¼ðb3ZtþZtþ1Þþðb3Ztþ1þZtþ2Þ

þ ���þðb3Ztþ l�1þZtþ lÞ ð8Þ

The above equation can be written in the following form:

Xl j ¼ 1

Y tþ j ¼ Xl j ¼ 1

Xnj i ¼ 1

cij3Ztþkij ð9Þ

where nj is the number of Ztþkij terms in each of fY tþjg lþ1 j ¼ 1

and is the corresponding coefficient for each Ztþkij. Further details about the coefficients are given by [18].

Based on (9), the conditional expected value of the aggregated INMA(1) process is

E Xl j ¼ 1

Y tþ j9Y t

0 @

1 A¼ Xl

j ¼ 1

X2 i ¼ 1

cij

0 @

1 Al¼ Xl

j ¼ 1

ð1þbÞ

0 @

1 Al¼ lð1þbÞl

ð10Þ

2.3. INARMA(p,q) models

This paper examines aggregation and forecasting of a general INARMA process over a forecast horizon. The INARMA(p,q) pro- cess is given by

Y t ¼ Xp i ¼ 1

ai3Y t�iþZtþ Xq j ¼ 1

bj3Zt�j ð11Þ

where a1,y,ap�1A[0,1],apA(0,1]; b1,y,bq�1A[0,1],bqA(0,1] and {Zt} is a sequence of i.i.d. non-negative integer-valued random variables, independent of Yt with mean mZ and finite variance s2Z . The thinning operations are defined as follows:

a3Y ¼ XY i ¼ 1

Xi ð12Þ

where {Xi} is a sequence of i.i.d. Bernoulli random variables with P(Xi¼1)¼a for i¼1,y,Y. This paper follows the approach of Du and Li [19] and McKenzie [20] regarding the binomial thinning mechanisms for INAR(p) and INMA(q), respectively. Therefore, it is assumed that the individual thinning operations ai3Y t�i for i¼1,y,p and bj3Zt�j for j¼1,y,q are performed independently not only from each other, but also from corresponding operations at previous times in (11).

M. Mohammadipour, J.E. Boylan / Omega 40 (2012) 703–712 705

The stationarity conditions of this process are the same as those of an INAR(p) process. Neal and Rao [21] suggest that the invertibility conditions for this process are the same as the those of an MA(q) process (

Pq j ¼ 1

bj o1). The MMSE one-step-ahead forecast for an INARMA(p,q) pro-

cess of (11) is

Ŷ T þ1 ¼a1Y T þ���þapY T�pþ1þlþb1ZT þ���þbqZT�qþ1 ð13Þ

The h-step ahead forecast when hrq is

Ŷ T þh ¼a1 Y T þh�1þ���þap Y T þh�pþlþbh ZT þ���þbqZT þh�q þlðb1þ. . .þbh�1Þ ð14Þ

where the Y values on the RHS of (14) may be either actual or forecast values. When h4q, the h-step ahead forecast becomes

Ŷ T þh ¼a1 Y T þh�1þ���þap Y T þh�pþl Xq j ¼ 0

bj ð15Þ

where again the Y values on the RHS of the above equation may be either actual or forecast values and b0¼1.

We next present two propositions regarding the aggregation and forecasting of an INARMA(p,q) process.

Proposition 1. Aggregation of an INARMA(p,q) process over a forecast horizon results in an INARMA(p,q) process.

Proof. For an INARMA(p,q) process of (11), the aggregated pro- cess over a forecast horizon can be written as

Xl j ¼ 1

Y tþj ¼ Xl j ¼ 1

Xp i ¼ 1

ai3Y tþ j�i

" # þZtþjþ

Xq i ¼ 1

bi3Ztþj�i

" #( )

¼ Xp i ¼ 1

ai3 Xl j ¼ 1

Y tþ j�i

2 4

3 5þ Xl

j ¼ 1

Ztþjþ Xq i ¼ 1

bi3 Xl j ¼ 1

Ztþ j�i

2 4

3 5 ð16Þ

Now, if we assume that Pl

j ¼ 1 Y tþj ¼Yt and Pl

j ¼ 1 Ztþ j ¼ Zt,

(16) can be written as

Yt ¼ Xp i ¼ 1

ai3Yt�iþZtþ Xq i ¼ 1

bi3Zt�i ð17Þ

which is also an INARMA(p,q) process. Therefore, aggregation of

an INARMA(p,q) process over a forecast horizon results in an

INARMA(p,q) process with the same INAR and INMA parameters

but with a different innovation parameter. If Zt�Po(l), Zt will be the sum of l independent Poisson variables; thus, Zt�Po(ll).

Proposition 2. The forecast horizon aggregated INARMA(p,q) process can be written in terms of the last p observations as follows:

Xl j ¼ 1

Y tþj ¼ Xl j ¼ 1

Xn1j i ¼ 1

c1ij3Y tþ Xl j ¼ 1

Xn2j i ¼ 1

c2ij3Y t�1þ���

þ Xl j ¼ 1

Xnpj i ¼ 1

cp ij 3Y t�pþ1þ

Xl j ¼ 1

Xnpþ1j i ¼ 1

cpþ1 ij

3Ztþkij ð18Þ

with the parameters as shown in Table 1 (see Appendix B for the proof).

The conditional expected value of the aggregated process given the p-previous observations can then be obtained from

E Xl j ¼ 1

Y tþj9Y t�pþ1 ,. . .,Y t�1 ,Y t

0 @

1 A¼ Xl

j ¼ 1

Xn1j i ¼ 1

c1ij

0 @

1 AY t

þ Xl j ¼ 1

Xn2j i ¼ 1

c2ij

0 @

1 AY t�1þ���þ Xl

j ¼ 1

Xnpj i ¼ 1

cp ij

0 @

1 AY t�pþ1

þ Xl j ¼ 1

Xnpþ 1j i ¼ 1

cpþ1 ij

0 B@

1 CAl ð19Þ

The above equation is then used to forecast the aggregated process. In the next section, the accuracy of such aggregated forecasts will be assessed for INAR(1), INMA(1) and INARMA(1,1) processes.

3. Simulation

In order to test the benefit of using the INARMA horizon- aggregated forecasts for short and long histories, Monte Carlo simulations are conducted. The results for the INARMA horizon- aggregated forecasts (hereafter abbreviated as INARMA-Agg) have been compared to the results of using cumulative h-step ahead forecasts over the horizon (hereafter abbreviated as INARMA-h). The latter is given by

Pl i ¼ 1 Ŷ tþ i, where Ŷ tþi is the i-step ahead

forecast. It should be emphasized that there are two approaches in the literature regarding h-step ahead forecasting. The first approach used by Brännäs and Hellström [22] is based on repeated substitution of the INARMA process. For example, the h-step ahead forecast of an INAR(1) process can be obtained from

Y T þh ¼ah3Y T þ Xh i ¼ 1

ah�i3ZT þi ð20Þ

It can be easily shown that aggregation of (20) over a horizon results in (7). Therefore, this approach results in the same aggregated forecast as that proposed by this study.

In the second approach, the Y value on the right hand side of the equation Ŷ T þh ¼aŶ T þh�1þl is a forecast [19,23]. The h-step ahead forecasts in this section are calculated based on this approach. For simulation purposes, it is assumed that the innova- tions, {Zt}, have a Poisson distribution with parameter l. Although this assumption is restrictive, a large number of data series of the empirical datasets used in this study met the above condition (see Section 4). This is consistent with the larger empirical study by Eaves [24] discussed in Section 4. The theoretical findings in this paper, however, are not based on any distributional assumptions and can be used as a framework for future studies based on other marginal distributions.

Three INARMA processes are considered in the simulation experiment: INAR(1), INMA(1) and INARMA(1,1). The aggregated forecasts for these models can be derived from Section 2. For an INAR(1) process the aggregated forecast is provided by (2). For the INMA(1) and INARMA(1,1) processes, these are given by

E Xl i ¼ 1

Y T þi9Y T

" # ¼ lð1þbÞl ð21Þ

E Xl i ¼ 1

Y T þi9Y T

" # ¼ að1�alÞ

1�a Y T þ

lð1þbÞ 1�a

l� Xl j ¼ 1

aj 2 4

3 5 ð22Þ

The control parameters to be varied are a, b, l and n (number of historical observations). The range of these parameters is given in Table 2. Different lengths of series are considered in the Monte Carlo simulations to test the sensitivity of the results to the length of history. In real cases, we are often restricted by short lengths of history (as will be seen in the empirical analysis). Hence, we use n¼24,36,48,96,500 to encompass short data histories as well as long. The forecast horizons considered are three, six and nine periods. The number of replications is set to 1000. This is

Table 1 Parameters of the forecast horizon aggregated INARMA(p,q) model.

for w¼1,y,p nwj ¼

Xp i ¼ 1

nwj�i

! þ1 jrp�ðw�1Þ

Xp i ¼ 1

nwj�i j4p�ðw�1Þ

8>>>>>< >>>>>:

cwij ¼

apc w iðj�pÞ i ¼1,. . .,n

w j�p

^ ^

a1c w iðj�1Þ i ¼n

w j�2 þ1,. . .,nw

j�2 þnw

j�1

ajþðw�1Þ i ¼nwj�1þ1

jrp�ðw�1Þ

apc w iðj�pÞ i¼1,. . .,n

w j�p

^ ^

a1c w iðj�1Þ i¼n

w j�2 þ1,. . .,nw

j�2 þnw

j�1

j4p�ðw�1Þ

2 664

2 666664

8>>>>>< >>>>>:

n pþ1 j ¼

Xp i ¼ 1

n pþ1 j�i

! þðqþ1Þ

cpþ1 ij ¼

apc pþ1 iðj�pÞ

i¼1,. . .,n pþ1 j�p

^ ^

a1c pþ1 iðj�1Þ

i¼n pþ1 j�2 þ1,. . .,n

pþ1 j�2 þn

pþ1 j�1

bq ,. . .,b1 ,1 i¼n pþ1 j�1 þ1,. . .,n

pþ1 j�1 þn

pþ1 j

8>>>>>>< >>>>>>:

kij ¼

fkiðj�pÞg i¼1,. . .,n pþ1 j�p

^ ^

fkiðj�1Þg i¼ Xp z ¼ 2

n pþ1 j�z þ1,. . .,ð

Xp z ¼ 2

n pþ1 j�z Þþn

pþ1 j�1

j�q,. . .,j�1,j i¼ Xp z ¼ 1

n pþ1 j�z þ1,. . .,n

pþ1 j

8>>>>>>>>>>< >>>>>>>>>>:

Table 2 Range of control parameters.

Number of observations n¼24,36,48,96,500 INAR(1) a¼0:1,l¼0:5 a¼0:1,l¼1 a¼0:1,l¼3 a¼0:1,l¼5

a¼0:5,l¼0:5 a¼0:5,l¼1 a¼0:5,l¼3 a¼0:5,l¼5

INMA(1) b¼0:1,l¼0:5 b¼0:1,l¼3 b¼0:1,l¼5 b¼0:5,l¼0:5 b¼0:5,l¼1 b¼0:5,l¼3 b¼0:5,l¼5 b¼0:9,l¼0:5 b¼0:9,l¼1 b¼0:9,l¼3 b¼0:9,l¼5

INARMA(1,1) a¼0:1,b¼0:1,l¼0:5 a¼0:1,b¼0:1,l¼1 a¼0:1,b¼0:1,l¼5 a¼0:1,b¼0:9,l¼0:5 a¼0:1,b¼0:9,l¼1 a¼0:1,b¼0:9,l¼5 a¼0:1,b¼0:5,l¼0:5 a¼0:5,b¼0:5,l¼1 a¼0:5,b¼0:5,l¼5

M. Mohammadipour, J.E. Boylan / Omega 40 (2012) 703–712706

consistent with other studies of INARMA processes, which used the same or fewer replications (e.g. [22,25–27]) and have been found to give reliable results when compared with findings known from theory.

The data series are divided into two periods: estimation period and performance period. Initialization and estimation of parameters are conducted in the estimation period and the forecasting accuracy is assessed in the performance period. If at least two non-zero values are observed in the estimation period, the first half of the observa- tions is assigned for the estimation period and the other half for the performance period. However, if fewer than two non-zero values are observed in the estimation period, this period will be extended until the second non-zero value is observed.

As an example consider the case of n¼ 24 and l¼3. Under this experimental scenario the length of the estimation period is 12 (if at least two non-zero values are observed, else it is extended until two such values are available). The forecast errors are then calculated in the performance block of periods from period t ¼ 13 (for periods 13, 14 and 15) to period t ¼22 (for periods 22, 23 and 24).

All coding is in MATLAB and random numbers are generated by the poissrnd function for the Poisson distribution and binornd function for the binomial thinning operations.

The autoregressive, moving average and innovation parameters are estimated using the Yule�Walker (YW) estimation method and the estimates are updated in each period. This is a simple estimation method and it has been shown that the YW estimators are asympto- tically equivalent to the Conditional Least Squares (CLS) estimators for an INAR(1) process [28]. Also, for INMA(1) and INARMA(1,1) pro- cesses, the forecasts produced by the two methods are close in terms of Mean Square Error [18].

The forecast accuracy of the two aggregated forecasting methods (INARMA-Agg and INARMA-h) is compared in terms of Mean Square Error (MSE). MSE is a widely used measure in the forecasting literature, is mathematically easy to handle and is a sensible measure for evaluating an individual time series. For methods that produce unbiased forecasts, the MSE coincides with the forecast error var- iance. In the empirical analysis presented in Section 4, the empirical bias properties of the INARMA methods are checked.

Table 3 compares the MSE of the two methods for an INAR(1) process when l¼3. The results for the cases where l¼6,9 are presented in Appendix C.

The results of Table 3 show that for an INAR(1) process, the INARMA-Agg method outperforms the INARMA-h method for all parameter ranges. The improvement increases with the value of l and the length of history. Also, the improvement is higher when the autoregressive parameter is higher because of the nonlinear- ity of the model. For large values of h, the INARMA-h forecasts converge to the unconditional mean of the INAR(1) process:

Ŷ T þh- l

1�a

The larger the value of a, the more variable the data, which makes this convergence less desirable. Some authors have sug- gested using different models for different horizons to improve forecast accuracy [29–31].

Next, the INARMA-Agg and INARMA-h methods are compared for an INMA(1) process when l¼3 (see Appendix C for the cases of l¼ 6,9). It can be seen from Table 4 that for an INMA(1) process, INARMA-Agg method has very slightly better forecasts in terms of MSE than the INARMA-h method. The former is based on (21) and

Table 3 MSEAgg/MSEhof aggregated forecasts for INAR(1) series when l¼3.

Parameters n¼24 n¼36 n¼48 n¼96 n¼500

a¼0:1,l¼0:5 0.9929 0.9937 0.9950 0.9938 0.9960 a¼0:5,l¼0:5 0.8398 0.8527 0.8367 0.8201 0.7885 a¼0:1,l¼1 0.9785 0.9887 0.9924 0.9935 0.9921 a¼0:5,l¼1 0.8278 0.7953 0.7906 0.7410 0.7216 a¼0:1,l¼3 0.9631 0.9604 0.9727 0.9757 0.9835 a¼0:5,l¼3 0.7015 0.6536 0.6328 0.5862 0.5395 a¼0:1,l¼5 0.9231 0.9311 0.9473 0.9614 0.9734 a¼0:5,l¼5 0.5954 0.5512 0.5184 0.4768 0.4397

Table 4 MSEAgg/MSEh of aggregated forecasts for INMA(1) series when l¼3.

Parameters n¼24 n¼36 n¼48 n¼96 n¼500

b¼0:1,l¼0:5 0.9995 0.9996 0.9998 0.9999 1.0000 b¼0:5,l¼0:5 0.9969 0.9981 0.9985 0.9991 0.9998 b¼0:9,l¼0:5 0.9944 0.9961 0.9968 0.9982 0.9996 b¼0:1,l¼1 0.9992 0.9997 0.9998 0.9999 1.0000 b¼0:5,l¼1 0.9971 0.9979 0.9982 0.9991 0.9998 b¼0:9,l¼1 0.9937 0.9958 0.9967 0.9982 0.9996 b¼0:1,l¼3

0.9961 0.9981 0.9982 0.9991 0.9998

b¼0:9,l¼3 0.9932 0.9947 0.9966 0.9981 0.9996 b¼0:1,l¼5 0.9996 0.9993 0.9997 0.9999 1.0000 b¼0:5,l¼5 0.9953 0.9975 0.9984 0.9989 0.9998 b¼0:9,l¼5 0.9929 0.9950 0.9964 0.9979 0.9996

Table 5 MSEAgg/MSEh of aggregated forecasts for INARMA(1,1) series when l¼3.

Parameters n¼24 n¼36 n¼48 n¼96 n¼500

a¼0:1,b¼0:1,l¼0:5 1.1288 1.0764 1.0565 1.0159 0.9976 a¼0:1,b¼0:9,l¼0:5 1.0045 0.9895 0.9661 0.9551 0.9718 a¼0:5,b¼0:5,l¼0:5 0.8688 0.8342 0.8209 0.8145 0.8129 a¼0:1,b¼0:1,l¼1 1.0946 1.0590 1.0474 1.0144 0.9969 a¼0:1,b¼0:9,l¼1 1.0402 0.9692 0.9641 0.9611 0.9735 a¼0:5,b¼0:5,l¼1 0.8646 0.8623 0.8402 0.8285 0.8134 a¼0:1,b¼0:1,l¼5 1.0425 1.0258 1.0292 1.0042 0.9959 a¼0:1,b¼0:9,l¼5 0.9790 0.9680 0.9514 0.9561 0.9639 a¼0:5,b¼0:5,l¼5 0.8652 0.8525 0.8481 0.8252 0.8013

1 This dataset is available from http://www.forecasters.org/ijf/data/Empiri

cal%20Data.xls.

M. Mohammadipour, J.E. Boylan / Omega 40 (2012) 703–712 707

the latter has the following expression:

Xl i ¼ 1

Ŷ T þ i ¼ Xl i ¼ 1

b̂T þiẐT þiþl̂T þi ð23Þ

Comparing (21) and (23) reveals that the only difference between the two forecasting methods is that the INARMA-h method uses the forecast of the innovation term ðẐT þiÞ, while the INARMA-Agg method uses l̂ as an estimate for the innovation term. This difference remains for large samples, because l̂ will converge to a constant (the true), but ẐT þ1, say, would remain a random variable. However, it is expected that the two methods should be very close and the results of Table 4 confirm this.

It can also be seen from Table 4 that with an increase in b, the MSE of the INARMA-Agg method slightly improves compared to that of an INARMA-h method. This could also be attributed to the fact that, for large values of h, the INARMA-h forecasts converge to the unconditional mean of the INMA(1) process l(1þb). Again, larger values of b produce more variable data; therefore, this convergence would result in less accurate forecasts.

Finally, Table 5 compares the MSE of INARMA-Agg and INARMA-h methods for an INARMA(1,1) process when l¼3. The

results confirm the above arguments that when the autoregres- sive parameter is high, the INARMA-Agg method has smaller MSE than the INARMA-h method and the improvement increases with the length of horizon. However, for small autoregressive and moving average parameters, the INARMA-h method is better than the INARMA-Agg. Again, when the number of observation increases the difference decreases.

The above results along with the results of Appendix C suggest that for an INAR(1) process, the INARMA-Agg method outper- forms the INARMA-h method in terms of MSE. The difference between the two methods is high when the autoregressive parameter is high. However, for an INMA(1) process, the two methods produce very close forecasts. With an increase in the moving average parameter the improvement of the INARMA-Agg method over the INARMA-h method slightly increases.

For an INARMA(1,1) process, when the length of horizon is short and the autoregressive and moving average parameters are small, the INARMA-h forecasts have smaller MSEs than INARMA- Agg forecasts. For all the other cases, the latter method beats the former method using MSE.

4. Empirical analysis

In this section, an empirical analysis is conducted to validate the findings on real data. The real demand data series for this research consists of the Royal Air Force (RAF) individual demand histories of 16,000 Stock Keeping Units (SKUs) over a period of 6 years (monthly observations). We have also used another dataset, which consists of 3000 real intermittent demand data series from the automotive industry1 (from [32]), which, unlike the previous one, has more occurrences of positive demand than zeros. This data series consists of demand histories of 3000 SKUs over a period of 2 years (24 months). These two datasets are called Dataset 1 and Dataset 2 from now on.

As previously mentioned, this paper has focused on INARMA processes with Poisson innovations. Although some of the theo- retical results are not based on a distributional assumption, whenever a specific distribution was needed, such as for estima- tion of parameters, a Poisson distribution was assumed.

Out of the four INARMA processes of this study, (INARMA(0,0), INAR(1), INMA(1), and INARMA(1,1)), three of them have a Poisson distribution when the innovation terms are Poisson. The only exception is the INARMA(1,1) process where

varðY tÞ

EðY tÞ ¼

1þaþbþ3ab 1þaþbþab

r1:5 for 0rar1, 0rbr1 ð24Þ

In order to remove the data series with highly variable demands, a Poisson dispersion test (also called the variance test) is needed for all processes except INARMA(1,1). Under the null hypothesis that X1,y,Xn are Poisson distributed, the test statistic:

T CC ¼ Xn i ¼ 1

ðXi�XÞ 2

X ð25Þ

has a chi-square distribution with (n�1) degrees of freedom. Therefore, H0 is rejected if T CC 4w2n�1;1�A.where A is the signifi- cance level. A revised statistic is used to allow for the difference between the mean and variance of an INARMA(1,1) process. The new test statistic is given by

T CCR ¼ T CC 1:5

ð26Þ

Table 6 Information about filtered 16,000 and 3000 datasets.

Min Mean Max Percentage of zeros

Dataset 1: 16,000 series (5,168 filtered series)

0 0.2177 14 85.62

Dataset 2: 3,000 series (1,943 filtered series)

0 2.0194 2 22.77

Table 7 Identification resultsn for Dataset 1 and Dataset 2.

Models INARMA(0,0) INAR(1) INMA(1) INARMA(1,1)

Dataset 1 98.12 0.66 1.04 0.17

Dataset 2 54.55 23.88 17.96 3.60

n In terms of percentage of series.

Table 8 MSEAgg/MSEh of aggregated forecasts for INARMA series when l¼3.

Models INAR(1) INMA(1) INARMA(1,1)

Dataset 1 0.9636 0.9907 1.5557

Dataset 2 0.9706 0.9841 1.1180

Table 9 MSEAgg/MSEh of aggregated forecasts for INARMA series when l¼6.

Models INAR(1) INMA(1) INARMA(1,1)

Dataset 1 0.9057 0.9789 1.6815

Dataset 2 0.9242 0.9573 1.2176

M. Mohammadipour, J.E. Boylan / Omega 40 (2012) 703–712708

The new statistic also has a chi-square distribution with (n�1) degrees of freedom.

The above filtering, with A¼0.05, results in the exclusion of some data series. Out of the 16,000 series, 12,800 series remained and out the 3000 series, 1943 series remained. As mentioned in Section 3, a high percentage of series used in this study can be modeled using the Poisson assumption. This is consistent with the study by Eaves [24], in which over 80% of series had lead-time demand fitting the Poisson distribution at the 5% significance level.

Further filtering of data was performed for series with fewer than two nonzero demands. Out of the 16,000 series, 5168 series met the above criteria and are therefore used for empirical analysis. The filtering of the 3000 series results in 1943 series. It can be seen that although a substantial number of series has the potential to benefit from PoINARMA models; for a large number of series these models are not appropriate. Other distributional assumptions would obviously result in different number of filtered series, which can be pursued as a further study. Relevant characteristics of the filtered datasets are summarized in Table 6.

The design of the empirical analysis follows the detailed simulation design of Section 3. The Yule�Walker estimation method has been used to estimate the parameters of the INARMA models. Two values for forecast horizon have been considered: l¼3,6.

The appropriate INARMA model needs to be identified among the four possible candidates. This is done using a two-stage identification procedure [18]. The first stage distinguishes between the INARMA(0,0) and the other INARMA models. The Ljung�Box statistic of

Q n ¼ nðnþ2Þ Xk j ¼ 1

r̂2j n�j

ð27Þ

is used for this reason. This is a standard test used for conven- tional ARMA models that are included in most software packages (including MATLAB, which is used in this paper) and, based on the argument by Latour [33], it can be used for INARMA models as well. The AIC, as calculated by the formula AIC� Nlogŝ2a þ2m, is then used for identification among the other INARMA models. This is again based on the argument of Latour [33] to use the standard programs for ARMA models for INARMA models. It should also be mentioned that the AIC of ARMA models has been used in the INARMA literature (e.g. [16]).

This identification procedure is applied on our empirical data and the results in terms of the percentage of each of the four INARMA models for each dataset are presented in Table 7.

As shown in Table 7, the great majority of series in Dataset 1 are identified as INARMA(0,0), which was expected due to high number of zeros. The INARMA-Agg and INARMA-h methods produce the same result for INARMA(0,0) series; therefore, these methods are only compared for the other three INARMA models (INAR(1), INMA(1) and INARMA(1,1)). The MSE results are com- pared in Tables 8 and 9 for l¼3 and l¼6, respectively. The bias, in terms of Mean Error, has been checked and found to be low (see [18]).

In order to compare the results with simulation results, the range of estimated parameters for each of the INARMA models is provided in Table 10.

It can be seen that the results for the INAR(1) and INMA(1) processes are in agreement with the simulation results. For the INAR(1) series, the results are comparable to the simulation results of Tables 3 and C1. For the INMA(1) series, this is comparable to the simulation results of Tables 4 and C3. Finally, for the INARMA(1,1) series the results of Dataset 2 are comparable to the simulation results of Tables 5 and C5, which suggest that INARMA-h produces better results than INARMA-Agg for those parameter ranges. It is worth mentioning that for Dataset 1, only a few series were identified as INARMA(1,1).

Therefore, the above results suggest that for an INAR(1) process, the INARMA-Agg method has lower MSE than the INARMA-h method and the improvement increases with the length of history. For an INMA(1) process, the two methods are very close. For INARMA(1,1) series, the empirical results of Dataset 2 confirm the simulation result that, for low autoregressive and moving average parameters, INARMA-h outperforms INARMA-Agg in terms of MSE.

5. Conclusions

This paper addresses forecast horizon aggregation in INARMA processes. The conditional mean of the aggregated process is obtained for the general INARMA(p,q) process, which can be used for forecasting. The purpose of the paper is not to propose new forecasting methods, but rather to compare the performance of alternative INARMA approaches. This has been achieved by simulation and empirical analysis.

It is shown that the aggregation of an INARMA process over a horizon results in an INARMA process. The conditional mean of the aggregated process is also derived as a basis for forecasting. The results of a simulation experiment are provided to assess the accuracy of the forecasts produced using the conditional mean of the aggregated process for three INARMA processes: INAR(1), INMA(1) and INARMA(1,1). The results are compared to the case

Table 10 Parameters’ estimates for Dataset 1 and Dataset 2.

Models INAR(1) INMA(1) INARMA(1,1)

Dataset 1 â is close to 0.2 (the average is 0.2460 and 52.94 percent is between 0.1 and 0.3)

b̂ is close to zero (the average is 0.0898 and 46.29 percent is between 0 and 0.1)

0:1oâo0:3 (the average is 0.2988 and 66.67 percent is between 0.05 and 0.35)

l̂ is around 0.5 (the average is 0.3562 and 97.06 percent is between 0 and 1)

l̂ is around 0.3 (the average is 0.3782 and 55.56 percent is between 0.2 and 0.4)

b̂ is close to zero (the average is 0.1405 and 77.78 percent is between 0 and 0.1)

l̂ is around 0.3 (the average is 0.3558 and 44.44 percent is between 0.2 and 0.5)

Dataset 2 â is close to 0.1 (the average is 0.1234 and 50.65 percent is between 0.05 and 0.15)

b̂ is close to zero (the average is 0.0374 and 79.94 percent is between 0 and 0.05)

0:1oâo0:3 (the average is 0.1907 and 54.28 percent is between 0.05 and 0.35)

l̂ is around 2 (the average is 2.5972 and 40.52 percent is between 1 and 3)

l̂ is between 2 and 3 (the average is 2.7357 and 43.55 percent is between 2 and 3)

b̂ is close to zero (the average is 0.0773 and 57.14 percent is between 0.01 and 0.1)

l̂ is around 2 (the average is 2.1996 and 67.14 percent is between 1 and 2.5)

M. Mohammadipour, J.E. Boylan / Omega 40 (2012) 703–712 709

where the forecasts are produced by adding up the h-step-ahead forecasts over the forecast horizon.

The simulation results suggest that, in most cases, the aggre- gation method generates forecasts with smaller MSEs than the cumulative h-step-ahead method. The difference is substantial when the autoregressive parameter is high. The only case in which the INARMA-h method is better than the INARMA-Agg method is for an INARMA(1,1) process with small autoregressive and moving average parameters and short length of forecast horizon.

The performance of these forecasts is also tested on empirical data of two real demand data series and the results generally confirm the simulation results.

As previously mentioned, this paper has focused on INARMA processes with Poisson innovations. Other discrete self-decom- posable distributions such as generalized Poisson and negative binomial distributions could be used as marginal distributions. Also, the findings of this paper are based on MSE. Other perfor- mance measures could be used to examine the accuracy of forecasts. In an inventory management context, this can be done by looking at inventory implication metrics such as service level and inventory level [34].

Finally, aggregated forecasts are becoming increasingly impor- tant for Decision Support Systems (DSS) in the area of production planning [35]. Further research into issues related to the applica- tion of aggregated forecasts in such a context should be very important both from academic and practitioner perspectives.

Acknowledgments

The authors would like to thank the three anonymous referees for their comments, which have assisted in sharpening the statements of motivation, scope and contribution of the paper.

Appendix A. Horizon forecasting for an INAR(1) model

In this appendix, it is shown how to derive the conditional second moment of a forecast horizon aggregated PoINAR(1) process ðmZ ¼s

2 Z ¼lÞ. The aggregated process over horizon l can be written as

Xl i ¼ 1

Y tþi ¼ d ða3Y tþZtþ1Þþða23Y tþa3Ztþ1þZtþ2Þþ ���

þðal3Y tþal�13Ztþ1þal�23Ztþ2þ���þa3Ztþ l�1þZtþ lÞ ðA:1Þ

where ¼ d means equal in distribution. It can be simplified as

Xl i ¼ 1

Y tþi ¼ d ða3Y tþa23Y tþ���þal3Y tÞþðZtþ1þa3Ztþ1

þ���þal�13Ztþ1ÞþðZtþ2þa3Ztþ2þ���þal�23Ztþ2Þ þ ���þðZtþ l�1þa3Ztþ l�1ÞþZtþ l ðA:2Þ

We know that

covðai3X,aj3XÞ¼aiajEðX2Þ�aiEðXÞajEðXÞ¼aiajvarðXÞ:

Hence, we have covðai3Zt ,aj3ZtÞ¼aiajl. The variance of (A.2) given Yt is

var Xl i ¼ 1

Y tþi9Y t

" # ¼varða3Y t9Y tÞþvarða23Y t9Y tÞþ �� �þvarðal3Y t9Y tÞ

þ2covða3Y t ,a23Y t9Y tÞþ2covða3Y t ,a33Y t9Y tÞ þ �� �þ2covða3Y t ,al3Y t9Y tÞþ2covða23Y t ,a33Y t9Y tÞ

þ2covða23Y t ,a43Y t9Y tÞþ �� �þ2covða23Y t ,al3Y t9Y tÞ

þ �� �þ2covðal�23Y t ,al�13Y t9Y tÞ þ2covðal�23Y t ,al3Y t9Y tÞþ2covðal�13Y t ,al3Y t9Y tÞ

þvarðZtþ1Þþvarða3Ztþ1Þþvarða23Ztþ1Þþ �� � þvarðal3Ztþ1Þþ2covðZtþ1,a3Ztþ1Þ þ2covðZtþ1,a23Ztþ1Þþ �� �þ2covðZtþ1,al3Ztþ1Þ

þ2covða3Ztþ1,a23Ztþ1Þþ2covða3Ztþ1,a33Ztþ1Þ

þ �� �þ2covða3Ztþ1,al3Ztþ1Þþ �� � þ2covðal�23Ztþ1,al�13Ztþ1Þ þ2covðal�23Ztþ1,al3Ztþ1Þþ2covðal�13Ztþ1,al3Zt

þ1ÞþvarðZtþ2Þþvarða3Ztþ2Þþ �� � þvarðal�13Ztþ2Þþ2covðZtþ2 ,a3Ztþ2Þ þ2covðZtþ2,a23Ztþ2Þþ �� �þ2covðZtþ2,al�13Ztþ2Þ

þ2covða3Ztþ2,a23Ztþ2Þþ2covða3Ztþ2,a33Ztþ2Þ þ �� �þ2covða3Ztþ2,al�13Ztþ2Þþ �� � þ2covðal�33Ztþ2,al�23Ztþ2Þ þ2covðal�33Ztþ2,al�13Ztþ2Þ þ2covðal�23Ztþ2,al�13Ztþ2Þþ �� � þvarðZtþ l�2Þþvarða3Ztþ l�2Þþvarða23Ztþ l�2Þ þ2covðZtþ l�2 ,a3Ztþ l�2Þ þ2covðZtþ l�2,a23Ztþ l�2Þþ2covða3Ztþ l�2,a23Ztþ l�2Þ

þvarðZtþ l�1Þþvarða3Ztþ l�1Þ

M. Mohammadipour, J.E. Boylan / Omega 40 (2012) 703–712710

þ2covðZtþ l�1 ,a3Ztþ l�1ÞþvarðZtþ lÞ ðA:3Þ

Since Yt is fixed, covðai3Y t ,aj3Y t9Y tÞ¼aiajvarðY t9Y tÞ¼0. Hence

var Xl i ¼ 1

Y tþi9Y t

" # ¼að1�aÞEðY t9Y tÞþa2ð1�a2ÞEðY t9Y tÞþ �� �

þalð1�alÞEðY t9Y tÞþlþ½a2lþað1�aÞl� þ½a4lþa2ð1�a2Þl�þ �� �þ½a2llþalð1�alÞl� þ2½aþa2þ���þal�lþ2½a3þa4þ���þalþ1�l

þ���þ2½a2l�3þa2l�2�lþ2½a2l�1�lþl þ½a2lþað1�aÞl�þ½a4lþa2ð1�a2Þl�þ �� � þ½a2l�2lþal�1ð1�al�1Þl�þ2½aþa2þ��� þal�1�lþ2½a3þa4þ���þal�lþ���þ2½a2l�5

þa2l�4�lþ2½a2l�3�lþ���þlþ½a2lþað1�aÞl�

þ½a4lþa2ð1�a2Þl�þ2½aþa2�lþ2½a3�lþl þ½a2lþað1�aÞl�þ2alþl ðA:4Þ

The above result can be summarized to

var Xl i ¼ 1

Y tþi9Y t

" # ¼Y t

Xl j ¼ 1

ajð1�ajÞþ l

1�a l� Xl j ¼ 1

aj 2 4

3 5

þ 2l

1�a Xl j ¼ 1

a2j�1 ðl�jÞ� að1�al�jÞ

1�a

� � ðA:5Þ

Appendix B. Horizon forecasting for an INARMA(p,q) model

In order to find the conditional mean of the forecast horizon aggregated process, we need to express the aggregated INAR- MA(p,q) process in terms of the last p observations (Yt�pþ1, Yt�pþ2,y,Yt�1,Yt). The aggregated process is given by

Xl j ¼ 1

Y tþj ¼Y tþ1þY tþ2þ. . .þY tþ l ðB:1Þ

To write the above equation in the form of (19), we need to know the following:

�

the number and the coefficient of fY t�wþ1g p w ¼ 1

and

�

the number, the coefficient and the subscript of Ztþkij

Table C1 MSEAgg/MSEh of aggregated forecasts for INAR(1) series when l¼6.

Parameters n¼24 n¼36 n¼48 n¼96 n¼500

a¼0:1,l¼0:5 0.9534 0.9750 0.9961 0.9979 0.9534 a¼0:5,l¼0:5

0.9515 0.9677 0.9929 0.9966 0.9515

a¼0:5,l¼1 0.8610 0.8545 0.8474 0.8302 0.8610 a¼0:1,l¼3 0.9444 0.9567 0.9811 0.9923 0.9444 a¼0:5,l¼3 0.7693 0.7365 0.6942 0.6742 0.7693 a¼0:1,l¼5 0.9211 0.9398 0.9677 0.9857 0.9211 a¼0:5,l¼5 0.9534 0.9750 0.9961 0.9979 0.9534

in (B.1). Each of these is discussed in the following subsections.

B.1. The number of fY t�wþ1g p w ¼ 1

Each of the fY tþjg l j ¼ 1

in the RHS of (B.1) needs to be expressed in terms of fY t�wþ1g

p w ¼ 1

by repeated substitution of Ytþj in the equation for the INARMA(p,q) model (11). Because the autore- gressive order of the process is p, Ytþj can be expressed in terms of p previous observations by utilizing the first component of the RHS of (11), namely, a13Y tþj�1þ���þap3Y tþ j�p.

Now, if jrp�(w�1), there is one Yt�(w�1) when we express the jth observation in the RHS of (B.1) (Ytþj) without any need for further substitution. Repeated substitution of (Ytþ1,y, Ytþp�(w�1)) by their p previous observations would result in obtaining more Yt�(w�1). Therefore, in total, the number of Yt�(w�1) in each of fY tþ jg

lþ1 j ¼ 1

when jrp�(w�1) is equal to the number of Yt�(w�1) in its p previous observations plus one.

However, when j4p�(w�1), each Ytþj from (B.1) should be substituted by (11) in order to reach Yt�(w�1), and the number of

Yt�(w�1) in each of fY tþ jg l j ¼ 1

would be equal to the number of Yt�(w�1) in its p previous observations.

B.2. The coefficient of fY t�wþ1g p w ¼ 1

For jrp�(w�1), the corresponding coefficient of Yt�(w�1) in the jth observation in the RHS of (B.1), Ytþj, is ajþ(w�1) because

Y tþ j ¼a13Y tþj�1þ���þajþðw�1Þ3Y t�ðw�1Þþ �� �þap3Y tþj�pþZtþ j

þ Xq i ¼ 1

bi3Ztþ j�i

For other Yt�(w�1) the coefficient in each of fY tþjg l j ¼ 1

is ai thinned the coefficient of Yt�(w�1) in the ith previous observation for i¼1,y,p.

For j4p�(w�1), again, the coefficient of Yt�(w�1) in each of fY tþ jg

l j ¼ 1

is ai thinned the coefficient of Yt�(w�1) in the ith previous observation for i¼1,y,p (the difference with the pre- vious case is that we do not have ajþ(w�1)).

B.3. The number of Ztþkij

Now we come back to (B.1) to find the Z terms in each of fY tþ jg

l j ¼ 1

in the RHS of the equation when they are expressed in terms of fY t�wþ1g

p w ¼ 1

. As the process has a moving average component of order q, each fY tþ jg

l j ¼ 1

has qþ1 innovation terms {Ztþj,Ztþj�1,y,Ztþj�q}. However, by repeated substitution, each fY tþ jg

l j ¼ 1

can be expressed in terms of p previous observations, each also with qþ1 innovation terms.

Therefore, the total number of innovation terms in each of fY tþ jg

l j ¼ 1

is equal to the number of innovation terms in the p previous observations, plus qþ1.

B.4. The coefficient of Ztþkij

The corresponding coefficients for the qþ1 terms {Ztþj,Ztþj�1,y,Ztþj�q} are {1,b1,y,bq}. For the innovation terms that come from the p previous observations, coefficients would be ak thinned the coefficient of Ztþkij in the kth previous observation for k¼1,y,p.

B.5. The subscript of Ztþkij

tþkij denotes the subscript of Z for each i,j (j¼1,y,l and i¼1,.. .,n

pþ1 j

). Each fY tþjg l j ¼ 1

has qþ1 innovation terms {Ztþj, Ztþj�1,y,Ztþj�q}. Therefore, the subscripts for the last qþ1 innova- tion terms in each fY tþ jg

l j ¼ 1

are {j�q,j�1,y,j}. This is shown in Table 1 by i¼n

pþ1 j�1 þ1,.. .,n

pþ1 j�1 þn

pþ1 j

. The other subscripts of innovation terms in each of fY tþjg

l j ¼ 1

simply are the subscripts of the innovation terms of p previous observations.

As a result, the aggregated process can be expressed as (18) with the associated parameters as defined in Table 1.

Table C6 MSEAgg/MSEh of aggregated forecasts for INARMA(1,1) series when l¼9.

Parameters n¼24 n¼36 n¼48 n¼96 n¼500

a¼0:1,b¼0:1,l¼0:5 0.8270 1.1104 1.1042 1.0444 0.9984

M. Mohammadipour, J.E. Boylan / Omega 40 (2012) 703–712 711

Appendix C. Simulation results for horizon-aggregated forecasts for INAR(1), INMA(1) and INARMA(1,1) models

In this appendix, the forecast accuracy of the two forecasting methods (INARMA-Agg and INARMA-h) is compared in terms of

Table C2 MSEAgg/MSEh of aggregated forecasts for INAR(1) series when l¼9.

Parameters n¼24 n¼36 n¼48 n¼96 n¼500

a¼0:1,l¼0:5 0.6497 0.8698 0.9284 0.9868 0.9983 a¼0:5,l¼0:5 0.6973 0.8425 0.8887 0.9251 0.9189 a¼0:1,l¼1 0.6418 0.8723 0.9305 0.9829 0.9960 a¼0:5,l¼1 0.6924 0.8362 0.9046 0.9056 0.8920 a¼0:1,l¼3 0.6657 0.8574 0.8979 0.9717 0.9939 a¼0:5,l¼3 0.6235 0.7602 0.8016 0.8048 0.7732 a¼0:1,l¼5 0.6055 0.8669 0.9119 0.9677 0.9904 a¼0:5,l¼5 0.6384 0.7262 0.7532 0.7052 0.6736

Table C3 MSEAgg/MSEh of aggregated forecasts for INMA(1) series when l¼6.

Parameters n¼24 n¼36 n¼48 n¼96 n¼500

b¼0:1,l¼0:5 0.9996 0.9996 0.9999 1.0000 1.0000 b¼0:5,l¼0:5 0.9982 0.9989 0.9992 0.9996 0.9999 b¼0:9,l¼0:5 0.9969 0.9978 0.9984 0.9992 0.9998 b¼0:1,l¼1 0.9997 0.9998 0.9999 1.0000 1.0000 b¼0:5,l¼1 0.9986 0.9988 0.9992 0.9995 0.9999 b¼0:9,l¼1 0.9966 0.9979 0.9983 0.9991 0.9998 b¼0:1,l¼3 0.9999 0.9996 1.0000 1.0000 1.0000 b¼0:5,l¼3 0.9974 0.9987 0.9992 0.9995 0.9999 b¼0:9,l¼3 0.9948 0.9975 0.9981 0.9991 0.9998 b¼0:1,l¼5 1.0002 0.9999 0.9998 1.0000 1.0000 b¼0:5,l¼5 0.9976 0.9988 0.9992 0.9995 0.9999 b¼0:9,l¼5 0.9945 0.9975 0.9980 0.9990 0.9998

Table C4 MSEAgg/MSEh of aggregated forecasts for INMA(1) series when l¼9.

Parameters n¼24 n¼36 n¼48 n¼96 n¼500

b¼0:1,l¼0:5 0.9995 0.9998 1.0000 1.0000 1.0000 b¼0:5,l¼0:5 0.9984 0.9992 0.9995 0.9997 0.9999 b¼0:9,l¼0:5 0.9973 0.9986 0.9989 0.9995 0.9999 b¼0:1,l¼1 0.9995 0.9999 0.9999 1.0000 1.0000 b¼0:5,l¼1 0.9984 0.9993 0.9995 0.9997 0.9999 b¼0:9,l¼1 0.9973 0.9985 0.9989 0.9994 0.9999 b¼0:1,l¼3 0.9991 0.9998 0.9999 1.0000 1.0000 b¼0:5,l¼3 0.9978 0.9994 0.9994 0.9997 0.9999 b¼0:9,l¼3 0.9973 0.9985 0.9990 0.9994 0.9999 b¼0:1,l¼5 1.0000 1.0000 0.9999 1.0000 1.0000 b¼0:5,l¼5 0.9975 0.9987 0.9995 0.9997 0.9999 b¼0:9,l¼5 0.9955 0.9976 0.9987 0.9994 0.9999

Table C5 MSEAgg/MSEh of aggregated forecasts for INARMA(1,1) series when l¼6.

Parameters n¼24 n¼36 n¼48 n¼96 n¼500

a¼0:1,b¼0:1,l¼0:5 1.1280 1.1290 1.1129 1.0340 0.9975 a¼0:1,b¼0:9,l¼0:5 0.9817 0.9902 0.9727 0.9676 0.9864 a¼0:5,b¼0:5,l¼0:5 0.9291 0.8821 0.9278 0.9009 0.9044 a¼0:1,b¼0:1,l¼1 1.0289 1.1009 1.0355 1.0200 0.9991 a¼0:1,b¼0:9,l¼1 0.9652 1.0028 0.9657 0.9736 0.9857 a¼0:5,b¼0:5,l¼1 0.8933 0.9272 0.8895 0.9066 0.9065 a¼0:1,b¼0:1,l¼5 0.8930 1.0278 1.0238 1.0005 0.9966 a¼0:1,b¼0:9,l¼5 0.9320 0.9472 0.9848 0.9668 0.9829 a¼0:5,b¼0:5,l¼5 0.8281 0.8787 0.8837 0.8904 0.9029

a¼0:1,b¼0:9,l¼0:5 0.7848 0.9304 0.9364 0.9652 0.9898 a¼0:5,b¼0:5,l¼0:5 0.8773 0.8642 0.9091 0.9123 0.9389 a¼0:1,b¼0:1,l¼1 0.7756 0.9952 1.0110 1.0206 0.9981 a¼0:1,b¼0:9,l¼1 0.6699 0.9117 0.9455 0.9603 0.9901 a¼0:5,b¼0:5,l¼1 0.5531 0.8320 0.8976 0.9279 0.9372 a¼0:1,b¼0:1,l¼5 0.3814 0.8981 0.9941 1.0063 0.9967 a¼0:1,b¼0:9,l¼5 0.3441 0.8315 0.9185 0.9654 0.9862 a¼0:5,b¼0:5,l¼5 0.3577 0.8164 0.8840 0.9123 0.9332

MSE. The results are for the cases where data is produced by an INAR(1), an INMA(1) or an INARMA(1,1) process and the forecast horizon is l¼6 or l¼9.

See Tables C1–C6.

References

[1] Quenouille MH. Discrete autoregressive schemes with varying time intervals. Metrika 1958;1:21–7.

[2] Amemiya T, Wu. RY. The effect of aggregation on prediction in the autoregressive model. Journal of the American Statistical Association 1972;67(339): 628–32.

[3] Brewer KRW. Some consequences of temporal aggregation and systematic sampling for ARMA and ARMAX models. Journal of Econometrics 1973;1(2): 133–54.

[4] Harvey AC, Pierse RG. Estimating missing observations in economic time series. Journal of the American Statistical Association 1984;79(385):125–31.

[5] Nijman TE, Palm FC. Predictive accuracy gain from disaggregate sampling in ARIMA models. Journal of Business & Economic Statistics 1990;8(4):405–15.

[6] Drost FC, Nijman TE. Temporal aggregation of GARCH processes. Econome- trica 1993;61(4):909–27.

[7] Marcellino M. Some consequences of temporal aggregation in empirical analysis. Journal of Business & Economic Statistics 1999;17(1):129–36.

[8] Teles P, Wei. WWS. The use of aggregate time series in testing for Gaussianity. Journal of Time Series Analysis 2002;23(1):95–116.

[9] Man KS. Linear prediction of temporal aggregates under model misspecifica- tion. International Journal of Forecasting 2004;20(4):659–70.

[10] Tang D. Managing finished-goods inventory under capacitated delayed differentiation. OMEGA: The International Journal of Management Science 2011;39(5):481–92.

[11] Hosoda T, Disney SM. On variance amplification in a three-echelon supply chain with minimum mean square error forecasting. OMEGA: The Interna- tional Journal of Management Science 2006;34(4):344–58.

[12] Thury G, Witt SF. Forecasting industrial production using structural time series models. OMEGA: The International Journal of Management Science 1998;26(6):751–67.

[13] Brännäs K, Hellström J, Nordström. J. A new approach to modelling and forecasting monthly guest nights in hotels. International Journal of Forecast- ing 2002;18(1):19–30.

[14] Freeland RK, McCabe BPM. Forecasting discrete valued low count time series. International Journal of Forecasting 2004;20:427–34.

[15] Gourieroux C, Jasiak J. Heterogeneous INAR(1) model with application to car insurance. Insurance: Mathematics and Economics 2004;34(2):177–92.

[16] Brännäs K, Quoreshi. S. Integer-valued moving average modelling of the number of transactions in stocks. Applied Financial Economics 2010;20(18): 1429–40.

[17] Sueutel FW, van Harn K. Discrete analogues of self-decomposability and stability. Annals of Probability 1979;7(5):893–9.

[18] Mohammadipour M. Intermittent demand forecasting with integer autore- gressive moving average models. Buckinghamshire New University, Brunel University; 2009.

[19] Du JG, Li. Y. The integer-valued autoregressive (INAR(p)) model. Journal of Time Series Analysis 1991;12(2):129–42.

[20] McKenzie E. Some ARMA models for dependent sequences of Poisson counts. Advances in Applied Probability 1988;20(4):822–35.

[21] Neal P, Rao TS. MCMC for integer-valued ARMA processes. Journal of Time Series Analysis 2007;28(1):92–110.

[22] Brännäs K, Hellström J. Generalized integer-valued autoregression. Econo- metric Reviews 2001;20(4):425–43.

[23] Jung RC, Tremayne AR. Coherent forecasting in integer time series models. International Journal of Forecasting 2006;22:223–38.

[24] Eaves AHC. Forecasting for the ordering and stock-holding of consumable spare parts. Lancaster University; 2002.

[25] Al-Osh MA, Alzaid AA. First order integer-valued autoregressive (INAR(1)) processes. Journal of Time Series Analysis 1987;8(3):261–75.

M. Mohammadipour, J.E. Boylan / Omega 40 (2012) 703–712712

[26] Bu R, McCabe B, Hadri K. Maximum likelihood estimation of higher-order integer-valued autoregressive processes. Journal of Time Series Analysis 2008;29(6):973–94.

[27] Silva ME, Oliveira VL. Difference equations for the higher-order moments and cumulants of the INAR(1) model. Journal of Time Series Analysis 2004;25(3): 317–33.

[28] Freeland RK, McCabe BPM. Asymptotic properties of CLS estimators in the Poisson AR(1) model. Statistics and Probability Letters 2005;73(2):147–53.

[29] Cox DR. Prediction by exponentially weighted moving averages and related

methods. Journal of the Royal Statistical Society: Series B 1961;23:414–22.

[30] Tiao GC, Xu D. Robustness of maximum likelihood estimates for multi-step

predictors: the exponential smoothing case. Biometrika 1993;80:623–41.

[31] Kang I-B. Multi-period forecasting using different models for different horizons: an application to U.S. economic time series data. International Journal of Forecasting 2003;19:387–400.

[32] Syntetos AA, Boylan. JE. The accuracy of intermittent demand estimates. International Journal of Forecasting 2005;21:303–14.

[33] Latour A. Existence and stochastic structure of a non-negative integer-valued autoregressive process. Journal of Time Series Analysis 1998;19(4):439–55.

[34] Teunter RH, Duncan L. Forecasting intermittent demand: a comparative study. Journal of the Operational Research Society 2009;60(3):321–9.

[35] da Silva CG, et al. An interactive decision support system for an aggregate production planning model based on multiple criteria mixed integer linear programming. OMEGA: The International Journal of Management Science 2006;34(2):167–77.