- Email: [email protected]
Testing causality using efficiently parametrized vector ARMA models
Testing Causality Using Efficiently Parametrized Vector ARMA Models Paul Newbold and Steven M. Hotopp Department of Economics University of Illinois, Urbana-Champaign Champaign, Illinois 61820
ABSTRACT We outline a practical methodology for fitting to data vector autoregressive moving-average models, and employ these models to test for Granger causality. The procedures are applied to data on employment and hourly wages in manufacturing in
the United States.
1.
INTRODUCTION Given observations through time on m series X{ = (Xi,, X,,, . . . , Xm(), we
wish
to
fit
a
member
of
the
vector
autoregressive
moving-average
[ARMA( p, q)] class of models
where the api and Oj are m X m matrices of fixed coefficients, vector white noise, so that
and et is
The model (1.1) is an obvious extension to the multivariate case of the univariate models discussed by Box and Jenkins [8]. For the single-series case, Box and Jenkins outlined a model-building methodology involving an iterative cycle of model selection, parameter estimation, and model checking. At the first stage, based on judgmental interpretation of statistics readily computed APPLIED
MATHEMATZCS AND COMPVTATZON 20: 329-348
(1986)
329
0 Elsevier Science Publishing Co., Inc., 1986 52 Vanderbilt Ave., New York, NY 10017
009fUOO3/86/$03.50
PAUL NEWBOLD AND STEVEN M. HOTOPP
330
from the data, a specific model from the general class is chosen for further analysis. Next, using statistically efficient methods, the fixed parameters of the selected model are estimated. Finally, diagnostic checks on the adequacy of representation of the model to the available data are applied. Any inadequacies revealed at this stage may suggest an alternative structure, and the cycle of selection, estimation, and checking is iterated until a satisfactory model is achieved. Some early work on the fitting to data of models of the form (1.1) is reported by Quenouille [43]. More recently, Jenkins and Alavi [33] and Tiao and Box [49] have discussed the extension of the Box-Jenkins approach to this problem. In similar spirit, we report in Section 2 of this paper on a model-building methodology that we have successfully applied to several small sets of economic time series. These procedures will be illustrated by reference to data on employment and hourly wages in manufacturing in the United States. Experience with other data sets is discussed in [32]. In Section 3 we discuss the use of fitted vector autoregressivemoving average models to testing for “causality,” in the sense of Granger [18]. Such tests are likely to have higher power than generally applied procedures based on heavily parametrized structures. 2.
MODEL
BUILDING
In this section we consider the three stages of the model-building cycle, when the objective is to fit to data a model from the class (1.1). Model Selection In practice, even for the case of a single time series, it is the initial selection of an appropriate model from the general class that causes most difficulty in the application of the Box-Jenkins methodology. This difficulty is considerably more acute in the multivariate case. For the analysis of a single series, Box and Jenkins base model selection on examination of the sample autocorrelations and partial autocorrelations of the series. Analogous statistics are also helpful in the multivariate case. Define the population autocovariance matrices of X, as
r,c=E(X,X;+,) Then it is well known that, for a pure moving-average process of order q [p = 0 in (l.l)], rk = 0
forall
k>q.
Testing Causality
Using Vector ARMA Models
331
This suggests that the corresponding sample autocovariance matrices
should be useful in identifying pure moving-average behavior. It is easier to interpret the corresponding autocorrelation matrices R,, whose (i, j) elements are given by
where ck; i, jis the (i, j) element of Ck. For checking for pure moving average behavior of order 4, since the sample autocorrelations have asymptotic normal distributions, comparison can be made with the asymptotic standard errors, given by Bartlett [S] as l/2
4
SE[Tk;i,j]
=
12-
1+2
1 11 =
Ph;i,iPh;j,j
1
I
1
where, in practice, the population autocorrelations in this expression are replaced by the corresponding sample quantities. As in the case of a single time series, pure autoregressive models [q = 0 in (l.l)] can be detected through the sample partial autocorrelations. Consider theregressionofX,onX,_j(j=1,2,...,k): x,
=
@klXf-l+
@k2Xt-2
+
”
’ +
@kkXt-k
+Uk,t
(2.1)
where the error term uk t has covariance matrix Qk, and also the regression of X, on Xt+i (j=l,2,...:k)
where uz t has covariance matrix Qt. I_& A, and A$ denote the symmetric square roots of the matrices at, and Dz. Then, Quenouille [43, p. 411 defines the multiple partial autocorrelation of order k + 1 as P k+l=
&‘E(uk
7tUf’t-k-l)A:-l. 1
332
PAUL NEWBOLD AND STEVEN M. HOTOPP
If the generating process is autoregressive of order k, then u k,t and uf, t k_ 1 will be uncorrelated, and hence P k+l equal to zero. Ledolter [35] and Ansley and Newbold [4] have extended the algorithm of Durbin [13] to the computation of estimates i),, i of the partial autocorrelation matrices. Hannan [24, p. 3981 showsnthat, for an autoregressive process of order k, the elements of the matrices Pk+ j (j > 1) have independent asymptotic normal distributions, with zero means and variances n- ‘. Thus, for such a process, the statistics =
(Ik+j
ntrbk+ji);+j
have asymptotic chi-square distributions with m2 degrees of freedom. Altematively, these statistics can be expressed as (2.3)
a k+jzfltr’k+j,k+j’k*+j,k+j9
where the &‘s are the estimates of the parameter matrices in (2.1) and (2.2). The quantities (2.3) can therefore be employed in the detection of pure autoregressive models. Much more difficult is the distinction, based on statistics conveniently computed from the data, among mixed ARMA(p, o) models. To see why this is so, consider the natural extension of the Box-Jenkins single series approach to model selection to the vector case. Postmultiplying through (1.1) by Xi_, and taking expectations yields
r, = l-k_@;
+ .
’ ’
+
r,_,@;
forall
k>q.
(2.4)
In the single series case, Box and Jenkins base model selection on the visual detection of patterns similar to (2.4) in the sample autocovariances. This can be quite tricky, and in the multivariate case it is practically impossible, as the patterns of interest in (2.4) are followed by m X m matrices of autocovariantes. Several authors, including Gray, Kelley, and McIntire [22], Beguin, Gourieroux, and Monfort [7], J en h ‘ns and Alavi [33], Woodward and Gray [54], and Glasbey [16], have used (2.4) in the development of statistics to aid model selection. However, as illustrated by Newbold and Bos [39] and Davies and Petruccelli [12], these statistics can have very heavy-tailed sampling distributions, making their practical interpretation extremely difficult. Our experience in analyzing real data sets suggests that a more useful approach to model selection can be based on a proposal of Hannan and Rissanen [26]. (For a recent extension, see [25].) We proceed in two stages. First, to obtain estimates of the innovations E, of (1.1) the ARMA process is
333
Testing Causality Using Vector ARMA Models
approximated by a pure autoregression of sufficiently high order. This can be accomplished by estimating (2.1) for all positive integers k, up to some maximum, which we generally fix at ten. This involves no additional computations, as the parameter estimates, and the estimated error covariance matrices fik, will already have been computed in the algorithm used to obtain the partial autocorrelations. The order of the approximating autoregression is chosen using the AIC criterion of Akaike [ 11. Thus, we choose the value of k for which
AIC( k) = log&]+
2km2
-
n
’
where n is the number of observations in each series, is smallest. The innovations of the underlying ARMA process are then estimated by 0, =x,
- fblx,_l- . . . - Qi_k >
(2.5)
where the 8ri are the estimated coefficient matrices of the approximating autoregression. It is well known that the AIC criterion is inconsistent, having in practice a tendency to pick over elaborate models. However, since the aim at this stage is to estimate the innovation series rather than find a “true” model, it may be desirable to employ here an order selection criterion that is, in this sense, conservative. At the second stage we regress, for various values of p and q, X t on XtPi (i = l,..., p) and &j (j=l,..., q), where the 2, are obtained from (2.5). Hence, we estimate models of the form
x, = @J_1
+ . . * + $X,_,
- @ii,_,
- . . . - O,O,_, +u~,~,~,
(2.6)
where up 9 t is an error term. The advantage of first deriving the estimated innovations’Q is, of course, that (2.6) is now easily estimated by ordinary least squares. Let G2,,, denote the estimated error convariance matrices from (2.6). Then, the BIC criterion of Akaike [2], Rissanen [45], and Schwarz [46] is used to suggest a possible model. Thus, for values of p and q (which, generally in our applied work have been taken not to exceed five) we compute logn BIC(p,q)
choosing that set (p,q)
=lod&,,I+(p
+ q)m2-
for which BIC is smallest.
n
’
(2.7)
334
PAUL NEWBOLD AND STEVEN M. HOTOPP
It is important to emphasize the use to which we put the statistics (2.7). The procedure is not treated as an orderestimation criterion, designed once and for all to give a unique model for subsequent efficient estimation. Rather, our approach is in the spirit of Box and Jenkins [8]. Borrowing these authors’ phraseology, we seek at this point a model to be “tentatively entertained,” but which may subsequently be abandoned if diagnostic checks suggest the advisability of such a course. The art of model selection, in this philosophy, is to develop inexpensively computed statistics whose sensible interpretation should lead the analyst to the choice of a model which constitutes a useful foundation on which to build. In this sense, we have found the statistics (2.7), used in conjunction with the sample autocorrelations and partial autocorrelations, to be extremely helpful. Of course, this approach also admits the possibility that other statistics might also profitably be taken into consideration. For example, recent proposals of Cooper and Wood [lo] and Tiao and Tsay [51] appear to merit further investigation. To illustrate, consider series of 135 quarterly seasonally adjusted observations on rates of change of employment in manufacturing (Xi) and hourly wages in manufacturing (X,) in the United States. The relationship between this pair of series has also been examined by Geary and Kennan [14]. The sample autocorrelation matrices contained large values, even for high lags, so that pure moving-average behavior is not indicated. The partial autocorrelation statistics (2.3), for lags up to 10, are displayed in Table 1. Comparison with tabulated values of the chi-square distribution with four degrees of freedom suggests that, if the generating process is pure autoregressive, the order should be at least two. However, the partial autocorrelation statistics for lags 3 and 4 are also moderately large.
TABLE 1 PARTIAL
AUTOCORRELATION
FOR EMPLOYMENT
k 1 2 3 4 5 6 7 8 9 10
AND
STATISTICS WAGES ak
121.99 13.94 6.21 8.61 2.24 2.62 2.97 2.64 7.43 0.53
DATA
Testing Causality Using Vector ARMA Models
335
Using the criterion (2.7), the ARMA(l, 1) model was selected. Since the partial autocorrelations suggest that, for a pure autoregressive structure, at least two parameter matrices would be required, the ARMA(l,l) model seems to be quite parsimonious, and it was decided to proceed with that model.
Parameter Estimation
Wilson [53] describes an estimation procedure that is essentially an extension of the least-squares algorithm employed by Box and Jenkins [8] in the single-series case. However, simulation evidence, reported by Ansley and Newbold [5], indicates a preference for full maximum likelihood in the estimation of the parameters of univariate ARMA models. It seems reasonable to conjecture that this preference should extend to the vector case. Direct derivations of the likelihood function are provided by Nicholls and Hall [40] and Hillmer and Tiao [27]. In our empirical work, we used a program based on the latter paper, written at the University of Wisconsin, and described in Tiao et al. [50]. The likelihood function is directly evaluated as a function of the unknown parameters, and maximized numerically. For the data on employment (X,) and wages (X,), we fitted the ARMA( 1,l) model
x, -
ipX,_l =a+Et-OE,_l,
where a is a vector of constants to allow for nonzero means. The parameter estimates (with estimated standard errors shown in parentheses) were 0.0087(0.0038) ’ =
0.0014(0.0008) &=
r
,.
1 [ >
’ =
0.43(0.12)
0.077(0.032)
- 0.28(0.12)
- 0.35(0.31)
0.069(0.051)
0.34(0.10)
- 0.53(0.25) O.SS(O.05)
1*
A computationally efficient alternative means of computing the likelihood function involves putting the model (1.1) into state-space form, and applying the Kalman filter. This approach is discussed by, for example, Ansley and Kohn [3], who also incorporate the possibility of missing observations. Model Checking
Following Box and Jenkins [8], we take the view that an initially chosen model should be subjected to checks on the adequacy with which it represents the data. Thus, for example, we do not see the criterion (2.7) as suggesting a
336
PAUL NEWBOLD AND STEVEN M. HOTOPP
final model, but simply as one of a number of statistics to be examined in selecting an initial model. The most commonly applied checks of time-series model adequacy are generally categorized in two distinct, though not unrelated, ways. First, as in [8], we might “over-fit’‘-that is, compare the originally chosen model with a more elaborate alternative. Again, since the innovations of a correctly specified model should be white noise, this can be checked through examination of the residual autocorrelations of the fitted model. Second, we might employ either tests of model adequacy designed to have high power against specific alternatives, or “portmanteau tests” which contemplate more general departures from model adequacy. However, as we will see, these distinctions are, in a sense, illusory. Suppose that the model (1.1) has been fitted to n sets of observations on m related times series. An obvious check on adequacy would be to test the null hypothesis of ARMA(p, 9) specification against the alternatives of a model involving K additional parameter matrices-say, an ARMA( p, 9 + K) model -for some positive integer K. One approach is to estimate the more elaborate model, and employ a likelihood-ratio test. However, estimation of heavily parametrized multivariate ARMA models can be computationally quite expensive. The estimation of the more general model can be avoided by using the Lagrange-multiplier, or “scores,” test of Rao [44] and Silvey [47]. This test is discussed in the present context for single-series models by Godfrey [17], Poskitt and Tremayne [41], and Hosking [28], and in the multivariate case by Hosking [30] and Poskitt and Tremayne [42]. Lagrange-multiplier tests are based on the derivatives of the log-likelihood function of the more general model, with respect to the parameters, evaluated at the maximum-likelihood estimates under the null hypothesis. For models of the form (Ll), the log likelihood can be written approximately, up to an additive constant, as
It then follows that l3L -
&;O-lvi,
= -
a9r.i.
j
j,t_,,
t
where +)r,i,j is the (i, j) element of @, and vi
.
j
3
t
-0v. 11,j.t
_l-
.
.
.
-‘q~,,j,t_q= -
Ei,jX,,
(2.8)
Testing Causality Using Vector ARMA Models where Ei j is the m X m matrix whose (i, j elements zero. Similarly, t3L ael,i.j
= -
337
) element is one, with all other
&;rlq,
j,t_l,
t
where
ui
j t
. -0~. 9
1.J.”
_ =E..E'J t’ 4
(2.9)
In practice, we obtain estimates G,,j, t and Q i, j, t by substituting in (2.8) and (2.9) the maximum-likelihood estimates under the null hypothesis Oj, and the residuals St from the fitted model. It is often possible to derive Lagrange-multiplier test statistics through fitting an auxiliary regression. Poskitt and Tremayne [42] show that, in the present instance, if the aim is to test the ARMA( p, 9) specification against an ARMA(p, 9 + K) alternative, this can be accomplished by considering the regression of 2’ = (E’,, E’,_ 1,. . . ) on
Under the null hypothesis that the ARMA(p, 9) specification is correct, the error covariance matrix for this regression is consistently estimated by I @ a, where fi is the estimated innovation covariance matrix derived from fitting (1.1) by maximum likelihood. The Lagrange-multiplier test statistic is then
Under the null hypothesis, this statistic has asymptotic chi-square distribution with Km’ degrees of freedom, that hypothesis being rejected for high values of (2.10). Poskitt and Tremayne further show that, asymptotically, an identical test results if the alternative specification is ARMA(p + K, Q), and that the test can be viewed as a pure significance test in the sense of Cox and Hinkley [ll, Chapter 31.
PAUL NEWBOLD AND STEVEN M. HOTOPP
338
To get some further insight into the form of the test statistic, we return to (2.8) and (2.9) and write ilL -
= -
&;wwyB)EijEt_I
ael,i, j
t
= - trO-‘(B)EijCEt-,E;~2’,
(2.11)
where
O(B)=&O,B-
.-.
-@,Bq
and B is the usual back-shift operator on the index of the time series. Let us write the residual autocovariance matrix from the fitted model (1.1) as
(2.12)
Then, it follows from (2.11) that the Lagrange multiplier test is essentially based on the first K of these matrices. We can therefore conclude, extending the corresponding result for the single-series case (shown for example in [38]), that if a check for model adequacy is to be based on the autocovariances of the residuals from the fitted ARMA( p, 4) model, the Lagrange-multiplier test provides an appropriate framework for such a check. Hosking [30] and Poskitt and Tremayne [42] also show how a portmanteau test of model adequacy, based on a moderately large number of residual autocovariance matrices (2.12), can be viewed as a Lagrange-multiplier test against a very general alternative specification. The test statistic is
(2.13) k-l
which, provided K is moderately large, has under the null hypothesis of correct model specification on asymptotic chi-square distribution with m2( K - p - q) degrees of freedom. This statistic was also derived directly from consideration of the asymptotic distribution of the residual autocovariantes by Hosking [29] and Li and McLeod [36]. Hosking [31] establishes the equivalence of alternative formulations of the test statistic.
Testing Causality
Using Vector ARMA Models
339
Following Hosking [29], we use the modified version of the portmanteau statistic
P*=
n2
2 (n - k)-‘tr(C?j$
(2.14)
k=l
In moderate-sized samples, simulation studies find that the empirical distribution of (2.14) is closer to chi-square than that of (2.13). As we have seen, for testing against alternative models, examination of individual residual autocorrelations is equivalent to employing Lagrange-multiplier tests. Nevertheless, we find it useful to look at these correlations, as their relative values can provide some insight into the possible sources of any model inadequacies. Thus, we compute, from the residual autocovariance matrices (2.12) the corresponding correlation matrices 8,. In addition to the checks just described, we have also found useful an informal procedure, where we compare the single-series models for each member of X, implied by the fitted vector ARMA model with those derived from analysis of the individual series alone. Although this comparison does not involve formal testing of hypotheses, we have found that it has the potential to pick up inadequacies of specification that might otherwise have been overlooked. To see how the single-series models can be derived from a multivariate ARMA( p, 9) model, write (1.1) as
O(B)X, = O(B)Et, where B is the back-shift operator. It then follows that
p(B) Ix, = @*(B)@(B)&,,
(2.15)
where IQ(B)] is the determinant and Q*(B) the adjoint matrix of Q(B). Now, the polynomial operator in B on the left-hand side of (2.15) is simply a scalar polynomial of degree mp in B. This, then, is the autoregressive operator of the single-series representation of each member of X,. Moreover, the matrix Q*(B)@(B) is of degree (m - 1)~ + 4 in B. As a consequence, it follows, using an argument given in Granger and Morris [20], that the model for each Xi t has a moving-average operator of order (m - 1)~ + 9. We can therefore write the ARMA[mp,(m - 1)~ + 91 univariate model for Xi t as
I@tB)
I’it
=
‘itBJait,
(2.16)
340
PAUL NEWBOLD AND STEVEN M. HOTOPP
where a it is white noise. The moving-average coefficients 13[(B) in (2.16) can be found as functions of the autoregressive and moving-average coefficients, and the innovation error variance, of the multivariate model (1.1). These results are employed as the basis of an approach to multivariate time-series model selection by Wallis [52] and Chan and Wallis [9]. In our own empirical work, the burden imposed on Equations (2.15) and (2.16) is rather less heavy. We appeal to them only to provide a supplementary informal check on the adequacy of representation of a fitted model. For the data on employment and wages, we fitted an ARMA(l,l) model. As a first step in checking the adequacy of this model, we show in Table 2 the residual autocorrelation matrices. As a crude check, we can compare the elements of these matrices with f 2/G = h 0.17. Viewed in this way, very nearly all of the 40 residual autocorrelations and cross-correlations reported in Table 2 are satisfactorily small. However, there is a very large correlation between 2,,, and gl,t_5. In principle, it is possible to add an additional term to our model to “get rid of” this awkward residual autocorrelation. However, in general, we do not favor such a course of action. Certainly, a better-fitting model could be achieved in this way, but, in our experience, such models, purpose-built to flatten residual autocorrelations at moderately high lags, often yield disappointingly poor forecasts. We feel that this attitude, which is rather more conservative than that adopted by many practitioners, is in the spirit of Box and Jenkins’s “parsimonious parametrization.” One possible explanation of this anomolous residual autocorrelation is that it is an artifact induced by seasonal adjustment of the data, as discussed for example by Grether and Nerlove [23]. There is also a large residual autocorrelation at lag 9, which again we are inclined to overlook in the interests of parsimony.
TABLE 2 RESIDUAL FROM
k 1 I
2
ARMA( 1,1)
AUTOCORRELATION
MODEL
MATRICES
TO EMPLOYMENT
k
[.oo 1 .02 .Ol [ .Ol
-.Ol .03 I .Ol - .04 I .02 .06
1 1 1
-
4
.06 [ - .02 .OO .06
5
--.28 .04
--.08 .07
3
F’ITI-ED
.03 .07
AND
WAGES
DATA
Testing Causality Using Vector ARMA Models
341
We computed the Lagrange-multiplier statistics (2.10) to test the null hypothesis of an ARMA(l, 1) specification against alternatives involving one or two additional moving-average parameter matrices. The statistics obtained were 2.52 and 7.23. Comparisons of these with tabulated values of the chi-square distribution with four and eight degrees of freedom provides little ground on which to question the adequacy of the fitted model. The modified portmanteau statistic (2.14) with K = 10,15,20, yielded 47.3, 62.1, and 77.2. The appropriate standard of comparison is the chi-square distribution with, respectively, 32, 52, and 72 degrees of freedom. The first of these test statistics is just significant at the 5-percent level. This reflects the large residual autocorrelations at lags five and nine. It is certainly feasible to produce a rather exotic modification of our model in response to this problem. However, as we have already indicated, our inclination is not to do so. As a further check, we use the parameters of the estimated multivariate ARMA( 1,l) model, together with the estimated innovation covariance matrix, to deduce, corresponding to (2.16), the univariate ARMA representations (ignoring the constant terms) (1 - 1.308B +0.417B2)X1,
= (1 - 0.599B - 0.212B2)a,,
(l-
= (l-
and 1.308B +0.417B2)X2,
0.728B +0.195B2)a2,.
Writing these in the form of infinite order autoregressions gives (1 - 0.709B
+0.204B2
- 0.028B3 +0.026B4
+O.OlOB’ + . . . )X1, = a,, (2.17)
and (1 - 0.580B
- 0.200B2
- 0.035B3 +0.015B4
+0.017B5
+ . . . )X2t = a21. (2.18)
Next, using the usual univariate time-series model-building methods, we independently constructed models for the two series, obtaining (l-
0.736B +0.216B2)X1,
= a,,
(2.19)
342
PAUL NEWBOLD AND STEVEN M. HOTOPP
and (1-
B)X,,
= (1 - 0.380B)a,,.
Transposing the second of these to the form of an infinite-order autoregression yields (1 - 0.620B
- 0.236B2 - 0.090B3 - 0.034B4
- 0.013B5 + . . . )Xst = ua1. (2.20)
The representation (2.17) is reassuringly similar to (2.19), and (2.18) is also very close to (2.20), suggesting no cause for alarm about the appropriateness of our multivariate model specification. We should add that the fifth sample autocorrelation of the series Xi was quite large, a factor we ignored, for reasons already discussed, in developing the univariate model. It might be said then that this check has also detected a possible source of model inadequacy, already revealed by examination of the residual autocorrelations. In this paper, we have illustrated a time-series model building methodology as applied to one pair of series. However, more extensive practical experience, reported in [32], leads us to believe that this general approach is likely to prove widely successful in the construction of vector ARMA models.
3.
CAUSALITY
TESTING
We now discuss the application of fitted vector ARMA models to the testing of causality, using the definition of causality proposed by Granger [ 181 and expanded upon in [21] and [19]. Granger’s work, together with that of Sims [48], stimulated considerable interest in causality testing in the economics literature, and many empirical studies have been published in the last few years. To introduce Granger’s definition of causality, let (Xi,, Xzt) be a purely nondeterministic process, on which a time series of observations is available. In principle, Xi, and Xzt could be vectors of time series, but for ease of exposition we treat them here as scalars. Consider, now, the problem of predicting Xi, n+ i, using information available at time n. Suppose that forecasts are to be based on the two information sets
ICI=
[Xl,n-j;
jaO]
Testing Causality
Using Vector ARMA Model.s
343
and
I,=
[Xi,n_j; jao,
i=1,2].
To obtain a testable definition of causality, attention is restricted to predictors that are linear functions of the members of these information sets, and forecast quality is judged in terms of expected squared error. Then, X, is said to cause Xi if and only if Xr++r is better predicted using information set I, than using I,. If both X, causes X, and X, causes X,, the pair of series are said to exhibit feedback. The most commonly applied tests for Granger causality are not based on parsimoniously parametrized models fitted to the available data. Rather, testing is generally carried out within the framework of a heavily parametrized structure. For example, what is sometimes called the Granger test, as it is implicit in [18], is based on fitting a high-order autoregressive model. Thus, the null hypothesis that X, does not cause Xi is checked through testing the hypothesis n, j = 0 ( j = 1,. . . , K) through the least-squares fitting of the regression K
K
‘1,=
C j=l
mljXl
‘-j+
’
C
1T2jx2
j=l
t-j+&,.
’
The autoregressive order K is often fixed arbitrarily at some moderately high value, with perhaps the application of a check that a larger value is not needed. The most commonly used tests of Granger causality are summarized and compared in [37] and [15]. Based on their simulation experiments, Nelson and Schwert [37] find marked gains in the power of the tests when they are based on the correct efficiently parametrized data-generating process rather than an overparametrized model. This encourages us to consider testing for causal ordering on the basis of multivariate ARMA models, fitted to the available data using the procedures described in the previous section of this paper. Let (Xi,, Xst ) be a pair of time series, generated by a multivariate ARMA(p, 4) model, which we now write as
[
h@)
@12(B)Xl,
@21(B)
@22(B)
11 X2t 1 =
40) [ fl2,u-Q
42(B) &2(B)
I[1 E1t
E2t
’
where the s@~~(B)and Bij(B) are scalar polynomials of respective degrees p
344
PAUL NEWBOLD AND STEVEN M. HOTOPP
and q in the back-shift operator B. Kang [34] has shown that a necessary and
sufficient condition for X, not to cause X,, in the Granger sense, is
Similarly, Xi does not Granger-cause X, if and only if
Notice, from (3.1), that the conditions
+,,(B) = 0,
MB) = 0
(3.2)
are, in general, sufficient but not necessary for noncausality. However, in the special case where the generating mechanism is either pure autoregressive or pure moving average, (3.1) and (3.2) are equivalent. At this point, our analysis diverges from that of Kang. We test directly for the restrictions implied by (3.1) on the parameters of the multivariate ARMA model. In effect this equation implies p + q restrictions on these parameters. One possible approach is to reestimate the model subject to these constraints, and apply a likelihood-ratio test. However, except in the special cases of pure autoregressive or pure moving-average models, our estimation program is not readily adapted to incorporate such constraints. An alternative approach, which is in any case computationally less expensive, is to employ the Wald test. This requires only information obtained from estimation of the unconstrained model. Let R denote the vector of autoregressive and moving-average parameters, and
f(P)=0
(3.3)
the p + q constraints on these parameters implied by (3.1). Further, let b denote the unconstrained maximum-likelihood estimators of R, with estimated covariance matrix Ea. With af/c?#l’ denoting the matrix of first partial derivatives from the left-hand side of (3.3), set
e,=(-$)&(-$)‘. Then, the Wald test statistic for the null hypothesis (3.3) is given by w=f(&‘e,lf@).
(3.4)
Testing Causality
Using Vector ARMA Models
345
Under the null hypothesis, the statistic (3.4) has an asymptotic chi-square distribution with p + 9 degrees of freedom, that hypothesis being rejected for high values of the statistic. We now apply this approach to the data on employment and wages. Writing the ARMA(l, 1) model so that
- h2B
I- h,B
l-&B
- %lB
Xl,
1 - 8,,B
- e,2B
I[ Xzt 1 =
- 8,,B
1 - e,,B
[
the necessary and sufficient conditions (3.1) for X, Granger sense, X, are -
I[1’ E1t Ezt
not to cause, in the
&B(l - %.-$)+ 4,B(l- +zzB)= 0
4, - h2 = 0,
h2e22- e12+,,= 0.
This pair of restrictions can be tested using the Wald statistic (3.4), which employs the estimated covariance matrix of the parameter estimators. This matrix is computed by our estimation program. For these data the computed value of the test statistic was 4.718. Comparison with tabulated values of the chi-square distribution for two degrees of freedom reveals that the null hypothesis that wages do not cause employment can be rejected at the lo-percent but not at the Spercent significance level. We further found that the null hypothesis that employment does not cause wages can be rejected at the 5-percent level. The computed test statistic (3.4) was 7.204. This research was partially funded under contract number SES-8011175.
by the National
Science Foundation,
REFERENCES 1 2
3 4
H. Akaike, Fitting autoregressive models for prediction, Ann. Inst. Statist. Math. 21:243-247 (1969). H. Akaike, On entropy maximization principle, in Applications of Statistics (P. R. Krishnaiah, Ed.), North Holland, Amsterdam, 1977, pp. 27-41. C. F. Ansiey and R. Kohn, Exact likelihood of vector autoregressive-moving average process with missing or aggregated data, Biometiku 70:275-278 (1983). C. F. Ansley and P. Newbold, Multivariate partial autocorrelations, Proc. Bus. Econom. Statist., Amer. Statist. Assoc., 1979, pp. 349-353.
346
8 9 10 11 12
13 14 15
16 17 18 19 20 21 22 23 24 25 26
PAUL NEWBOLD
AND STEVEN M. HOTOPP
C. F. Ansley and P. Newbold, Finite sample properties of estimators for autoregressive moving average models, J. Econometrics 13:159-183 (1980). M. S. Bartlett, On the theoretical specification of sampling properties of autocorrelated time series, I. Roy. Statist. Sot. Ser. I? 8:27-41 (1946). J. M. Beguin, C. Gourieroux, and A. Monfort, Identification of a mixed autoregressive-moving average process: The corner method, Time Series (0. D. Anderson, Ed.), North Holland, Amsterdam, 1980, pp. 423436. G. E. P. Box and G. M. Jenkins, Time Series Analysis, Forecasting and Control, Holden Day, San Francisco, 1970. W. Y. T. Chan and K. F. WaIlis, Multiple time series modelling: Another look at the mink-muskrat interaction, Appl. Statist. 27: 168- 175 (1978). D. M. Cooper and E. F. Wood, Identifying multivariate time series models, J. Time Ser. Anal. 3:153-164 (1982). D. R. Cox and D. V. Hinkley, Theoretical Statistics, Chapman and Hall, London, 1974. N. Davies and J. D. Petruccelli, On the use of the general partial autocorrelation function for order determination in ARMA(p, 9) processes, J. Amer. Statist. Assoc. 79:374-377 (1984). J. Durbin, The fitting of time series models, Reu. Inst. Internat. Statist. 28:233-244 (1960). P. T. Geary and J. Kennan, The employment-real wage relationship: An intemational study, J. Polit. Econ. 90:85&871 (1982). J. Geweke, R. Meese, and W. Dent, Comparing alternative tests of causality in temporal systems: Analytic results and experimental evidence, J. Econometrics 21:161-194 (1983). C. A. Glasbey, A generalization of partial autocorrelations useful in identifying ARMA models, Technometrics 241223-228 (1982). L. G. Godfrey, Testing the adequacy of a time series model, Biometrika 66:67-72 (1979). C. W. J. Granger, Investigating causal relations by econometric models and cross-spectral methods, Econometrica 37:424-438 (1969). C. W. J. Granger, Testing for causality: A personal viewpoint, J. Econom. Dynamics Control 2:329-352 (1980). C. W. J. Granger and M. J. Morris, Time series modelling and interpretation, J. Roy. Statist. Sot. Ser. A 139:246-257 (1976). C. W. J. Granger and P. Newbold, Forecasting Economic Time Series, Academic, New York, 1977. H. L. Gray, G. D. KelIey, D. D. McIntire, A new approach to ARMA modelling, Comm. Statist. B-Simulation Comput. 7:1-77 (1978). D. M. Grether and M. Nerlove, Some properties of “optimal” seasonal adjustment, Econumetrica 38:682-703 (1970). E. J. Hannan, Multiple Time Series, Wiley, New York, 1970. E. J. Hannan and L. Kavalieris, A method for autoregressive-moving average estimation, Biometrika 71:273-280 (1984). E. J. Hannan and J. Rissanen, Recursive estimation of mixed autoregressive-moving average order, Biometrika 69:81-94 (1982); Correction, 70:303.
Testing Causality Using Vector ARMA Models 27 28 29 30 31 32 33 34 35 36
37
38 39 40 41 42 43 44
45 46 47 48
347
S. C. Hilhner and G. C. Tiao, Likelihood function of stationary multiple autoregressive moving average models, J. Amer. Statist. Assoc. 74:652-660 (1979). J. R. M. Hosking, Lagrange multiplier tests of time series models, J. Roy. Statist. Sot. Ser. B 42:170-181 (1980). J. Il. M. Hosking, The multivariate portmanteau statistic, J. Amer. Statist. Assoc. 75:602-608 (1980). J. R. M. Hosking, Lagrange multiplier tests of multivariate time series models, J. Roy. Statist. Sot. Ser. B 43:219-230 (1981). J. R. M. Hosking, Equivalent forms of the multivariate portmanteau statistic, J. Roy. Statist. Sot. Ser. B 43:261-262 (1981). S. M. Hotopp, Practical methods in multivariate time series analysis, Ph.D. Thesis, Dept. of Economics, Univ. of Illinois, Champaign-Urbana, 1985. G. M. Jenkins and A. S. Alavi, Some aspects of modelling and forecasting multivariate time series, J. Time Ser. Anal. 2:1-47 (1981). H. Kang, Necessary and sufficient conditions for causality testing in multivariate ARMA models, J. Time Ser. AnaE. 2:95-101 (1981). J. Ledolter, The analysis of multivariate time series applied to problems in hydrology, J. Hydrol. 36:327-352 (1978). W. K. Li and A. I. McLeod, Distribution of the residual autocorrelations in multivariate ARMA time series models. J. Roy. Statist. Sot. Ser. B 43:231-239 (1981). C. R. Nelson and G. W. Schwert, Tests for predictive relationships between time series variables: A Monte Carlo investigation, J. Amer. Statist. Assoc. 77:11-18 (1982). P. Newbold, The equivalence of two tests for time series model adequacy, Biometrika 67:463-465 (1980). P. Newbold and T. Bos, On q-conditioned partial correlations, J. Time Ser. Anal. 45%55 (1983). D. F. Nicholls and A. D. HalI, The exact likelihood function of multivariate autoregressive-moving average models, Biumetriku, 66:259-264 (1979). D. S. Poskitt and A. R. Tremayne, Testing the specification of a fitted ARMA model, Biumetrika 67:359-363 (1980). D. S. Poskitt and A. R. Tremayne, Diagnostic tests for multiple time series models, Ann. Statist. 10:114-120 (1982). M. H. QuenouiIle, The Analysis of Multiple Time Series, Griffin, London, 1957. C. R. Rao, Large sample tests of statistical hypotheses concerning several parameters with applications to problems of estimation, Proc. Cambridge Philos. Sot. 44:50-59 (1948). J. Rissanen, Modeling by shortest data description, Autumatica 14:465-471 (1978). G. Schwarz, Estimating the dimension of a model, Ann. Statist. 6:461-464 (1978). S. D. Silvey, The Lagrangian multiplier test, Ann. Math. Statist. 30:389-407 (1959). C. A. Sims, Money, income and causality, Am. Econom. Rev. 62:540-552 (1972).
348
PAUL NEWBOLD
49
G. C. Tiao and G. E. P. Box, Modeling multiple time series with applications, J. Amer. Statist. Assoc. 76802-816 (1981). G. C. Tiao, G. E. P. Box, M. R. Grupe, G. B. Hudak, W. R. Bell, and I. Chang, The Wisconsin multiple time series (WMTS-1) program: A preliminary guide, Dept. of Statistics, Univ. of Wisconsin, Madison, 1979. G. C. Tiao and R. S. Tsay, Multiple time series modeling and extended sample cross-correlations, J. Bus. and Ecunon. Statist. 1:43-56 (1983). K. F. WaIIis, Multiple time series analysis and the final form of econometric models, Eccmumetrica 45: 1481- 1497 (1977). G. T. Wilson, The estimation of parameters in multivariate time series models, J. Roy. Statist. Sot. Ser. B 35:76-85 (1973). W. A. Woodward and H. L. Gray, On the relationship between the S array and the Box-Jenkins method of ARMA model identification, J. Amer. Statist. Assoc. 76:579-587 (1981).
59
51 52 53 54
AND STEVEN M. HOTOPP
ABSTRACT We outline a practical methodology for fitting to data vector autoregressive moving-average models, and employ these models to test for Granger causality. The procedures are applied to data on employment and hourly wages in manufacturing in
the United States.
1.
INTRODUCTION Given observations through time on m series X{ = (Xi,, X,,, . . . , Xm(), we
wish
to
fit
a
member
of
the
vector
autoregressive
moving-average
[ARMA( p, q)] class of models
where the api and Oj are m X m matrices of fixed coefficients, vector white noise, so that
and et is
The model (1.1) is an obvious extension to the multivariate case of the univariate models discussed by Box and Jenkins [8]. For the single-series case, Box and Jenkins outlined a model-building methodology involving an iterative cycle of model selection, parameter estimation, and model checking. At the first stage, based on judgmental interpretation of statistics readily computed APPLIED
MATHEMATZCS AND COMPVTATZON 20: 329-348
(1986)
329
0 Elsevier Science Publishing Co., Inc., 1986 52 Vanderbilt Ave., New York, NY 10017
009fUOO3/86/$03.50
PAUL NEWBOLD AND STEVEN M. HOTOPP
330
from the data, a specific model from the general class is chosen for further analysis. Next, using statistically efficient methods, the fixed parameters of the selected model are estimated. Finally, diagnostic checks on the adequacy of representation of the model to the available data are applied. Any inadequacies revealed at this stage may suggest an alternative structure, and the cycle of selection, estimation, and checking is iterated until a satisfactory model is achieved. Some early work on the fitting to data of models of the form (1.1) is reported by Quenouille [43]. More recently, Jenkins and Alavi [33] and Tiao and Box [49] have discussed the extension of the Box-Jenkins approach to this problem. In similar spirit, we report in Section 2 of this paper on a model-building methodology that we have successfully applied to several small sets of economic time series. These procedures will be illustrated by reference to data on employment and hourly wages in manufacturing in the United States. Experience with other data sets is discussed in [32]. In Section 3 we discuss the use of fitted vector autoregressivemoving average models to testing for “causality,” in the sense of Granger [18]. Such tests are likely to have higher power than generally applied procedures based on heavily parametrized structures. 2.
MODEL
BUILDING
In this section we consider the three stages of the model-building cycle, when the objective is to fit to data a model from the class (1.1). Model Selection In practice, even for the case of a single time series, it is the initial selection of an appropriate model from the general class that causes most difficulty in the application of the Box-Jenkins methodology. This difficulty is considerably more acute in the multivariate case. For the analysis of a single series, Box and Jenkins base model selection on examination of the sample autocorrelations and partial autocorrelations of the series. Analogous statistics are also helpful in the multivariate case. Define the population autocovariance matrices of X, as
r,c=E(X,X;+,) Then it is well known that, for a pure moving-average process of order q [p = 0 in (l.l)], rk = 0
forall
k>q.
Testing Causality
Using Vector ARMA Models
331
This suggests that the corresponding sample autocovariance matrices
should be useful in identifying pure moving-average behavior. It is easier to interpret the corresponding autocorrelation matrices R,, whose (i, j) elements are given by
where ck; i, jis the (i, j) element of Ck. For checking for pure moving average behavior of order 4, since the sample autocorrelations have asymptotic normal distributions, comparison can be made with the asymptotic standard errors, given by Bartlett [S] as l/2
4
SE[Tk;i,j]
=
12-
1+2
1 11 =
Ph;i,iPh;j,j
1
I
1
where, in practice, the population autocorrelations in this expression are replaced by the corresponding sample quantities. As in the case of a single time series, pure autoregressive models [q = 0 in (l.l)] can be detected through the sample partial autocorrelations. Consider theregressionofX,onX,_j(j=1,2,...,k): x,
=
@klXf-l+
@k2Xt-2
+
”
’ +
@kkXt-k
+Uk,t
(2.1)
where the error term uk t has covariance matrix Qk, and also the regression of X, on Xt+i (j=l,2,...:k)
where uz t has covariance matrix Qt. I_& A, and A$ denote the symmetric square roots of the matrices at, and Dz. Then, Quenouille [43, p. 411 defines the multiple partial autocorrelation of order k + 1 as P k+l=
&‘E(uk
7tUf’t-k-l)A:-l. 1
332
PAUL NEWBOLD AND STEVEN M. HOTOPP
If the generating process is autoregressive of order k, then u k,t and uf, t k_ 1 will be uncorrelated, and hence P k+l equal to zero. Ledolter [35] and Ansley and Newbold [4] have extended the algorithm of Durbin [13] to the computation of estimates i),, i of the partial autocorrelation matrices. Hannan [24, p. 3981 showsnthat, for an autoregressive process of order k, the elements of the matrices Pk+ j (j > 1) have independent asymptotic normal distributions, with zero means and variances n- ‘. Thus, for such a process, the statistics =
(Ik+j
ntrbk+ji);+j
have asymptotic chi-square distributions with m2 degrees of freedom. Altematively, these statistics can be expressed as (2.3)
a k+jzfltr’k+j,k+j’k*+j,k+j9
where the &‘s are the estimates of the parameter matrices in (2.1) and (2.2). The quantities (2.3) can therefore be employed in the detection of pure autoregressive models. Much more difficult is the distinction, based on statistics conveniently computed from the data, among mixed ARMA(p, o) models. To see why this is so, consider the natural extension of the Box-Jenkins single series approach to model selection to the vector case. Postmultiplying through (1.1) by Xi_, and taking expectations yields
r, = l-k_@;
+ .
’ ’
+
r,_,@;
forall
k>q.
(2.4)
In the single series case, Box and Jenkins base model selection on the visual detection of patterns similar to (2.4) in the sample autocovariances. This can be quite tricky, and in the multivariate case it is practically impossible, as the patterns of interest in (2.4) are followed by m X m matrices of autocovariantes. Several authors, including Gray, Kelley, and McIntire [22], Beguin, Gourieroux, and Monfort [7], J en h ‘ns and Alavi [33], Woodward and Gray [54], and Glasbey [16], have used (2.4) in the development of statistics to aid model selection. However, as illustrated by Newbold and Bos [39] and Davies and Petruccelli [12], these statistics can have very heavy-tailed sampling distributions, making their practical interpretation extremely difficult. Our experience in analyzing real data sets suggests that a more useful approach to model selection can be based on a proposal of Hannan and Rissanen [26]. (For a recent extension, see [25].) We proceed in two stages. First, to obtain estimates of the innovations E, of (1.1) the ARMA process is
333
Testing Causality Using Vector ARMA Models
approximated by a pure autoregression of sufficiently high order. This can be accomplished by estimating (2.1) for all positive integers k, up to some maximum, which we generally fix at ten. This involves no additional computations, as the parameter estimates, and the estimated error covariance matrices fik, will already have been computed in the algorithm used to obtain the partial autocorrelations. The order of the approximating autoregression is chosen using the AIC criterion of Akaike [ 11. Thus, we choose the value of k for which
AIC( k) = log&]+
2km2
-
n
’
where n is the number of observations in each series, is smallest. The innovations of the underlying ARMA process are then estimated by 0, =x,
- fblx,_l- . . . - Qi_k >
(2.5)
where the 8ri are the estimated coefficient matrices of the approximating autoregression. It is well known that the AIC criterion is inconsistent, having in practice a tendency to pick over elaborate models. However, since the aim at this stage is to estimate the innovation series rather than find a “true” model, it may be desirable to employ here an order selection criterion that is, in this sense, conservative. At the second stage we regress, for various values of p and q, X t on XtPi (i = l,..., p) and &j (j=l,..., q), where the 2, are obtained from (2.5). Hence, we estimate models of the form
x, = @J_1
+ . . * + $X,_,
- @ii,_,
- . . . - O,O,_, +u~,~,~,
(2.6)
where up 9 t is an error term. The advantage of first deriving the estimated innovations’Q is, of course, that (2.6) is now easily estimated by ordinary least squares. Let G2,,, denote the estimated error convariance matrices from (2.6). Then, the BIC criterion of Akaike [2], Rissanen [45], and Schwarz [46] is used to suggest a possible model. Thus, for values of p and q (which, generally in our applied work have been taken not to exceed five) we compute logn BIC(p,q)
choosing that set (p,q)
=lod&,,I+(p
+ q)m2-
for which BIC is smallest.
n
’
(2.7)
334
PAUL NEWBOLD AND STEVEN M. HOTOPP
It is important to emphasize the use to which we put the statistics (2.7). The procedure is not treated as an orderestimation criterion, designed once and for all to give a unique model for subsequent efficient estimation. Rather, our approach is in the spirit of Box and Jenkins [8]. Borrowing these authors’ phraseology, we seek at this point a model to be “tentatively entertained,” but which may subsequently be abandoned if diagnostic checks suggest the advisability of such a course. The art of model selection, in this philosophy, is to develop inexpensively computed statistics whose sensible interpretation should lead the analyst to the choice of a model which constitutes a useful foundation on which to build. In this sense, we have found the statistics (2.7), used in conjunction with the sample autocorrelations and partial autocorrelations, to be extremely helpful. Of course, this approach also admits the possibility that other statistics might also profitably be taken into consideration. For example, recent proposals of Cooper and Wood [lo] and Tiao and Tsay [51] appear to merit further investigation. To illustrate, consider series of 135 quarterly seasonally adjusted observations on rates of change of employment in manufacturing (Xi) and hourly wages in manufacturing (X,) in the United States. The relationship between this pair of series has also been examined by Geary and Kennan [14]. The sample autocorrelation matrices contained large values, even for high lags, so that pure moving-average behavior is not indicated. The partial autocorrelation statistics (2.3), for lags up to 10, are displayed in Table 1. Comparison with tabulated values of the chi-square distribution with four degrees of freedom suggests that, if the generating process is pure autoregressive, the order should be at least two. However, the partial autocorrelation statistics for lags 3 and 4 are also moderately large.
TABLE 1 PARTIAL
AUTOCORRELATION
FOR EMPLOYMENT
k 1 2 3 4 5 6 7 8 9 10
AND
STATISTICS WAGES ak
121.99 13.94 6.21 8.61 2.24 2.62 2.97 2.64 7.43 0.53
DATA
Testing Causality Using Vector ARMA Models
335
Using the criterion (2.7), the ARMA(l, 1) model was selected. Since the partial autocorrelations suggest that, for a pure autoregressive structure, at least two parameter matrices would be required, the ARMA(l,l) model seems to be quite parsimonious, and it was decided to proceed with that model.
Parameter Estimation
Wilson [53] describes an estimation procedure that is essentially an extension of the least-squares algorithm employed by Box and Jenkins [8] in the single-series case. However, simulation evidence, reported by Ansley and Newbold [5], indicates a preference for full maximum likelihood in the estimation of the parameters of univariate ARMA models. It seems reasonable to conjecture that this preference should extend to the vector case. Direct derivations of the likelihood function are provided by Nicholls and Hall [40] and Hillmer and Tiao [27]. In our empirical work, we used a program based on the latter paper, written at the University of Wisconsin, and described in Tiao et al. [50]. The likelihood function is directly evaluated as a function of the unknown parameters, and maximized numerically. For the data on employment (X,) and wages (X,), we fitted the ARMA( 1,l) model
x, -
ipX,_l =a+Et-OE,_l,
where a is a vector of constants to allow for nonzero means. The parameter estimates (with estimated standard errors shown in parentheses) were 0.0087(0.0038) ’ =
0.0014(0.0008) &=
r
,.
1 [ >
’ =
0.43(0.12)
0.077(0.032)
- 0.28(0.12)
- 0.35(0.31)
0.069(0.051)
0.34(0.10)
- 0.53(0.25) O.SS(O.05)
1*
A computationally efficient alternative means of computing the likelihood function involves putting the model (1.1) into state-space form, and applying the Kalman filter. This approach is discussed by, for example, Ansley and Kohn [3], who also incorporate the possibility of missing observations. Model Checking
Following Box and Jenkins [8], we take the view that an initially chosen model should be subjected to checks on the adequacy with which it represents the data. Thus, for example, we do not see the criterion (2.7) as suggesting a
336
PAUL NEWBOLD AND STEVEN M. HOTOPP
final model, but simply as one of a number of statistics to be examined in selecting an initial model. The most commonly applied checks of time-series model adequacy are generally categorized in two distinct, though not unrelated, ways. First, as in [8], we might “over-fit’‘-that is, compare the originally chosen model with a more elaborate alternative. Again, since the innovations of a correctly specified model should be white noise, this can be checked through examination of the residual autocorrelations of the fitted model. Second, we might employ either tests of model adequacy designed to have high power against specific alternatives, or “portmanteau tests” which contemplate more general departures from model adequacy. However, as we will see, these distinctions are, in a sense, illusory. Suppose that the model (1.1) has been fitted to n sets of observations on m related times series. An obvious check on adequacy would be to test the null hypothesis of ARMA(p, 9) specification against the alternatives of a model involving K additional parameter matrices-say, an ARMA( p, 9 + K) model -for some positive integer K. One approach is to estimate the more elaborate model, and employ a likelihood-ratio test. However, estimation of heavily parametrized multivariate ARMA models can be computationally quite expensive. The estimation of the more general model can be avoided by using the Lagrange-multiplier, or “scores,” test of Rao [44] and Silvey [47]. This test is discussed in the present context for single-series models by Godfrey [17], Poskitt and Tremayne [41], and Hosking [28], and in the multivariate case by Hosking [30] and Poskitt and Tremayne [42]. Lagrange-multiplier tests are based on the derivatives of the log-likelihood function of the more general model, with respect to the parameters, evaluated at the maximum-likelihood estimates under the null hypothesis. For models of the form (Ll), the log likelihood can be written approximately, up to an additive constant, as
It then follows that l3L -
&;O-lvi,
= -
a9r.i.
j
j,t_,,
t
where +)r,i,j is the (i, j) element of @, and vi
.
j
3
t
-0v. 11,j.t
_l-
.
.
.
-‘q~,,j,t_q= -
Ei,jX,,
(2.8)
Testing Causality Using Vector ARMA Models where Ei j is the m X m matrix whose (i, j elements zero. Similarly, t3L ael,i.j
= -
337
) element is one, with all other
&;rlq,
j,t_l,
t
where
ui
j t
. -0~. 9
1.J.”
_ =E..E'J t’ 4
(2.9)
In practice, we obtain estimates G,,j, t and Q i, j, t by substituting in (2.8) and (2.9) the maximum-likelihood estimates under the null hypothesis Oj, and the residuals St from the fitted model. It is often possible to derive Lagrange-multiplier test statistics through fitting an auxiliary regression. Poskitt and Tremayne [42] show that, in the present instance, if the aim is to test the ARMA( p, 9) specification against an ARMA(p, 9 + K) alternative, this can be accomplished by considering the regression of 2’ = (E’,, E’,_ 1,. . . ) on
Under the null hypothesis that the ARMA(p, 9) specification is correct, the error covariance matrix for this regression is consistently estimated by I @ a, where fi is the estimated innovation covariance matrix derived from fitting (1.1) by maximum likelihood. The Lagrange-multiplier test statistic is then
Under the null hypothesis, this statistic has asymptotic chi-square distribution with Km’ degrees of freedom, that hypothesis being rejected for high values of (2.10). Poskitt and Tremayne further show that, asymptotically, an identical test results if the alternative specification is ARMA(p + K, Q), and that the test can be viewed as a pure significance test in the sense of Cox and Hinkley [ll, Chapter 31.
PAUL NEWBOLD AND STEVEN M. HOTOPP
338
To get some further insight into the form of the test statistic, we return to (2.8) and (2.9) and write ilL -
= -
&;wwyB)EijEt_I
ael,i, j
t
= - trO-‘(B)EijCEt-,E;~2’,
(2.11)
where
O(B)=&O,B-
.-.
-@,Bq
and B is the usual back-shift operator on the index of the time series. Let us write the residual autocovariance matrix from the fitted model (1.1) as
(2.12)
Then, it follows from (2.11) that the Lagrange multiplier test is essentially based on the first K of these matrices. We can therefore conclude, extending the corresponding result for the single-series case (shown for example in [38]), that if a check for model adequacy is to be based on the autocovariances of the residuals from the fitted ARMA( p, 4) model, the Lagrange-multiplier test provides an appropriate framework for such a check. Hosking [30] and Poskitt and Tremayne [42] also show how a portmanteau test of model adequacy, based on a moderately large number of residual autocovariance matrices (2.12), can be viewed as a Lagrange-multiplier test against a very general alternative specification. The test statistic is
(2.13) k-l
which, provided K is moderately large, has under the null hypothesis of correct model specification on asymptotic chi-square distribution with m2( K - p - q) degrees of freedom. This statistic was also derived directly from consideration of the asymptotic distribution of the residual autocovariantes by Hosking [29] and Li and McLeod [36]. Hosking [31] establishes the equivalence of alternative formulations of the test statistic.
Testing Causality
Using Vector ARMA Models
339
Following Hosking [29], we use the modified version of the portmanteau statistic
P*=
n2
2 (n - k)-‘tr(C?j$
(2.14)
k=l
In moderate-sized samples, simulation studies find that the empirical distribution of (2.14) is closer to chi-square than that of (2.13). As we have seen, for testing against alternative models, examination of individual residual autocorrelations is equivalent to employing Lagrange-multiplier tests. Nevertheless, we find it useful to look at these correlations, as their relative values can provide some insight into the possible sources of any model inadequacies. Thus, we compute, from the residual autocovariance matrices (2.12) the corresponding correlation matrices 8,. In addition to the checks just described, we have also found useful an informal procedure, where we compare the single-series models for each member of X, implied by the fitted vector ARMA model with those derived from analysis of the individual series alone. Although this comparison does not involve formal testing of hypotheses, we have found that it has the potential to pick up inadequacies of specification that might otherwise have been overlooked. To see how the single-series models can be derived from a multivariate ARMA( p, 9) model, write (1.1) as
O(B)X, = O(B)Et, where B is the back-shift operator. It then follows that
p(B) Ix, = @*(B)@(B)&,,
(2.15)
where IQ(B)] is the determinant and Q*(B) the adjoint matrix of Q(B). Now, the polynomial operator in B on the left-hand side of (2.15) is simply a scalar polynomial of degree mp in B. This, then, is the autoregressive operator of the single-series representation of each member of X,. Moreover, the matrix Q*(B)@(B) is of degree (m - 1)~ + 4 in B. As a consequence, it follows, using an argument given in Granger and Morris [20], that the model for each Xi t has a moving-average operator of order (m - 1)~ + 9. We can therefore write the ARMA[mp,(m - 1)~ + 91 univariate model for Xi t as
I@tB)
I’it
=
‘itBJait,
(2.16)
340
PAUL NEWBOLD AND STEVEN M. HOTOPP
where a it is white noise. The moving-average coefficients 13[(B) in (2.16) can be found as functions of the autoregressive and moving-average coefficients, and the innovation error variance, of the multivariate model (1.1). These results are employed as the basis of an approach to multivariate time-series model selection by Wallis [52] and Chan and Wallis [9]. In our own empirical work, the burden imposed on Equations (2.15) and (2.16) is rather less heavy. We appeal to them only to provide a supplementary informal check on the adequacy of representation of a fitted model. For the data on employment and wages, we fitted an ARMA(l,l) model. As a first step in checking the adequacy of this model, we show in Table 2 the residual autocorrelation matrices. As a crude check, we can compare the elements of these matrices with f 2/G = h 0.17. Viewed in this way, very nearly all of the 40 residual autocorrelations and cross-correlations reported in Table 2 are satisfactorily small. However, there is a very large correlation between 2,,, and gl,t_5. In principle, it is possible to add an additional term to our model to “get rid of” this awkward residual autocorrelation. However, in general, we do not favor such a course of action. Certainly, a better-fitting model could be achieved in this way, but, in our experience, such models, purpose-built to flatten residual autocorrelations at moderately high lags, often yield disappointingly poor forecasts. We feel that this attitude, which is rather more conservative than that adopted by many practitioners, is in the spirit of Box and Jenkins’s “parsimonious parametrization.” One possible explanation of this anomolous residual autocorrelation is that it is an artifact induced by seasonal adjustment of the data, as discussed for example by Grether and Nerlove [23]. There is also a large residual autocorrelation at lag 9, which again we are inclined to overlook in the interests of parsimony.
TABLE 2 RESIDUAL FROM
k 1 I
2
ARMA( 1,1)
AUTOCORRELATION
MODEL
MATRICES
TO EMPLOYMENT
k
[.oo 1 .02 .Ol [ .Ol
-.Ol .03 I .Ol - .04 I .02 .06
1 1 1
-
4
.06 [ - .02 .OO .06
5
--.28 .04
--.08 .07
3
F’ITI-ED
.03 .07
AND
WAGES
DATA
Testing Causality Using Vector ARMA Models
341
We computed the Lagrange-multiplier statistics (2.10) to test the null hypothesis of an ARMA(l, 1) specification against alternatives involving one or two additional moving-average parameter matrices. The statistics obtained were 2.52 and 7.23. Comparisons of these with tabulated values of the chi-square distribution with four and eight degrees of freedom provides little ground on which to question the adequacy of the fitted model. The modified portmanteau statistic (2.14) with K = 10,15,20, yielded 47.3, 62.1, and 77.2. The appropriate standard of comparison is the chi-square distribution with, respectively, 32, 52, and 72 degrees of freedom. The first of these test statistics is just significant at the 5-percent level. This reflects the large residual autocorrelations at lags five and nine. It is certainly feasible to produce a rather exotic modification of our model in response to this problem. However, as we have already indicated, our inclination is not to do so. As a further check, we use the parameters of the estimated multivariate ARMA( 1,l) model, together with the estimated innovation covariance matrix, to deduce, corresponding to (2.16), the univariate ARMA representations (ignoring the constant terms) (1 - 1.308B +0.417B2)X1,
= (1 - 0.599B - 0.212B2)a,,
(l-
= (l-
and 1.308B +0.417B2)X2,
0.728B +0.195B2)a2,.
Writing these in the form of infinite order autoregressions gives (1 - 0.709B
+0.204B2
- 0.028B3 +0.026B4
+O.OlOB’ + . . . )X1, = a,, (2.17)
and (1 - 0.580B
- 0.200B2
- 0.035B3 +0.015B4
+0.017B5
+ . . . )X2t = a21. (2.18)
Next, using the usual univariate time-series model-building methods, we independently constructed models for the two series, obtaining (l-
0.736B +0.216B2)X1,
= a,,
(2.19)
342
PAUL NEWBOLD AND STEVEN M. HOTOPP
and (1-
B)X,,
= (1 - 0.380B)a,,.
Transposing the second of these to the form of an infinite-order autoregression yields (1 - 0.620B
- 0.236B2 - 0.090B3 - 0.034B4
- 0.013B5 + . . . )Xst = ua1. (2.20)
The representation (2.17) is reassuringly similar to (2.19), and (2.18) is also very close to (2.20), suggesting no cause for alarm about the appropriateness of our multivariate model specification. We should add that the fifth sample autocorrelation of the series Xi was quite large, a factor we ignored, for reasons already discussed, in developing the univariate model. It might be said then that this check has also detected a possible source of model inadequacy, already revealed by examination of the residual autocorrelations. In this paper, we have illustrated a time-series model building methodology as applied to one pair of series. However, more extensive practical experience, reported in [32], leads us to believe that this general approach is likely to prove widely successful in the construction of vector ARMA models.
3.
CAUSALITY
TESTING
We now discuss the application of fitted vector ARMA models to the testing of causality, using the definition of causality proposed by Granger [ 181 and expanded upon in [21] and [19]. Granger’s work, together with that of Sims [48], stimulated considerable interest in causality testing in the economics literature, and many empirical studies have been published in the last few years. To introduce Granger’s definition of causality, let (Xi,, Xzt) be a purely nondeterministic process, on which a time series of observations is available. In principle, Xi, and Xzt could be vectors of time series, but for ease of exposition we treat them here as scalars. Consider, now, the problem of predicting Xi, n+ i, using information available at time n. Suppose that forecasts are to be based on the two information sets
ICI=
[Xl,n-j;
jaO]
Testing Causality
Using Vector ARMA Model.s
343
and
I,=
[Xi,n_j; jao,
i=1,2].
To obtain a testable definition of causality, attention is restricted to predictors that are linear functions of the members of these information sets, and forecast quality is judged in terms of expected squared error. Then, X, is said to cause Xi if and only if Xr++r is better predicted using information set I, than using I,. If both X, causes X, and X, causes X,, the pair of series are said to exhibit feedback. The most commonly applied tests for Granger causality are not based on parsimoniously parametrized models fitted to the available data. Rather, testing is generally carried out within the framework of a heavily parametrized structure. For example, what is sometimes called the Granger test, as it is implicit in [18], is based on fitting a high-order autoregressive model. Thus, the null hypothesis that X, does not cause Xi is checked through testing the hypothesis n, j = 0 ( j = 1,. . . , K) through the least-squares fitting of the regression K
K
‘1,=
C j=l
mljXl
‘-j+
’
C
1T2jx2
j=l
t-j+&,.
’
The autoregressive order K is often fixed arbitrarily at some moderately high value, with perhaps the application of a check that a larger value is not needed. The most commonly used tests of Granger causality are summarized and compared in [37] and [15]. Based on their simulation experiments, Nelson and Schwert [37] find marked gains in the power of the tests when they are based on the correct efficiently parametrized data-generating process rather than an overparametrized model. This encourages us to consider testing for causal ordering on the basis of multivariate ARMA models, fitted to the available data using the procedures described in the previous section of this paper. Let (Xi,, Xst ) be a pair of time series, generated by a multivariate ARMA(p, 4) model, which we now write as
[
h@)
@12(B)Xl,
@21(B)
@22(B)
11 X2t 1 =
40) [ fl2,u-Q
42(B) &2(B)
I[1 E1t
E2t
’
where the s@~~(B)and Bij(B) are scalar polynomials of respective degrees p
344
PAUL NEWBOLD AND STEVEN M. HOTOPP
and q in the back-shift operator B. Kang [34] has shown that a necessary and
sufficient condition for X, not to cause X,, in the Granger sense, is
Similarly, Xi does not Granger-cause X, if and only if
Notice, from (3.1), that the conditions
+,,(B) = 0,
MB) = 0
(3.2)
are, in general, sufficient but not necessary for noncausality. However, in the special case where the generating mechanism is either pure autoregressive or pure moving average, (3.1) and (3.2) are equivalent. At this point, our analysis diverges from that of Kang. We test directly for the restrictions implied by (3.1) on the parameters of the multivariate ARMA model. In effect this equation implies p + q restrictions on these parameters. One possible approach is to reestimate the model subject to these constraints, and apply a likelihood-ratio test. However, except in the special cases of pure autoregressive or pure moving-average models, our estimation program is not readily adapted to incorporate such constraints. An alternative approach, which is in any case computationally less expensive, is to employ the Wald test. This requires only information obtained from estimation of the unconstrained model. Let R denote the vector of autoregressive and moving-average parameters, and
f(P)=0
(3.3)
the p + q constraints on these parameters implied by (3.1). Further, let b denote the unconstrained maximum-likelihood estimators of R, with estimated covariance matrix Ea. With af/c?#l’ denoting the matrix of first partial derivatives from the left-hand side of (3.3), set
e,=(-$)&(-$)‘. Then, the Wald test statistic for the null hypothesis (3.3) is given by w=f(&‘e,lf@).
(3.4)
Testing Causality
Using Vector ARMA Models
345
Under the null hypothesis, the statistic (3.4) has an asymptotic chi-square distribution with p + 9 degrees of freedom, that hypothesis being rejected for high values of the statistic. We now apply this approach to the data on employment and wages. Writing the ARMA(l, 1) model so that
- h2B
I- h,B
l-&B
- %lB
Xl,
1 - 8,,B
- e,2B
I[ Xzt 1 =
- 8,,B
1 - e,,B
[
the necessary and sufficient conditions (3.1) for X, Granger sense, X, are -
I[1’ E1t Ezt
not to cause, in the
&B(l - %.-$)+ 4,B(l- +zzB)= 0
4, - h2 = 0,
h2e22- e12+,,= 0.
This pair of restrictions can be tested using the Wald statistic (3.4), which employs the estimated covariance matrix of the parameter estimators. This matrix is computed by our estimation program. For these data the computed value of the test statistic was 4.718. Comparison with tabulated values of the chi-square distribution for two degrees of freedom reveals that the null hypothesis that wages do not cause employment can be rejected at the lo-percent but not at the Spercent significance level. We further found that the null hypothesis that employment does not cause wages can be rejected at the 5-percent level. The computed test statistic (3.4) was 7.204. This research was partially funded under contract number SES-8011175.
by the National
Science Foundation,
REFERENCES 1 2
3 4
H. Akaike, Fitting autoregressive models for prediction, Ann. Inst. Statist. Math. 21:243-247 (1969). H. Akaike, On entropy maximization principle, in Applications of Statistics (P. R. Krishnaiah, Ed.), North Holland, Amsterdam, 1977, pp. 27-41. C. F. Ansiey and R. Kohn, Exact likelihood of vector autoregressive-moving average process with missing or aggregated data, Biometiku 70:275-278 (1983). C. F. Ansley and P. Newbold, Multivariate partial autocorrelations, Proc. Bus. Econom. Statist., Amer. Statist. Assoc., 1979, pp. 349-353.
346
8 9 10 11 12
13 14 15
16 17 18 19 20 21 22 23 24 25 26
PAUL NEWBOLD
AND STEVEN M. HOTOPP
C. F. Ansley and P. Newbold, Finite sample properties of estimators for autoregressive moving average models, J. Econometrics 13:159-183 (1980). M. S. Bartlett, On the theoretical specification of sampling properties of autocorrelated time series, I. Roy. Statist. Sot. Ser. I? 8:27-41 (1946). J. M. Beguin, C. Gourieroux, and A. Monfort, Identification of a mixed autoregressive-moving average process: The corner method, Time Series (0. D. Anderson, Ed.), North Holland, Amsterdam, 1980, pp. 423436. G. E. P. Box and G. M. Jenkins, Time Series Analysis, Forecasting and Control, Holden Day, San Francisco, 1970. W. Y. T. Chan and K. F. WaIlis, Multiple time series modelling: Another look at the mink-muskrat interaction, Appl. Statist. 27: 168- 175 (1978). D. M. Cooper and E. F. Wood, Identifying multivariate time series models, J. Time Ser. Anal. 3:153-164 (1982). D. R. Cox and D. V. Hinkley, Theoretical Statistics, Chapman and Hall, London, 1974. N. Davies and J. D. Petruccelli, On the use of the general partial autocorrelation function for order determination in ARMA(p, 9) processes, J. Amer. Statist. Assoc. 79:374-377 (1984). J. Durbin, The fitting of time series models, Reu. Inst. Internat. Statist. 28:233-244 (1960). P. T. Geary and J. Kennan, The employment-real wage relationship: An intemational study, J. Polit. Econ. 90:85&871 (1982). J. Geweke, R. Meese, and W. Dent, Comparing alternative tests of causality in temporal systems: Analytic results and experimental evidence, J. Econometrics 21:161-194 (1983). C. A. Glasbey, A generalization of partial autocorrelations useful in identifying ARMA models, Technometrics 241223-228 (1982). L. G. Godfrey, Testing the adequacy of a time series model, Biometrika 66:67-72 (1979). C. W. J. Granger, Investigating causal relations by econometric models and cross-spectral methods, Econometrica 37:424-438 (1969). C. W. J. Granger, Testing for causality: A personal viewpoint, J. Econom. Dynamics Control 2:329-352 (1980). C. W. J. Granger and M. J. Morris, Time series modelling and interpretation, J. Roy. Statist. Sot. Ser. A 139:246-257 (1976). C. W. J. Granger and P. Newbold, Forecasting Economic Time Series, Academic, New York, 1977. H. L. Gray, G. D. KelIey, D. D. McIntire, A new approach to ARMA modelling, Comm. Statist. B-Simulation Comput. 7:1-77 (1978). D. M. Grether and M. Nerlove, Some properties of “optimal” seasonal adjustment, Econumetrica 38:682-703 (1970). E. J. Hannan, Multiple Time Series, Wiley, New York, 1970. E. J. Hannan and L. Kavalieris, A method for autoregressive-moving average estimation, Biometrika 71:273-280 (1984). E. J. Hannan and J. Rissanen, Recursive estimation of mixed autoregressive-moving average order, Biometrika 69:81-94 (1982); Correction, 70:303.
Testing Causality Using Vector ARMA Models 27 28 29 30 31 32 33 34 35 36
37
38 39 40 41 42 43 44
45 46 47 48
347
S. C. Hilhner and G. C. Tiao, Likelihood function of stationary multiple autoregressive moving average models, J. Amer. Statist. Assoc. 74:652-660 (1979). J. R. M. Hosking, Lagrange multiplier tests of time series models, J. Roy. Statist. Sot. Ser. B 42:170-181 (1980). J. Il. M. Hosking, The multivariate portmanteau statistic, J. Amer. Statist. Assoc. 75:602-608 (1980). J. R. M. Hosking, Lagrange multiplier tests of multivariate time series models, J. Roy. Statist. Sot. Ser. B 43:219-230 (1981). J. R. M. Hosking, Equivalent forms of the multivariate portmanteau statistic, J. Roy. Statist. Sot. Ser. B 43:261-262 (1981). S. M. Hotopp, Practical methods in multivariate time series analysis, Ph.D. Thesis, Dept. of Economics, Univ. of Illinois, Champaign-Urbana, 1985. G. M. Jenkins and A. S. Alavi, Some aspects of modelling and forecasting multivariate time series, J. Time Ser. Anal. 2:1-47 (1981). H. Kang, Necessary and sufficient conditions for causality testing in multivariate ARMA models, J. Time Ser. AnaE. 2:95-101 (1981). J. Ledolter, The analysis of multivariate time series applied to problems in hydrology, J. Hydrol. 36:327-352 (1978). W. K. Li and A. I. McLeod, Distribution of the residual autocorrelations in multivariate ARMA time series models. J. Roy. Statist. Sot. Ser. B 43:231-239 (1981). C. R. Nelson and G. W. Schwert, Tests for predictive relationships between time series variables: A Monte Carlo investigation, J. Amer. Statist. Assoc. 77:11-18 (1982). P. Newbold, The equivalence of two tests for time series model adequacy, Biometrika 67:463-465 (1980). P. Newbold and T. Bos, On q-conditioned partial correlations, J. Time Ser. Anal. 45%55 (1983). D. F. Nicholls and A. D. HalI, The exact likelihood function of multivariate autoregressive-moving average models, Biumetriku, 66:259-264 (1979). D. S. Poskitt and A. R. Tremayne, Testing the specification of a fitted ARMA model, Biumetrika 67:359-363 (1980). D. S. Poskitt and A. R. Tremayne, Diagnostic tests for multiple time series models, Ann. Statist. 10:114-120 (1982). M. H. QuenouiIle, The Analysis of Multiple Time Series, Griffin, London, 1957. C. R. Rao, Large sample tests of statistical hypotheses concerning several parameters with applications to problems of estimation, Proc. Cambridge Philos. Sot. 44:50-59 (1948). J. Rissanen, Modeling by shortest data description, Autumatica 14:465-471 (1978). G. Schwarz, Estimating the dimension of a model, Ann. Statist. 6:461-464 (1978). S. D. Silvey, The Lagrangian multiplier test, Ann. Math. Statist. 30:389-407 (1959). C. A. Sims, Money, income and causality, Am. Econom. Rev. 62:540-552 (1972).
348
PAUL NEWBOLD
49
G. C. Tiao and G. E. P. Box, Modeling multiple time series with applications, J. Amer. Statist. Assoc. 76802-816 (1981). G. C. Tiao, G. E. P. Box, M. R. Grupe, G. B. Hudak, W. R. Bell, and I. Chang, The Wisconsin multiple time series (WMTS-1) program: A preliminary guide, Dept. of Statistics, Univ. of Wisconsin, Madison, 1979. G. C. Tiao and R. S. Tsay, Multiple time series modeling and extended sample cross-correlations, J. Bus. and Ecunon. Statist. 1:43-56 (1983). K. F. WaIIis, Multiple time series analysis and the final form of econometric models, Eccmumetrica 45: 1481- 1497 (1977). G. T. Wilson, The estimation of parameters in multivariate time series models, J. Roy. Statist. Sot. Ser. B 35:76-85 (1973). W. A. Woodward and H. L. Gray, On the relationship between the S array and the Box-Jenkins method of ARMA model identification, J. Amer. Statist. Assoc. 76:579-587 (1981).
59
51 52 53 54
AND STEVEN M. HOTOPP