Let z_t be a time series process. For a stationary process, its mean, E(z_t) = μ, and variance, γz(0)=E(zt-μ)2=σz2, are constant. Also, in this case, its autocovariance function (ACF) between z_t and z_t+k, γ_z(k) = E(z_t − μ)(z_t+k − μ) = E(ż_tż_t+k), and autocorrelation function, ρ_z(k) = γ_z(k)/γ_z(0), are functions of only the time difference. The partial autocorrelation function (PACF) is defined as φ_kk = Corr(z_t, z_t+k|z_t+1, . . . , z_t+k−1). Some commonly used time series processes or models are:

• Autoregressive process of order p (AR(p) model)

  (1.1)z˙t-φ1z˙t-1-⋯-φpz˙t-p=at,(1-φ1B-⋯-φpBp)z˙t=at,φp(B)z˙t=at.

• Moving average process of order q (MA(q) model)

  (1.2)z˙t=at-θ1at-1-θ2at-2-⋯-θqat-p=(1-θ1B-θ1B2-⋯-θqBq) at=θq(B)at.

• Autoregressive moving average process (ARMA(p, q) model)

  (1.3)φp(B)z˙t=θq(B)at.

• Autoregressive integrated moving average process (ARIMA(p, d, q) model)

  (1.4)φp(B)(1-B)dzt=θ0+θq(B)at.

The model is stationary if the roots of its associated AR polynomial are all outside the unit circle, and the important characteristics of stationary models can be summarized in the following table:

Time series are used in many studies either for model building or inference. We must be careful when choosing what kind of time series data is used in the analysis. Since many time series variables like rainfall, industrial production, and sales exist only in some aggregated forms, we will begin with the issue related to the temporal aggregation effect on the model form. Given a time series z_t, let Z_T = (1 + B + · · · + B^m−1)z_mT . For example, with m = 3, Z₁ = (1 + B + B²)z₃ = z₁ + z₂ + z₃, Z₂ = (1 + B + B²)z₆ = z₄ + z₅ + z₆, etc. We will call z_t as non-aggregate series and Z_T as aggregate series. For m = 3 if z_t is a monthly series, then Z_T will be a quarterly series. To make inference, should we use non-aggregate series, z_t or aggregate series, Z_T? Do they make any difference? Are the time series models for z_t and for Z_T the same?

The first published papers on aggregation effects on ARIMA models were by Tiao (1972) and Amemiya and Wu (1972), and they led to many other studies on the topic including my Ph.D. dissertation and life time research in the area.

To answer the above questions, we need to study the relationship of autocovariances between the non-aggregate and aggregate series z_t and Z_T, or more generally between w_t = (1 − B)^dz_t and U_T = (1 − ℬ)^dZ_T . Define the m-period overlapping sum,

and note that Z_T = ζ_mT, (1 − ℬ)Z_T = Z_T − Z_T−1 = ζ_mT − ζ_m(T−1) = (1 − B^m)ζ_mT,

Hence, we have

Now let us consider an MA(2) model for z_t, z_t = (1 − θ₁B − θ₂B²)a_t. If m = 3, what is the model for Z_T? For this MA(2) model, we have

and

Note that for d = 0 and m = 3, from (2.1), we have

where

Hence, Z_T is a MA(1) process, Z_T = (1 − Θℬ)A_T where

and

More generally, we have the following results from Stram and Wei (1986):

• Temporal aggregation of the AR(p) process

  Suppose that the non-aggregate series z_t follows a stationary AR(p) process,

  (1-φ1B-⋯-φpBp) zt=at.

  Let φ_p(B) = (1 − φ₁B − · · · − φ_pB^p)z_t and δi-1 for i = 1, . . . , p^* be the distinct roots of φ_p(B), each with multiplicity s_i such that ∑i=1p*si=p. For any given value of m, let b equal the number of distinct values δim for i = 1, . . . , p^*. Furthermore, partition the numbers s_i for i = 1, . . . , p^* into b distinct sets A_i such that s_k and s_j ∈ A_i if and only if δkm=δjm. Then the mth order aggregate series, Z_T follows an ARMA(M, N₁) model,

  (1-α1B-⋯-αMBM) ZT=(1-β1B-⋯-βN1BN1) ET,

  where M=∑i=1bmax Ai, max A_i = the largest element in A_i, N₁ = [p+1−(p+1)/m]−(p−M) = [M + 1 − (p + 1)/m], the E_T are white noise with mean 0 and variance σE2, and α_i, β_j, and σE2 are functions of φ_k’s and σa2.

• If z_t ~ ARMA(p, q) model, then Z_T ~ ARMA(M, N₂), where N₂ = [p+1+(q−p−1)/m]−(p−M).

• If z_t ~ ARIMA(p, d, q) model, then Z_T ~ ARIMA(M, d, N₃), where N₃ = [p + d + 1 + (q − p − d − 1)/m] − (p − M).

The limiting behavior of aggregates was studied by Tiao (1972) and he showed that given z_t ~ ARIMA(p, d, q) model, the limiting model for the aggregates, Z_T exists, and as m → ∞, Z_T → IMA(d, d).

When a variable is a stock variable and we observe only every mth value of the variable, i.e., given z₁, z₂, z₃, . . . but we observe only Z_T = z_mT. For example, for m = 3, Z₁ = z₃, Z₂ = z₆, . . . . We have the following interesting result from Wei (1981).

Given z_t ~ IMA(d, q): (1 − B)^dz_t = (1 − θ₁B − θ₁B² − · · · − θ_qB^q)a_t. Let Z_T = z_mT. Then, as m → ∞, Z_T → IMA(d, d − 1). Theoretically, every ARIMA(p, d, q) process can be approximated by an IMA(d, q) process. Thus, if z_t is an ARIMA(p, 1, q), then as m → ∞, Z_T approaches IMA(1, 0), which could very likely explain why most daily stock prices follow a random walk model. In building an underlying time series model, we need to be aware of the effect of the use of aggregate series. In addition to the above cited references, we refer readers to some other useful references including Brewer (1973), Wei (1982), Ansley and Kohn (1983), Weiss (1984), Wei and Stram (1988), Marcellino (1999), Shellman (2004), and Sbrana and Silvestrini (2013).

Given the AR(1) model

where a_t is N(0,σa2) white noise process and t = 1, 2, . . . , n. To test a unit root, H₀: φ = 1 vs H₁: φ < 1, since

Dickey and Fuller (1979) suggested using the following test statistic and showed that under H₀: φ = 1,

where W(t) is a Wiener process (also known as Brownian motion process). We reject H₀ if the value of the test statistic is too small (negative).

In practice, aggregate data, Z_T = (1 + B + · · · + B^m−1)z_mT are often used. It has been shown by Teles et al. (2008) that under H₀: φ = 1, (1 − B)z_t = a_t, the corresponding model for the aggregate series is

where the E_T−1 are independent and identically distributed variables with zero mean and variance σE2 and the parameters Θ and σE2 are determined as follows:

• If m = 1, then Θ = 0; σE2=σa2;

• If m ≥ 2, then

  Θ=-2m2+1m2-1+{(2m2+1m2-1)2-1}12,σE2=σa2m (2m2+1)3 (1+Θ2).

  As a result, the corresponding test statistic becomes

  (3.4)N(φ^-1)=N-1∑T=2NZT-1ETN-2∑T=2NZT-12-ΘN-1∑T=2NZT-1ET-1N-2∑T=2NZT-12→D12{[W(1)]2-1}∫01[W(x)]2dx+m2-16m2∫01[W(x)]2dx.

  Comparing with (3.2), we see that the limiting distribution of the test statistic for the aggregate time series depends on the order of aggregation m. Since m ≥ 2, the distribution of the test statistic is shifted to the right, and the shift increases with the order of aggregation. Aggregation leads to empirical significance levels lower than the nominal level and significantly reduces the power of the test.

Example 1

In this example, we simulated a time series of 240 observations from the model z_t = 0.95z_t−1 + a_t where the a_t are i.i.d. N(0, 1). This series was then aggregated with m = 3.

• Test a unit root based on non-aggregate series, z_t (240 observations). Based on the sample autocorrelation function and partial autocorrelation function, we have an AR(1) model. The least squares estimation leads to the following result

  z^t=0.9603(.0182) zt-1.

  To test the hypothesis of a unit root, the value of the test statistic is

  n(φ^-1)=240(0.9603-1)=-9.53.

  At α = 5%, the critical point from Dickey and Fuller (1979) is between −8.0 and −7.9 for n = 240. Thus, the hypothesis of a unit root is rejected, and we conclude that the underlying model is stationary. This is consistent with the underlying simulated model.

• Test a unit root with aggregate series, Z_T, with m = 3 (80 observations). The sample autocorrelation function and partial autocorrelation function suggest an AR(1) model (an ARMA(1,1) model was also considered but its MA parameter was not significant). The least squares estimation leads to the following result

  Z^t=0.9205(.0448) zt-1.

  To test the hypothesis of a unit root, the value of the test statistic is

  n(φ^-1)=80(0.9205-1)=-6.36.

  Again, at α = 5%, the critical point from Dickey and Fuller (1979) is between −7.9 and −7.7 for n = 80. Thus, the hypothesis of a unit root is not rejected, and we conclude that the underlying model is nonstationary. This leads to a wrong conclusion. However, if we use the adjusted critical value given in Teles et al. (2008) based on the adjusted test statistic of (3.4), which is between −5.45 and −5.40 for n = 80, we will reject the null hypothesis of a unit root and leads to a consistent conclusion.

When aggregate time series are used in modeling and testing, we need to make sure to use a proper adjusted table for the test of its significance.

Time series are often used in regression analysis, which is possible the most commonly used statistical method. So we will also consider the consequence of the use of aggregate series in a regression model. Let us consider the simple regression model,

which is a one-sided causal relationship. If x_t−1 is also stochastic, for example, if it follows a MA(1) process, we can also write the joint system as

or

where a_t and e_t are independent N(0,σa2) and N(0,σe2), respectively. Let YT=(∑j=0m-1Bj)ymT,XT=(∑j=0m-1Bj)xmT,WT=(∑j=0m-1Bj)xmT-1, and ET=(∑j=0m-1Bj)emT. Equation (4.1) implies that

For m = 3, W₁ = x₀ + x₁ + x₂, W₂ = x₃ + x₄ + x₅, etc., which are not available, and the available data are X₁ = x₁ + x₂ + x₃, X₂ = x₄ + x₅ + x₆, etc. A natural way to estimate W_T is to consider its projection on X_T. Specifically, we let Z = [W_T, X_T]′ and compute its covariance matrix generating function

with

and hence

and

It follows that

It is interesting to note that Equation (4.6) can be rewritten as

which implies that the estimate of Ŵ_T is the weighted average of X_T and X_T−1 with weights (m−1)/m and 1/m, respectively. This is clearly reasonable. The aggregate model then becomes

where U_T = α(W_T − Ŵ_T ) + E_T = αV_T + E_T, GU(B)=α2GV(B)+mσe2, GV(B)=G11(B)-G12(B)[G22(B)]-1G21(B)=σa2{[m(1-θ)2]+2θ-θ(B+F)-[(1-θ)4[(m-1)+B][(m-1)+F]]/[m(1-θ)2+2θ-θ(B+F)]}, and F = B⁻¹. Thus, temporal aggregation turns a one-sided causal relationship into a two-sided feedback system. It is important to note that after proposing an underlying model for a study, one should use the same time unit in the hypothesis and data collection for modeling and testing. An improper use of time unit could lead to a very misleading conclusion. For a more detailed description, we refer readers to Tiao and Wei (1976). Other useful references include Wei (1978), and Lütkepohl (1987).

5.1. Representation of vector time series models

In studying the relationship of variables, other than the regression model, we often consider vector time series models. Although the basic procedures of model building between univariate time series and vector time series are the same, there are some important phenomena which are unique to vector time series models. We now discuss some special issues of vector time models.

First, let us review some results from univariate time series models. It is well known that we can always write a stationary process as a MA representation

or

such that ∑j=0∞∣ψj∣<∞. Similarly, we can write an invertible process as an AR representation

or

such that ∑j=0∞∣πj∣<∞. From these two representations, we have the well-known dual relationship between AR(p) and MA(q) models in the univariate time series processes. That is, a finite order AR process corresponds to an infinite order MA process, and a finite order MA process corresponds to an infinite order AR process. For example, an AR(1) model, (1 − φB)Ż_t = a_t, corresponds to an infinite order MA process, Ż_t = (1 − φB)⁻¹a_t = (1 + φB + φ²B² + · · · )a_t, and a MA(1) model, Ż_t = (1 − θB)a_t, corresponds to an infinite order AR process, (1 + θB + θ²B² + · · · )Ż_t = a_t.

Let Z_t = [Z_1,t, Z_2,t, . . . , Z_m,t]′ be the m-dimensional vector time series. Some commonly used vector time series models are VAR(p), VMA(q), and VARMA(p, q) processes. Again, we can write these vector processes in a moving average representation

where Ż = Z − μ and a_t is a m-dimensional vector white noise process N(0, ∑). We can also write it in an autoregressive representation

It follows that we can express a VAR(1) process,

in the following MA representation

where Φ₀ = I. A natural question to ask: is the VMA representation always an infinite order? People often think that the univariate model is a special case of the vector model with dimension equal 1, and so the answer to the question is obviously yes. However, let us consider the following 3-dimensional vector VAR(1) model

or equivalently

where Φ=[0.8000.5000]. Since Φ² ≠ 0 and Φ^j = 0 for j > 2, Equation (5.4) actually represents a VMA(2) model. In fact,

Thus, the inverse of a non-degenerate VAR(1) matrix polynomial (i.e., Φ(B) ≠ I) will be of a finite order if the determinant |Φ(B)| is independent of B. For more detailed discussion, we refer readers to Tiao and Tsay (1989), and Shen and Wei (1995).

5.2. Representation of multiplicative seasonal vector autoregressive moving average models

Given a univariate seasonal ARMA model

where

and a_t is a Gaussian white noise process with mean 0 and a constant variance σa2. The model is often denoted as ARMA(p, q)×(P, Q)_s. To facilitate our discussion, we will use the order of the polynomials appearing in the equation and denote it as ARMA(P)_s(p)(q)(Q)_s. When x_t = [x_1,t, . . . , x_k,t]′ is a k-dimensional vector, the natural extension is the following multiplicative vector autoregressive moving average VARMA(P)_s(p)(q)(Q)_s model,

where

are matrix polynomials. The matrix I is the k-dimensional identity matrix, the Φs, φs, θs, and Θs are k × k parameter matrices, and a_t is a vector Gaussian white noise process with mean vector 0 and E(atat′)=Ω.

Note that as expected, the vector model reduces to the univariate model when k = 1. Moreover, in such a case, the ARMA(P)_s(p)(q)(Q)_s model can also be written as the following ARMA(p)(P)_s(Q)_s(q) model

As a result, a multiplicative seasonal ARMA model is also traditionally written in the form of ARMA (P)_s(p)(q)(Q)_s. This traditional representation has been adopted by many researchers for both univariate and vector time series.

When k > 1 and x_t is a vector process, can we really extend the above operation and write the VARMA(P)_s(p)(q)(Q)_s model, Φ_P(B^s)φ_p(B)x_t = θ_q(B)Θ_Q(B^s)a_t, as the following VARMA(p)(P)_s(Q)_s (q) model?

Let us consider two simple seemingly equivalent bivariate VAR(1)₄(1) and VAR(1)(1)₄ representations with the following parameters

and the associated noise a_t is a vector Gaussian white noise process with mean zero and covariance matrix Ω. For the VAR(1)₄(1) representation,

For the VAR(1)(1)₄ representation,

The different implications between the two seemingly equivalent representations are clear and cannot be ignored especially when a policy related decision is to be made. Please see Yozgatligil and Wei (2009) for details.

In addition to temporal aggregation discussed in earlier sections, there is another commonly used aggregation. For example, the total money supply is the aggregate of demand deposits and currency in circulation. The total housing start is the aggregate of housing starts in the north east, north central, south, and west regions, which again are the subaggregates of housing starts in different states. The total sales of a company is the aggregate of the sales achieved by all of its branches throughout the country or countries.

Let z_1,t, z_2,t, . . . , z_m,t be the m component time series and Yt=∑i=1mzi,t be the corresponding series of aggregates. Suppose we are interested in forecasting the future aggregate Y_t+l for some l, based on the knowledge of the available time series up to the time t. Clearly, such forecasts can be obtained through the following three methods.

• Method 1: based on a model using the aggregate series Y_t and its l-step ahead forecast, Ŷ_t(l).

• Method 2: based on individual component models using the non-aggregate series z_i,t and the sum of the forecasts from all component models, i.e., Y^^t(l)=∑i=1mz^i,t(l).

• Method 3: based on a joint multiple time series model and the forecast from the joint multiple model, Y^^^t(l).

Question: what are the relative efficiencies among the three methods in terms of the minimum mean square error forecast? The answers are given below.

• E[Yt+l-Y^^^t(l)]2≤E[Yt+l-Y^t(l)]2; E[Yt+l-Y^^^t(l)]2≤E[Yt+l-Y^^t(l)]2; and equality holds when z_1,t, z_2,t, . . . , z_m,t are orthogonal to each other.

• The comparison between methods 1 and 2 depends on the model structure; there is no definite winner between methods 1 and 2.

The answer (2) above could be surprising to some people because aggregation normally will cause information loss. For the proof, we refer readers to Wei and Abraham (1981).

Next, let us consider the m-dimensional models related to both time and space where we write it as

The commonly used model to describe (6.1) is the following space-time autoregressive moving average (STARMA(p, q)) model,

where W(l)=⌊wi,j(l)⌋ are m × m the spatial weight matrices, ∑j=1mwi,j(l)=1,