- Email: [email protected]
Causality in Hydrologic Systems
Copyright © IFAC, Water and Related Land Resource Systems Cleveland, Ohio 1980
CAUSALITY IN HYDROLOGIC SYSTEMS A. Ramachandra Rao and R. L. Kashyap Electn'cal Eng£neering, Engineering, Purdue University, Schools of Civil and Electrical West L afayette, IN IN 47907, USA Lafayette,
Abstract. Cross correlation techniques have been used in the past to detect causal connections in hydrologic systems. However, cross correlation coefficients are weak indicators of causal dependence of time series. Often, spurious but statistically significant cross correlation coefficient values result even for series which are not causally related. In other instances, even though the series are weakly related, cross correlation coefficients indicate stronger relationships. If individual time series whose cross correlation characteristics charact eristics are being analyzed are non-white, these problems become quite severe. Many hydrologic time series are highly correlated, especially for smaller time scales. scales . Consequently better methods are needed to investigate the causal connections between autocorrelated time series such as those commonly encountered in hydrology. One of the recently proposed procedures to investigate the causal connections between pairs of variables whose time series are available is to fit valid models to these time series and to obtain residuals of these models. The cross correlation coefficients between these residuals are tested for statistical significance and for inferences about causal connections between the series . An alternative method is to develop transfer function models of time series. models . Obviously, the transfer function the residuals and examine these models. models of the residuals can be used to generate transfer function models of original origina l time series also. The basic objective of the paper is to investigate these methods. Some of these methods have been investigated by using weekly and monthly rainfall and runoff series from small and intermediate sized watersheds in New ew Zealand. The rainfall and runoff series used were strongly autocorrelated especially especia l ly at smaller lags and were also weakly periodic. The analysis of these data indicated that the weekly and monthly rainfall-runoff series were strongly but "instantaneously" causally connected. Implications of this finding to modeling rainfall-runoff processes have been discussed. Keywords.
I.
Correlation methods, discrete systems, linear systems, models, water resources.
IINTRODUCTION TRODUCTIO between which is represented by PxyCk) PxX(k) where pxCk), PyCk) k is the lag. Let Px(k), Py(k) ana and PxyCk) Pxy(k ) be the correlation functions of x, y and the cross correlation between x and y. Let rxCk), rx(k), ryCk), ry(k ) , rxy(k) be estimates of px(k), px(k ) , pyCk) py(k) and Pxy(kj. Pxy(k). Let N be the number of available observations of x(t) and yet). The asymptotic covariance between cross correlation estimates rxy(k) and rxy(k+£), rxy(k+ £), under the assumptions x( · ) and YC·) y( . ) are covariance stationary that xC-) and normally d~stributed d~stributed has been derived b by Bartlett (1966) as shown in eq. 1.
Cross correlation coefficients and cross spectra have been used to investigate causal relationships between time series. In hydrology, relationships between rainfall and runoff, between sunspots and precipitation, and other pairs of variables have been investigated by using cross correlation analysis. (RodriguezIturbe and Yevjevich, 1968; Kisiel, 1969) Interpretation of significance of cross correlation coefficients is difficult when the individual series which are used to compute the cross correlation coefficients are autoauto correlated. correlated . This aspect can be best discussed by considering the covariance of cross correlation functions. functions . Let us consider two series XCt) x(t) and yet), the cross correlation function
cov{r (k), r (k+ (k+£)} £)} = xy xy 00
_k)-l (N_k)-l 219
IL
.._ J=-OO J--OO
[p (j)p Cj)p (j+£) (j+ £) + P C-j)p (-j)p Cj+2k+£) (j+2k+£) x yY xy xy
A. R. Rao and R. L. Kashyap
220
l y2(j)} l x22 (j) + ~2l (k+t){p2 (j) + ~2l + p (k)p (k+,q,){p2 y2(J.)} xy xy xy "::;i>2
"::;i>2
(j+k+,q,) + P p C-j)p (j+k+,q,)} (-j)p (j+k+t)} - P (k){p (j)p (j+k+t) xy x xy xy y PxyC-j)PyCj+k)}] - Pxy(k+t){Px(j)p Pxy Ck+i){p x (j)p xy (j+k) + P xy (-j)py (j+k)}] (1) Cl)
It is clear from eq. 1 that the covariance rxy(k+,q,) depends on px(·), px(o), between rxy(k) and rxy(k+t) Py(~) Py(~) . and Pxy(·) PxyCo) . . Consequent1y Consequently testing emp1r1ca1 correlat1on coefficients for staemp1r1cal tistical significance and hence for causal connections or other relationships is difficult. Even if x(·) xCo) and y(.) yCo) are independent of each other but are autocorrelated, autocorre1ated, the cross correlation coefficients are correlated as shown in eq. 2. 00
cov{r
xy
(k)r
xy
(N_k)-l (N - k) - 1
(k+1)} (k+i)}
Cj+i) L p (j)p (j) p (j + 1) y j=-oo x Y (2)
In this case, the covariance between rxy(k) rxyCk) and rrxyCk+i) PxCo) and Py(·) Py(o) and xy (k+1) depend on Px(·) hence the significance of individual cross correlation coefficients cannot be estimated xCt) and Yet) yet) are both easily. However, if x(t) white noise series then the asymptotic variance of the autocorrelation coefficients is given by eq. 3. var{r vadr
xy
(k)r (k) r
xy
(k+1)} (k+i)}
= =
(N _ k)-1 k)-l
(3)
Equation 3 can be used to test the significance of cross correlation coefficients. Details of such tests are found in Jenkins and Watts (1968). Jenkins and Watts (1968) have also presented an example in which spuriously large cross correlations exist between two unrelated first order processes which points to the difficulties in interpretation of significance of cross correlation coefficients. In view of these considerations better methods are needed to investigate relationships between hydrologic variables such as the causal connections between hydrologic time series. In order to investigate the causal connections between two variables, we use the uncorrelated residuals from separate valid models fitted to the two corresponding time series. Cross correlation coefficients of the two residual sequences are tested for statistical significance. If these cross correlation coefficients are not significant then the null hypothesis that the two series are causally related is rejected. The main objective of the present paper is to determine the nature and strength of causal relationships between rainfall and runoff time series by testing the corresponding residual cross correlation relationships.
Several tests whi ch h have been Kashyap and Rao (1976) and by are used to investigate these ships.
proposed by Haugh (1976) causal relation-
The paper is organized as follows. In the nex~ section the data used in the study, nex~ the1r cross correlation properties are discussed. Models fitted to the data and validation tests are discussed in Section Ill. The cross correlation characteristics of the residual sequences and the tests on these cross correlation coefficients are discussed next. A discussion of results and their implication to development of models of hydrologic time series are presented in the last section. 11. (a)
DATA USED IN THE STUDY AND THEIR CHARACTERISTICS
Data Used in the Study and Some of Their Elementary Characteristics
Weekly and monthly runoff data from seven watersheds, four of which are in North Island and three in the South Island of New Zealand watershed~ were analyzed. Locations of these watershed~ are shown in Fig. 1 and they are listed in Table 1. The watershed areas range from 222 to 6350 km 22 . The data durations range from thirteen to twenty-eight years. These seven watersheds are in very different climatic regions. These range from rain shadow regions (Wairau) to very wet watersheds (Buller). The geomorphology and geology of these watersheds are also considerably different from each other. The terrain of these watersheds ranges from relatively flat to mountainous, and the channel slope ranges from 0.0046 to 0.014. The soil types and vegetation in these watersheds also vary considerably. The pumice soils in Rangitaiki and Whirinaki watersheds have considerably higher infiltration rates and storage than soils in other watersheds. Runoff from these watersheds and hence the statistical characteristics of runo~f a:e affected by this factor. The vegetat10n 1n these watersheds is mostly varieties tIon of the native grass and forest. Flows in these streams are not controlled. The flow records are rated good or better. better . Weekly and monthly rainfall from ·rain rain gauges either within the watershed or closest to it were also analyzed. The names and the New Zealand Meteorological Service numbers of these rain gauges are given in Table 1 and the location of these rain gauges are shown in Fig. 1. The mean, standard deviation, skewness coefficient and kurtosis of the weekly and monthly data were examined. The rainfall standard deviations were quite large indicating considerable variability in monthly and weekly rainfall. Similar variability was evident in runoff data also, except for runoff from
Causality in Hydrologic Systems
TABLE 1. RIVER NAME
STATION NUMBER
221
Details about Watersheds, Runoff Record Lengths and Climates SITE NAME
RAINFALL STATION
CATCHMENT AREA (km 22 ) Ckm
PERIOD OF RECORD from
CLIMATE
to
Awanui
School Cut
Kaitaia Aero
222
Jan 58
31 Dec 75
Subtropical, average annual precipitation approximately 1600 mm.
15408
Rangitaiki
Murupara
Waimihia Forest
1184
1 Jan 49
31 Dec 75
Average annual precipitation about 1500 mm. Occasional snow in winter.
15410
Whirinaki
Galatea
Minginui Forest
534
Jan 53
31 Dec 75
Average annual precipitation about 1600 mm.
32732
Moawhango
Waiouru
Waiouru
285
Jan 61
31 Dec 75
Average annual precipitation about 1800 mm. Kaimanawa Ranges give rain shadow effect.
60114
Wairau
Dip Flat
Hanmer
505
1 Jan 52
31 Dec 75
Average annual precipitation about 1200 mm with considerable areal variability. In rain shadow Alps . Signifof Southern Alps. icant winter snowfall.
93203
Buller
Te Kuha
Reef ton Reefton
6350
1 Jan 64
31 Dec 75
Average annual precipitation about 2700 mm mainly orographic. Some seasonal snow.
93206
Inangahua
Landing
Reefton Reef ton
100
1 Jan 64
31 Dec 75
Average annual precipitation about 3600 mm, mainly orographic. Some seasonal snow.
1316
Rangitaiki and Whirinaki watersheds for which small . the standard deviations were quite small. All the data were distributed with large skewness. ness . The skewness coefficients of weekly runoff were larger than those of monthly monthl y runoff.
The estimated cross correlation function be{x(t)} and {yet)} is given tween the series {xCt)} by Eq. 7. (k) (k) = =
rr
xy (b)
Cross Correlation Properties of the Data
The cross correlation coefficients of rainfall and runoff data were computed and examined.. The estimates CxyCk) ined Cxy(k) of cross covariYxy(k) between {xCt)} {x(t)} and {Yet)} {yet)} is ance YxyCk) given by Eqs. 4 and 5. N-k C
xy
(k)=N\ (k) = N\
-
L [x(t)-x] [y(t+k)-y]
(4)
t=l
O
N-k C (k) = Nlk [y (t) -y] [x (t+k) -x] N\ L [y(t)-y][x(t+k)-x] yx - t=l
CS) (5)
The lag k varies between zero and kmax .' It is easy to verify that Cxy(O) = Cyx(O) and that the relationships given by Eqs. 6 are valid. = C (-k) Cyx (k) (k)=C xy
(6)
C
x..:-y xy
(k)
Ic xx (0) IC
C C
yy
_ (0)
C (k) ~ -~ (J (J
(7)
°xOy x y
The absolute value of the cross correlation function is less than or equal to one, and the function obeys the relationship given in Eq. 88.. r
xy
=r (k) =
yx
(-k) C -k)
(8)
The cross correlation function rxyCk) rxy(k) at lag k indicates the relationship between present {y(e)} {x(e)} and ryx(k) the oppo{y(o)} and past {x(o)} site. In causally related variables, if yet) x(t) the variable is the effect variable and xCt) indicating cause, then rxyCk) rxy(k) should be larger than rryx(k). xCk). Furthermore, rxy(k) significant . The should be statistically significant. estimator of cross correlation coefficient in Eq. 4 is consistant and asymptotically normal under some generally acceptable conditions on xCt) and yet) (Hannon, 1970). x(t) Examples of cross correlation coefficients and the corresponding two-standard error limits are shown in Fig. 2 for weekly and
222
A. R. Rao and R. L. Kashyap
monthly data. In these figures the correlation coefficients computed with runoff leading rainfall [Cxy(k), Eq. 4] are shown with positive lags and those with rainfall leading runoff [Cxy(k), Eq. 5] are shown with negative A rough measure of the "memory" in lags. A watershed systems is the time, starting from lag zero, during which the cross correlation coefficients corresponding to positive lag are significantly different from zero. zero . The memory of the watershed system indicates the effect of present rainfall on future stream flows. The watershed memory is related to watershed storage. If the watershed storage is large, then the watershed system memory is large and vice versa. The cross correlation coefficients of weekly data (Fig. 2) indicate that large storages are present in Rangitaiki and Whirinaki watersheds in comparison with the storage in other watersheds. The memory of Rangitaiki watershed is about 25 weeks, and that of Whirinaki watershed about 10 weeks. The memory of other watersheds is of the order of four to weeks . A weak seasonal effect is also five weeks. apparent in the cross correlation plots. The lag-zero cross correlation coefficient--also called "instantaneous" correlation coefficient in the statistical literature--was quite large in Awanui, Whirinaki and Buller watersheds, weak in Rangitaiki Rangi taiki and Moawhango ~loawhango watersheds, and is smallest in Wairau watershed. As the "instantaneous" correlation coefficient indicates the rate at which watersheds respond to rainfall input, we may conclude that the Awanui, Buller and Whirinaki watersheds respond much faster to rainfall input than other watersheds. The monthly cross correlation coefficients (Fig. 2) behave similar to the weekly week ly correlation coefficients with the exception that the seasonal pattern may be interpreted more easily. Ill.
parameters and the residual sequence from which is white is selected. A discussion of the conditional maximum likelihood method and whiteness tests on residuals are found nqt in Kashyap and Rao (1976) and are nQt repeated here. Examples of the final models for weekly and monthly data are given in Table 2. The annual and semi-annual periodicities in the data were represented by sine and cosine terms. Terms corresponding to annual and six-month periodicities are generally needed to characterize the changing mean in both weekly and monthly models.
mode led The monthly rainfall series can be modeled by very simple models and the order of the autoregressive part did not exceed two. The monthly flow models are, in general, more complex than rainfall models, especially for Rangitaiki and Whirinaki rivers. The complex correlation structure of monthly flows in these rivers requires terms such as y(k-9) to be inCluded included in order to obta~n obta~n valid models. The structure of weekly and monthly models are approximately analogous to each other. Terms such as yy(k (k-5), -5), Y(k-7), y(k-7), etc., are sometimes needed in the weekly models to account for the strong correlations present in the data at lags corresponding to these terms, but the physical significance of these terms is not easy to explain. Complex storage mechanisms present in Rangitaiki and Whirinaki watersheds are responsible for at least some of these strong correlations and hence the presence of many of these terms in these models. Weekly rainfall models are more complex than runoff models. The residuals from the models such as those given in Table 2 were tested and found to be white . Examples of correlograms (with two white. standard error limits) and cumulative periodograms (with 95% Kolmogorov-Smirnov Limits) of residuals are shown in Fig. Fig . 3.
MODELS FITTED TO THE DATA
Models were fitted to the rainfall {x(t)} and runoff {Yet)} {yet)} data, t = 1, 2, ... , N. Generalized autoregressive models such as those shown shOlYll in eq. 9 were fitted fit ted to the data. data .
The residuals from these models were cross correlated and tested to determine the nature of relationship between rainfall and runoff as discussed in the next section. IV.
+
t-l ) 82211J1j!2( 2 (t-l)
+ ••• +
8£, 1lJ£ 1j!£, (t-l) (t-l ) 8£ 2
+
wet)
2
(9)
In eq. 9 1j!j( llJj(t-l) t-l ) are the deterministic trend terms which are introduced to account for periodicities commonly present in hydrologic hyd rologic data. The parameters in eq. 9 were estimated by the conditional maximum likelihood method after fixing the number of terms m and £2. £'2 ' The residual sequences are tested for whiteness. The model ~with ith the minimum number of
TESTS 0ON RESIDUALS
The uncorrelated residual sequences from models fitted to the "cause" x(t) and "effect" yet) ye t) variables were tested to determine the causal connections between x(t) and yet). y(tj. Let w1(t) be the residual sequence from a validated model fitted to runoff (effect) (e ffect) series and W2(t) the corresponding sequence from a validated model fitted to rainfall (cause) series. If x(t) x (t) causes the effect Yet) ye t ) then there exists at least one j and a corresponding aj ~ 0 so as to satisfy eq. 10. w1 (t) = a.w 2 (t-j) J
+
n(t), j
~
1
(10)
Systems Causality in Hydrologic Syste~
TABLE 2.
223
Examples of Models Fitted to Weekly and Monthly Rainfall and Runoff Data WEEKLY DATA
x(t) == 0.464 + 0.087x(t-l) - 0.014x(t-2) (0.002) (0.002) (0.002)
REEFTON: REEFTO
0.089x(t-5) - 0.085 sinA 33 (t-l) (t-l ) (0.002) (0.001)
+
- 0.0712 cosA 33 (t-l) - 0.06 sinA 44 (t-l) - 0.073 cosA 44 (t-l) (0.001) (0.001) (0.001) BULLER:
yet) Yet)
0 . 33 + 0.343y(t-l) 0.33 (0.001) (0.002)
+
0.017y(t-2) (0.002)
+
+
w2 (t)
0.049y(t-4) - 0.09 sinA 33 (t-l) (0.002) (0.001)
- 0.045 cosA 33 (t-l) - 0.065 sinA 44 (t-l) (0.001) (0.001)
0 . 065 cosA 44 (t-l) 0.065 (0.001)
+
+
wj(t) w1(t)
MONTHLY DATA MINGINUI:
= 0.389
x(t)
+
(0.002) +
0.135x(t-l) - 0.059x(t-2) - 0.032 sinA 2 (t-l) - 0.022 cosA2(t-l) (0 . 001) (0.004) (0.004) (0.001) (0.001)
w2(t) W2(t)
yet) == 0.176 + 0.45y(t-l) - 0.065y(t-2) 0 . 065y(t-2) Yet) (0 . 004) (0.002) (0.004) (0.004)
WHIRINAKI:
+
0.063y(t-7) (0.004)
+
0.057y(t-9) (0.003)
0.087y(t-12) - 0.126y(t-13) - 0.066 sinA1(t-l) sinAj(t-l) - 0.023 COSA1(t-l) cosAj(t-l) - O.087y(t-12) (0.001) (0.003) (0.003) (0.001) +
, 1\2
0.017 sinA 2 (t-l) (0.001)
+
0.0083 cosA 2 (t-l) (0.001)
+
wj(t) w1(t)
_ 21T 6
-
Numbers in parentheses are the standard errors of parameters. ~(,) is a zero mean lID sequence. In eq. 10, nee) If the coefficients aj in eq. 10 can be shown to be significantly dIfferent from zero for anyone value of j, j ~ 1 then we can say that x(t) is the cause and yet) is the effect. In order to test whether aj is significantly different from zero, we may use one of the three tests given below (Kashyap and Rao, 1976). 1976) .
Test 1 The modified likelihood values Ji' i = = 1, 2 for class 1 (aj == 0) and class 2 (aj 1 F 0) are computed by uSIng residuals w W (t) (t ) as shown in eqs. 11.
~ rI IL W~(t~ w~ (t~~ lB
2
=-
"2
In
rI L ~ (t)(I - r lE t=l W
Test 2
pr operty that under und er The test is based on the property the null hypothesis aj = 0, the statistic K(j) given in (12) obeys the F distribution distributi on with parameters 1 and v. v. kk(j) (j) = (N j -1) r ~ 2 (j ) // [[1-r~ 1- r ~ 22((jj)1 )] ~- F ( :-J --j-l)r~2(j) F(( 1l,, v)
((12) 12 )
Choose a threshold Fa from the th e tables tabl e s of F distribution corresponding to the th e prespecipre sp ec ified level of significance a. Then the test t est is as follows: follo ws:
=;'Reject null hypothheesi kk(j)~Fa (j ) ~F a ~ Re j ect tthe he nu 11 hypot si s
= - In lnrI J 1j =-lli t= 1 2 J
is accepted iiff J 1j > J 22'.
k(j)
12
(j
~ j
- I
( 11 ) (11)
We may emphasize that null hypothesis hypothesi s here her e means that rainfall and runoff are not not causally causall y related at lag j. Test 3
In eq. 11, N is the number of residuals w1(t), wj(t), r 12 j2 (j) is the cross correlation coefficient between Wl(t) Wj(t) and W2(t) W2(t ) at lag j with Wl(e) Wj(') leading w2 (e). ('). The cross correlation coefficients r12(j) r j 2(j) and r21(j) r2j(j) are computed by using eqs. 4-7. The hypothesis that that aj == 0
Based on the fact that the set of lagged l agged cross correlations between univariate uni variate residual r es i dual series, which are obtained by fitting fittin g ARMA AR}~ models to two independent series, seri e s, follow f ollow a normal distribution, Haugh (1976) has
224
A. R. Rao and R. L. Kashyap
indepen developed a chi-square test for the independence of the series. series . This is given below as test 3. The test statistic computed in this case is S(j) which is given in eq. eq . 13 where N is the number of observations S(j) S (j)
N
!i
r~ (k) r~2(k)
(13)
2
k=-j
The statistic S(j) is compared to a critical value from the chi-square distribution with (2j+l) degrees of freedom to obtain a significance test for independence of the two series. If Ikl is large compared to N, a different statistic S*(j) is used
S* (j) 5*(j)
= N
2
!
k=-j
(N-(N
IIkl)-lr~2(k) I - 1 r 2I2 (k) k)
(14 ) (14)
Haugh (1976) recommends using S*(j) instead of S(j) for j > N/lO. S*(j) would be compared dis to a critical value from the chi-square distribution with 2j+l degrees of freedom. Further details and properties of S(j) and S*(j) are found in Haugh (1976). Typical cross correlograms of residual sequences are shown in fig. 4. Comparing the correlograms of weekly Waimihia forestforest correlograrns Rangitaiki and Minginui forest-Whirinaki data in fig. 2 with the corresponding cross correlograms of residuals in fig. 3, we see that instead of the first 25 cross correlation coefficients being significant only the first two residual cross correlations are significant. A similar comparison of cross correlaReef ton and Buller tion of monthly data of Reefton with the corresponding residual cross correlation coefficients shows that only the lagzero residual cross correlation coefficient is significant. These conclusions are supported by the results of the three tests mentioned above. By using test 1, the lags at which Uj were not zero were selected. The statistIcs K(j) and S(j) at these lags were estimated by eqs. 12 and 13. Typical K(j) and S(j) values along with the corresponding critical F- and x2-statistics are given in Table 3. The Fand x2-values are given at 5% 5% significance level. The residual cross correlation coefficients and the corresponding two standard error limits are also shown in Table 3. The results presented in Table 3 show that the test 2 is more sensitive than test 3 in the following sense. If the magnitude of a cross correlation coefficient is larger than the corresponding 2-standard error value, the corresponding K(j) statistic is in general greater than the critical F-value. In other 2-test based on the S(j) statiswords, the xX2-test tic is liable to underestimate the significance of cross correlation values in comparison with the F-test.
Another important feature of the results in Table 3 is that the lag-zero cross correlation coefficients are very high compared to others, especially for weekly data. Lag-one weekly cross correlation coefficients are also high for watersheds with large storages. Only the lag-zero and lag-one monthly residual cross correlation coefficients for Whirinaki and Rangitaiki data and lag-zero cross correlation coefficient for Buller data are large. The residual cross correlation coefficients and the results of tests based on them thus indicate that the rainfall and runoff series are instantaneously correlated. Watersheds with large storages have smaller, but significant residual cross correlations at low lags also. The memory in these watersheds is certainly smaller than that indicated by the cross correlograms of rainfall and runoff data. V.
DEVELOPMENT OF TRANSFER FUNCTION MODELS
w1 ( · ) and w2 (·) The residual sequences wI(-) (-) may be used to develop transfer function models between the x and y sequences. In the first stage, models between wI(-) wl( ' ) and w2 (-) (· ) sequences are developed. The initial choice of terms to be included in these models is governed by the results of tests on residuals discussed in the last section. The parameters of these models, because w1l (-) ( · ) and w2(-) W2( ' ) sequences are white, may be estimated by the least squares method. Characteristics of residuals from the models which result are examined for whiteness. The standard errors of parameters are also examined to determine whether the parameters are significant. As a result of these tests models such as those given in Table 4 are derived. It is interesting to note that some of these models require moving average terms as well. By substituting back for w1l (-) (-) ( · ) and w2 (') sequences from the original models, such as those given in Table 3, models relating original variables such as monthly runoff and monthly rainfall may be derived. A monthly model thus derived for the BullerReefton Reef ton data is given below. Models such as those given in eq. 15 are good candidate transfer function models. In deriving eq. 15 some of the terms, the coefficients of y(t+l) = 0.511 + 0.196y(t) - 0.078y(t-l) 0.07Sy(t-l) + 0.165y(t-2) O. 165Y(t - 2) - 0.14y(t-7) O. 14Y(t - 7) + 0.10y(t-14) O. 1Oy (t -14 ) + 0.337x(t+l) 0 . 337x(t+l) + 0.019x(t-l) - 0.309x(t-5) - 0.017x(t-7) 0.017ic(t-7) - 0.081 O.OSl sinAlt sinA1t - 0.069coSA1t O. 069cosA It - 0.1034 sinA 2t - 0.004 0 . 004 COSA COSA 2t + n(t+l)
(15)
Causality in Hydrologic Systems
TABLE 3.
Results from Tests 2 and 3, Cross Correlation Values and the Critical F- and xL-Values xZ-Va1ues (S.E. denotes standard error) WEEKLY DATA 2 S.E. of r (j)
SERIES
LAG j
r (j)
Bu11erReefton Reef ton
0 1 33 51 60 0 1 4 29 51 67 76 82 150 0 1 4 51 SI 67 76 82 112 122 123
0.764 0.172 -0.082 0.101 -0.098 -0 . 098 0.524 0.356 -0.068 -0.066 -0.089 -0.078 0.094 -0.075 0.078 0.430 0.389 0.055 -0.087 -0.057 0.056 -0.061 -0.064 0.074 0.074
0.081
0 6 0 1 0 1 43
0.865 -0.166 0.691 0.187 0.54 0.306 0.153
0.176
WhirinakiMinginui
RangitaikiWaimihia
Bu11erReefton Reef ton WhirinakiMinginui RangitaikiWaimihia
TABLE 4.
225
0.059
0.054
852 18.6 18 . 6 3.86 5.68 5.33 4410 169 5.64 5.00 8.88 6.65 9.62 6.09 6.26 313 246 4.20 10.00 4.32 4.06 4.89 4.38 5.77 5.00
3.84
3.84
3.84
MONTHLY DATA 379 3.41 237 0.124 9.42 0.115 122 30.5 4.93
3.84 3.92 3.84 3.84
2
S (j)
/ (2j+1) X a
N
355.5 18.42 48.12 20.35 49.21 321 148 7 51 47 44 24 17 141 255 210 10 105 42 23 13 20 20 10
3.84 7.42 83 83.54 . 54 59.9 30.1 3.84 7.82 16.91 16 . 91 67.51 61.64 47.37 30.14 22.36 >124 .34 >124.34 3.84 7.82 16.92 118.75 47.37 30.14 22.36 80.23 32.67 7.81
609
96.6 19.4 125.1 43.9 87.1 29.90 84
3.84 12.59 3.84 77.81 .81 3.84 7.82 105.00
Fa (l,N-j) a
K(j)
1168
1381
129 262 300
Models Fitted to Residual Sequences Derived from Runoff-Rainfall Data
STATIONS
MODEL 0. 309w 2(t-5) + n(t+1) (0.003)
Bu11er-Reefton
Wj
wl (t+1)
0.337w2(t+1) 0. 337w 2(t+1) (0.003)
Whirinaki-Minginui
w Wj1 (t+1)
0.103w2 (t) 103w 2(t) 0.360w2(t+1) + 0. (0.001) (0.001)
0.308n(t) + n(t+1) (0.004)
Rangitaiki-Waimihia
Wj
w1 (k+1)
0. 051w 2(t+1) + O. 028w2 (t) (0.0003) (0.0003)
0.OO41w O. 0041 w22 (t-11) (t -11 ) (0.0002)
n(t+1) + nCt+1) W j (.) : w1(e):
Residuals Residua1s from monthly runoff model
e ):: w2 W2C(.)
Residua1s from monthly rainfall model
Numbers in parentheses are the standard errors of parameters.
0.296n(t) 0 . 296n(t) (0.003)
226
A. R. Rao and R. L. Kashyap
which were very small have been ignored. This model may be further refined as desired. Some of the theoretical bases of this type of transfer function model building have been discussed by Haugh and Box (1977). Another important use of models such as those given in Table 4 is in simulation. simulation . As W2(e) w2(-) and nee) sequences in the models given in Table 4 are white noise sequences, they may w1 (-) values be relatively easily simulated. w1(e) thus obtained may be used to generate the sequence of effect variables yet). Thus instead of constructing the usual transfer function models we may use the residual transfer function models to generate sequences of effect variables. Some preliminary work in this direction involving groundwater and precipitation data (Bathala, Rao and Spooner, 1977) has indicated the superiority of this approach over the univariate models. VI.
CONCLUSIONS
On the basis of results presented in this paper we may conclude that the relationships between rainfall and runoff processes may not be as strong as they appear by an analysis of their cross correlation characteristics. The cross correlation characteristics of residuals from validated models fitted to the rainfall and runoff processes and results of tests based on them appear to be better indicators of the causal connections between the processes. These residual sequences may be used to get good candidate transfer function models off the rainfall and runoff processes. ACKNOWLEDGEMENTS Part of A. R. Rao's work reported herein was conducted while he was a Senior Research Fellow of the National Research Advisory Council of New Zealand. I am grateful to the Council for this support.
REFERENCES Bartlett, M. S. (1966). An Introduction to Stochastic Processes, 2nd ed. Cambridge University Press, Cambridge, England. pp. 283-289. Bathala, C. T., A. R. Rao and J. A. Spooner (1977). Applications of Linear Systems Analysis to Groundwater Evaluation Studies. Tech. Rept. No. 91, Water Resources Research Center, Purdue University, W. Lafayette, IN 47907, pp. 69-110. Hannan, E. J. (1970). Multiple Time Series, John Wiley and Sons, Inc., New York, NY. Haugh, L. D. (1976). Checking the Independence of Two Covariance Stationary Time Series: A Uni-Variate Residual Cross-Correlation Approach. J. Am. Stat. Assoc., Vol. 71, 378-385. Haug~ Haug~ D. and G. E. P. Box (1977). Identification of Dynamic Regression (Distributed Lag) Models Connecting Two Time Series. J. Am. Stat. Assoc., Vol. 72, pp. 121-130. Kashyap, R. L. and A. R. Rao (1976). Dynamic Stochastic Models from Empirical Data. Academic Press, New York, pp. 222-226. (1969) . Time Series Analysis Kisiel, C. C. (1969). (Ed.), of Hydrologic Data. In V. T. Chow CEd.), Advances in Hydro Science, Vol. 5, AcaPress . New York, pp. 1-119. demic Press. Rodriguez-Iturbe, I. and V. Yevjevich (1968). The Investigation of Relationship between Hydrologic Time Series and Sunspot Numbers. Colorado State University Hydrology Papers No. 26. Ft. Collins, Colorado.
Causality in Hydrologic Systems
""
227
i
I
I
I
BUtLER AT TEKUHA
I \
~,..,,1.1
(m~JAJCJ;--A_AT_LAN _ _G-.I--I---. .
ANGJTAJKI AT ItNRUPARA
O'}
~
FOIfEST ""'11
WHRINAJ(J
AT GALA TEA
f&10}
FOREST
(BU67f/'
NOMHAN(JO AT MMJ
(31131}
MM~
(Ef&
Fig. 1.
Location of Watersheds
10,.....-------------. 1Or----------------------------,. .
101r---------------------------~ MINGINUI fOREST - WH:Q INAKI [ WJ
WA. IHIHIA FORESTFQRE ST- RANGITAIKI RANGITA IKI[W WAIHIHIA [W
REEF 0 ... RE~FTO
- BULLER !W] BUllE R [loll
~ 06
z ~ ~ ~
~2
0
~
~ -06 ~-06
0 ~
10,00 -80
1~--------------------~
o
~ IM IHI A
--60
10
~O
60
00
10' r---------------------------~ MINGINUI
FORES T- RANGllAl(ilH]
-Ob l AG [ WEEKS)
[ WEEKS]
10)------------------------------
REEF TO'" - BUlLHII'1~
FCREST- WHRINA,l(1
IMl 06 z
g ;:;
01
0 ~
- - --
z
Z 0
et ~
-~
-01
-
--, - - ~
-
-
- - -
--.
~ ~
0
-0 -06 VI
~-O
LAG
MONTHS
u -1 {}-,O:---;';t.O~-3""O--""'20=--""1O,--,O!c---,\1O~2bO,.-,3!"O-I..J,O,......,J"O ~O
Fig. 2.
~
0 '-'
:.;;
! ''''
~
.... S'
.....
-' O_=---..l:-........-,!~--:!_b",....---,\· ~.--+----,L,---="zt~-="3c,....-~ SC
Examples of Cross Correlation Coefficients of Weekly eekly (Top Row) and ~lonthly onthly (Bottom Row) Data
228
A. R. Rao and R. L. Kashyap
10
10r-------.-----..--------, 10,---------------- - - , BULLERIW BULLERI
·Or-------.-----..--------, 1o .-------------------~ RANGITAlKI RAN GI TAlKII[W) WI
r------,----,.-----r----.--_
/
/ /
E
L.::: --- --------.: : -
o -- - ---- - --- -- -
C QI
U
~
C
QI
t"
- - - -- ----.-, - - --0 -\--"7',------
QI
- --- - - --
- --
8'
o
u
'~weeksl
10
"
'() 10
10
/
/
E /
.~
~04
/
'"o
/
~<)~/
"U
~
/
/
~~6
o
Cl
i~
// /
QI
'"
c. 0.
/
/ /
/
/
~~V}/
0·1 02 ~3 03 04 04 Frequency Frequen cy [(PMI I(PMI
Fig. 3.
.
// . / //
/
//
/
...... ~2--003 .....3 --0....... 0 44-....J0.5 05 02 [(py II '~requency requenc y I(PY
01 0·1 0·2 OJ 04 0·4 02 03 Frequency. (CPWJ Frequ ency I(PW}
04 0·4
0-5 OS
'" QI
- - -
--
- /\-
- - - - - -
cC
--
--- - ---- - -----
/
- - - -,:: - - -. 7\ ---
QI ~
~'?-
QI ~
o
- -":' - - - - - ~
/ /
/
/
o
o
u
u LJ
01 0·2 02 0OJ3 Frequency Freq uenc y [( I( PMl PM }
0·4 04
OS 05
-10 . . - - - - - - :fA-----.b:-------.!30 -, 0 L L-----~'hr-O----~20Or-------.1 1O Log [months) Imon th s}
-, 0);0 ------;'-:;-0----~20<;;--------.I -100~--~1'A0---::2~0-----"J30 Log Imonths} (months)
Examples of Correlograms and Cumulative Periodograms of Residuals _..
u: u.: l.L u. w W
o uU
u.i W
Q: 0:: Q: 0:: 0 D-
U
tn Vl tn Vl
oQ:
0::
II
...
OS 05
100 r-------.-----B-U,....LL-E-R-[-H-I---, 1 ,---~--BU~ LL-E-R...,.IM-:I--,
C c
/
/
w
/
// /
/
/
/ / /
/
/
00·66
3" 0·2
/
/
,
,
/
E
/
02
,
// /
/
E
/
d
E
.3
1·0.-------,----..--------, RANGITAIKI IM} rH)
/
~/
~o
/
/
//
/
&
/
/
/
/ //
~\. / //, ~~\9
'"
0
//
o000~---"'0'-1 0 1-
30
BULLER [Ml//
OB 0·8
/
//
~o
/
/
Log [weeks)} Lo g Iweeks
1o.---------,r---R-r- G-I-TA"'T"'-KI-[-MY"T'""-:"/----., AN
0·8
20
/
/
/
/ /
/
/
/
/
LJ W
30
20 20
Log [weeks)
/
02 02
::::l
~
/
//
/
Cj 0, /
/
E E
o
o
u
// / /
/
/ h\o/ ~\o/ / // /
06 04
i:
~
o
LJ
h
//
BULLER (Wl // IW I /// /
/
/
0·8 08 III
10.---------,-----..--..----~10
RA GITAlKI [WI / / RANGI TA IKllwl
D.-
...
.ll
¥
-
... ...
"m
(WEEKS ) LAG (WEEKS)
<~
+
U
-1
u: u. u. l.L
w W
0o
U
l:l
W u.i Q:
0:: Q: 0:: 0D U (J) Vl (J) Vl
o0Q:
+
LAG (WEEKS)
0::: U
I
u: u.: l.L u. W w o u U
W W
--60-------- - -4-------.0.-"A LI.
B_,
(f) Vl (J) Vl
o
Q: 0:: U
6___6_~".
.t._m__-- ~ ~P~~~_ILL
---------------------llllll
""
_ _ A -I>..I:.-~~-- ~~~.~~ ~ __
LAG (MONTHS)
"-A-
--:-;'5
+
Fig. 4. Cross Correlograms of Residuals from Weekly Weekl y and Monthly Monthl y Models: Rangitaiki-Waimihia (weekly, top), Whirinaki-Minginui (weekly, middle), middle ) , and BUller-Reefton Buller-Reefton (monthly, bottom). Positive and negative lags respectively respectiv ely indicate runoff and rainfall as leading variables. vari ables.
CAUSALITY IN HYDROLOGIC SYSTEMS A. Ramachandra Rao and R. L. Kashyap Electn'cal Eng£neering, Engineering, Purdue University, Schools of Civil and Electrical West L afayette, IN IN 47907, USA Lafayette,
Abstract. Cross correlation techniques have been used in the past to detect causal connections in hydrologic systems. However, cross correlation coefficients are weak indicators of causal dependence of time series. Often, spurious but statistically significant cross correlation coefficient values result even for series which are not causally related. In other instances, even though the series are weakly related, cross correlation coefficients indicate stronger relationships. If individual time series whose cross correlation characteristics charact eristics are being analyzed are non-white, these problems become quite severe. Many hydrologic time series are highly correlated, especially for smaller time scales. scales . Consequently better methods are needed to investigate the causal connections between autocorrelated time series such as those commonly encountered in hydrology. One of the recently proposed procedures to investigate the causal connections between pairs of variables whose time series are available is to fit valid models to these time series and to obtain residuals of these models. The cross correlation coefficients between these residuals are tested for statistical significance and for inferences about causal connections between the series . An alternative method is to develop transfer function models of time series. models . Obviously, the transfer function the residuals and examine these models. models of the residuals can be used to generate transfer function models of original origina l time series also. The basic objective of the paper is to investigate these methods. Some of these methods have been investigated by using weekly and monthly rainfall and runoff series from small and intermediate sized watersheds in New ew Zealand. The rainfall and runoff series used were strongly autocorrelated especially especia l ly at smaller lags and were also weakly periodic. The analysis of these data indicated that the weekly and monthly rainfall-runoff series were strongly but "instantaneously" causally connected. Implications of this finding to modeling rainfall-runoff processes have been discussed. Keywords.
I.
Correlation methods, discrete systems, linear systems, models, water resources.
IINTRODUCTION TRODUCTIO between which is represented by PxyCk) PxX(k) where pxCk), PyCk) k is the lag. Let Px(k), Py(k) ana and PxyCk) Pxy(k ) be the correlation functions of x, y and the cross correlation between x and y. Let rxCk), rx(k), ryCk), ry(k ) , rxy(k) be estimates of px(k), px(k ) , pyCk) py(k) and Pxy(kj. Pxy(k). Let N be the number of available observations of x(t) and yet). The asymptotic covariance between cross correlation estimates rxy(k) and rxy(k+£), rxy(k+ £), under the assumptions x( · ) and YC·) y( . ) are covariance stationary that xC-) and normally d~stributed d~stributed has been derived b by Bartlett (1966) as shown in eq. 1.
Cross correlation coefficients and cross spectra have been used to investigate causal relationships between time series. In hydrology, relationships between rainfall and runoff, between sunspots and precipitation, and other pairs of variables have been investigated by using cross correlation analysis. (RodriguezIturbe and Yevjevich, 1968; Kisiel, 1969) Interpretation of significance of cross correlation coefficients is difficult when the individual series which are used to compute the cross correlation coefficients are autoauto correlated. correlated . This aspect can be best discussed by considering the covariance of cross correlation functions. functions . Let us consider two series XCt) x(t) and yet), the cross correlation function
cov{r (k), r (k+ (k+£)} £)} = xy xy 00
_k)-l (N_k)-l 219
IL
.._ J=-OO J--OO
[p (j)p Cj)p (j+£) (j+ £) + P C-j)p (-j)p Cj+2k+£) (j+2k+£) x yY xy xy
A. R. Rao and R. L. Kashyap
220
l y2(j)} l x22 (j) + ~2l (k+t){p2 (j) + ~2l + p (k)p (k+,q,){p2 y2(J.)} xy xy xy "::;i>2
"::;i>2
(j+k+,q,) + P p C-j)p (j+k+,q,)} (-j)p (j+k+t)} - P (k){p (j)p (j+k+t) xy x xy xy y PxyC-j)PyCj+k)}] - Pxy(k+t){Px(j)p Pxy Ck+i){p x (j)p xy (j+k) + P xy (-j)py (j+k)}] (1) Cl)
It is clear from eq. 1 that the covariance rxy(k+,q,) depends on px(·), px(o), between rxy(k) and rxy(k+t) Py(~) Py(~) . and Pxy(·) PxyCo) . . Consequent1y Consequently testing emp1r1ca1 correlat1on coefficients for staemp1r1cal tistical significance and hence for causal connections or other relationships is difficult. Even if x(·) xCo) and y(.) yCo) are independent of each other but are autocorrelated, autocorre1ated, the cross correlation coefficients are correlated as shown in eq. 2. 00
cov{r
xy
(k)r
xy
(N_k)-l (N - k) - 1
(k+1)} (k+i)}
Cj+i) L p (j)p (j) p (j + 1) y j=-oo x Y (2)
In this case, the covariance between rxy(k) rxyCk) and rrxyCk+i) PxCo) and Py(·) Py(o) and xy (k+1) depend on Px(·) hence the significance of individual cross correlation coefficients cannot be estimated xCt) and Yet) yet) are both easily. However, if x(t) white noise series then the asymptotic variance of the autocorrelation coefficients is given by eq. 3. var{r vadr
xy
(k)r (k) r
xy
(k+1)} (k+i)}
= =
(N _ k)-1 k)-l
(3)
Equation 3 can be used to test the significance of cross correlation coefficients. Details of such tests are found in Jenkins and Watts (1968). Jenkins and Watts (1968) have also presented an example in which spuriously large cross correlations exist between two unrelated first order processes which points to the difficulties in interpretation of significance of cross correlation coefficients. In view of these considerations better methods are needed to investigate relationships between hydrologic variables such as the causal connections between hydrologic time series. In order to investigate the causal connections between two variables, we use the uncorrelated residuals from separate valid models fitted to the two corresponding time series. Cross correlation coefficients of the two residual sequences are tested for statistical significance. If these cross correlation coefficients are not significant then the null hypothesis that the two series are causally related is rejected. The main objective of the present paper is to determine the nature and strength of causal relationships between rainfall and runoff time series by testing the corresponding residual cross correlation relationships.
Several tests whi ch h have been Kashyap and Rao (1976) and by are used to investigate these ships.
proposed by Haugh (1976) causal relation-
The paper is organized as follows. In the nex~ section the data used in the study, nex~ the1r cross correlation properties are discussed. Models fitted to the data and validation tests are discussed in Section Ill. The cross correlation characteristics of the residual sequences and the tests on these cross correlation coefficients are discussed next. A discussion of results and their implication to development of models of hydrologic time series are presented in the last section. 11. (a)
DATA USED IN THE STUDY AND THEIR CHARACTERISTICS
Data Used in the Study and Some of Their Elementary Characteristics
Weekly and monthly runoff data from seven watersheds, four of which are in North Island and three in the South Island of New Zealand watershed~ were analyzed. Locations of these watershed~ are shown in Fig. 1 and they are listed in Table 1. The watershed areas range from 222 to 6350 km 22 . The data durations range from thirteen to twenty-eight years. These seven watersheds are in very different climatic regions. These range from rain shadow regions (Wairau) to very wet watersheds (Buller). The geomorphology and geology of these watersheds are also considerably different from each other. The terrain of these watersheds ranges from relatively flat to mountainous, and the channel slope ranges from 0.0046 to 0.014. The soil types and vegetation in these watersheds also vary considerably. The pumice soils in Rangitaiki and Whirinaki watersheds have considerably higher infiltration rates and storage than soils in other watersheds. Runoff from these watersheds and hence the statistical characteristics of runo~f a:e affected by this factor. The vegetat10n 1n these watersheds is mostly varieties tIon of the native grass and forest. Flows in these streams are not controlled. The flow records are rated good or better. better . Weekly and monthly rainfall from ·rain rain gauges either within the watershed or closest to it were also analyzed. The names and the New Zealand Meteorological Service numbers of these rain gauges are given in Table 1 and the location of these rain gauges are shown in Fig. 1. The mean, standard deviation, skewness coefficient and kurtosis of the weekly and monthly data were examined. The rainfall standard deviations were quite large indicating considerable variability in monthly and weekly rainfall. Similar variability was evident in runoff data also, except for runoff from
Causality in Hydrologic Systems
TABLE 1. RIVER NAME
STATION NUMBER
221
Details about Watersheds, Runoff Record Lengths and Climates SITE NAME
RAINFALL STATION
CATCHMENT AREA (km 22 ) Ckm
PERIOD OF RECORD from
CLIMATE
to
Awanui
School Cut
Kaitaia Aero
222
Jan 58
31 Dec 75
Subtropical, average annual precipitation approximately 1600 mm.
15408
Rangitaiki
Murupara
Waimihia Forest
1184
1 Jan 49
31 Dec 75
Average annual precipitation about 1500 mm. Occasional snow in winter.
15410
Whirinaki
Galatea
Minginui Forest
534
Jan 53
31 Dec 75
Average annual precipitation about 1600 mm.
32732
Moawhango
Waiouru
Waiouru
285
Jan 61
31 Dec 75
Average annual precipitation about 1800 mm. Kaimanawa Ranges give rain shadow effect.
60114
Wairau
Dip Flat
Hanmer
505
1 Jan 52
31 Dec 75
Average annual precipitation about 1200 mm with considerable areal variability. In rain shadow Alps . Signifof Southern Alps. icant winter snowfall.
93203
Buller
Te Kuha
Reef ton Reefton
6350
1 Jan 64
31 Dec 75
Average annual precipitation about 2700 mm mainly orographic. Some seasonal snow.
93206
Inangahua
Landing
Reefton Reef ton
100
1 Jan 64
31 Dec 75
Average annual precipitation about 3600 mm, mainly orographic. Some seasonal snow.
1316
Rangitaiki and Whirinaki watersheds for which small . the standard deviations were quite small. All the data were distributed with large skewness. ness . The skewness coefficients of weekly runoff were larger than those of monthly monthl y runoff.
The estimated cross correlation function be{x(t)} and {yet)} is given tween the series {xCt)} by Eq. 7. (k) (k) = =
rr
xy (b)
Cross Correlation Properties of the Data
The cross correlation coefficients of rainfall and runoff data were computed and examined.. The estimates CxyCk) ined Cxy(k) of cross covariYxy(k) between {xCt)} {x(t)} and {Yet)} {yet)} is ance YxyCk) given by Eqs. 4 and 5. N-k C
xy
(k)=N\ (k) = N\
-
L [x(t)-x] [y(t+k)-y]
(4)
t=l
O
N-k C (k) = Nlk [y (t) -y] [x (t+k) -x] N\ L [y(t)-y][x(t+k)-x] yx - t=l
CS) (5)
The lag k varies between zero and kmax .' It is easy to verify that Cxy(O) = Cyx(O) and that the relationships given by Eqs. 6 are valid. = C (-k) Cyx (k) (k)=C xy
(6)
C
x..:-y xy
(k)
Ic xx (0) IC
C C
yy
_ (0)
C (k) ~ -~ (J (J
(7)
°xOy x y
The absolute value of the cross correlation function is less than or equal to one, and the function obeys the relationship given in Eq. 88.. r
xy
=r (k) =
yx
(-k) C -k)
(8)
The cross correlation function rxyCk) rxy(k) at lag k indicates the relationship between present {y(e)} {x(e)} and ryx(k) the oppo{y(o)} and past {x(o)} site. In causally related variables, if yet) x(t) the variable is the effect variable and xCt) indicating cause, then rxyCk) rxy(k) should be larger than rryx(k). xCk). Furthermore, rxy(k) significant . The should be statistically significant. estimator of cross correlation coefficient in Eq. 4 is consistant and asymptotically normal under some generally acceptable conditions on xCt) and yet) (Hannon, 1970). x(t) Examples of cross correlation coefficients and the corresponding two-standard error limits are shown in Fig. 2 for weekly and
222
A. R. Rao and R. L. Kashyap
monthly data. In these figures the correlation coefficients computed with runoff leading rainfall [Cxy(k), Eq. 4] are shown with positive lags and those with rainfall leading runoff [Cxy(k), Eq. 5] are shown with negative A rough measure of the "memory" in lags. A watershed systems is the time, starting from lag zero, during which the cross correlation coefficients corresponding to positive lag are significantly different from zero. zero . The memory of the watershed system indicates the effect of present rainfall on future stream flows. The watershed memory is related to watershed storage. If the watershed storage is large, then the watershed system memory is large and vice versa. The cross correlation coefficients of weekly data (Fig. 2) indicate that large storages are present in Rangitaiki and Whirinaki watersheds in comparison with the storage in other watersheds. The memory of Rangitaiki watershed is about 25 weeks, and that of Whirinaki watershed about 10 weeks. The memory of other watersheds is of the order of four to weeks . A weak seasonal effect is also five weeks. apparent in the cross correlation plots. The lag-zero cross correlation coefficient--also called "instantaneous" correlation coefficient in the statistical literature--was quite large in Awanui, Whirinaki and Buller watersheds, weak in Rangitaiki Rangi taiki and Moawhango ~loawhango watersheds, and is smallest in Wairau watershed. As the "instantaneous" correlation coefficient indicates the rate at which watersheds respond to rainfall input, we may conclude that the Awanui, Buller and Whirinaki watersheds respond much faster to rainfall input than other watersheds. The monthly cross correlation coefficients (Fig. 2) behave similar to the weekly week ly correlation coefficients with the exception that the seasonal pattern may be interpreted more easily. Ill.
parameters and the residual sequence from which is white is selected. A discussion of the conditional maximum likelihood method and whiteness tests on residuals are found nqt in Kashyap and Rao (1976) and are nQt repeated here. Examples of the final models for weekly and monthly data are given in Table 2. The annual and semi-annual periodicities in the data were represented by sine and cosine terms. Terms corresponding to annual and six-month periodicities are generally needed to characterize the changing mean in both weekly and monthly models.
mode led The monthly rainfall series can be modeled by very simple models and the order of the autoregressive part did not exceed two. The monthly flow models are, in general, more complex than rainfall models, especially for Rangitaiki and Whirinaki rivers. The complex correlation structure of monthly flows in these rivers requires terms such as y(k-9) to be inCluded included in order to obta~n obta~n valid models. The structure of weekly and monthly models are approximately analogous to each other. Terms such as yy(k (k-5), -5), Y(k-7), y(k-7), etc., are sometimes needed in the weekly models to account for the strong correlations present in the data at lags corresponding to these terms, but the physical significance of these terms is not easy to explain. Complex storage mechanisms present in Rangitaiki and Whirinaki watersheds are responsible for at least some of these strong correlations and hence the presence of many of these terms in these models. Weekly rainfall models are more complex than runoff models. The residuals from the models such as those given in Table 2 were tested and found to be white . Examples of correlograms (with two white. standard error limits) and cumulative periodograms (with 95% Kolmogorov-Smirnov Limits) of residuals are shown in Fig. Fig . 3.
MODELS FITTED TO THE DATA
Models were fitted to the rainfall {x(t)} and runoff {Yet)} {yet)} data, t = 1, 2, ... , N. Generalized autoregressive models such as those shown shOlYll in eq. 9 were fitted fit ted to the data. data .
The residuals from these models were cross correlated and tested to determine the nature of relationship between rainfall and runoff as discussed in the next section. IV.
+
t-l ) 82211J1j!2( 2 (t-l)
+ ••• +
8£, 1lJ£ 1j!£, (t-l) (t-l ) 8£ 2
+
wet)
2
(9)
In eq. 9 1j!j( llJj(t-l) t-l ) are the deterministic trend terms which are introduced to account for periodicities commonly present in hydrologic hyd rologic data. The parameters in eq. 9 were estimated by the conditional maximum likelihood method after fixing the number of terms m and £2. £'2 ' The residual sequences are tested for whiteness. The model ~with ith the minimum number of
TESTS 0ON RESIDUALS
The uncorrelated residual sequences from models fitted to the "cause" x(t) and "effect" yet) ye t) variables were tested to determine the causal connections between x(t) and yet). y(tj. Let w1(t) be the residual sequence from a validated model fitted to runoff (effect) (e ffect) series and W2(t) the corresponding sequence from a validated model fitted to rainfall (cause) series. If x(t) x (t) causes the effect Yet) ye t ) then there exists at least one j and a corresponding aj ~ 0 so as to satisfy eq. 10. w1 (t) = a.w 2 (t-j) J
+
n(t), j
~
1
(10)
Systems Causality in Hydrologic Syste~
TABLE 2.
223
Examples of Models Fitted to Weekly and Monthly Rainfall and Runoff Data WEEKLY DATA
x(t) == 0.464 + 0.087x(t-l) - 0.014x(t-2) (0.002) (0.002) (0.002)
REEFTON: REEFTO
0.089x(t-5) - 0.085 sinA 33 (t-l) (t-l ) (0.002) (0.001)
+
- 0.0712 cosA 33 (t-l) - 0.06 sinA 44 (t-l) - 0.073 cosA 44 (t-l) (0.001) (0.001) (0.001) BULLER:
yet) Yet)
0 . 33 + 0.343y(t-l) 0.33 (0.001) (0.002)
+
0.017y(t-2) (0.002)
+
+
w2 (t)
0.049y(t-4) - 0.09 sinA 33 (t-l) (0.002) (0.001)
- 0.045 cosA 33 (t-l) - 0.065 sinA 44 (t-l) (0.001) (0.001)
0 . 065 cosA 44 (t-l) 0.065 (0.001)
+
+
wj(t) w1(t)
MONTHLY DATA MINGINUI:
= 0.389
x(t)
+
(0.002) +
0.135x(t-l) - 0.059x(t-2) - 0.032 sinA 2 (t-l) - 0.022 cosA2(t-l) (0 . 001) (0.004) (0.004) (0.001) (0.001)
w2(t) W2(t)
yet) == 0.176 + 0.45y(t-l) - 0.065y(t-2) 0 . 065y(t-2) Yet) (0 . 004) (0.002) (0.004) (0.004)
WHIRINAKI:
+
0.063y(t-7) (0.004)
+
0.057y(t-9) (0.003)
0.087y(t-12) - 0.126y(t-13) - 0.066 sinA1(t-l) sinAj(t-l) - 0.023 COSA1(t-l) cosAj(t-l) - O.087y(t-12) (0.001) (0.003) (0.003) (0.001) +
, 1\2
0.017 sinA 2 (t-l) (0.001)
+
0.0083 cosA 2 (t-l) (0.001)
+
wj(t) w1(t)
_ 21T 6
-
Numbers in parentheses are the standard errors of parameters. ~(,) is a zero mean lID sequence. In eq. 10, nee) If the coefficients aj in eq. 10 can be shown to be significantly dIfferent from zero for anyone value of j, j ~ 1 then we can say that x(t) is the cause and yet) is the effect. In order to test whether aj is significantly different from zero, we may use one of the three tests given below (Kashyap and Rao, 1976). 1976) .
Test 1 The modified likelihood values Ji' i = = 1, 2 for class 1 (aj == 0) and class 2 (aj 1 F 0) are computed by uSIng residuals w W (t) (t ) as shown in eqs. 11.
~ rI IL W~(t~ w~ (t~~ lB
2
=-
"2
In
rI L ~ (t)(I - r lE t=l W
Test 2
pr operty that under und er The test is based on the property the null hypothesis aj = 0, the statistic K(j) given in (12) obeys the F distribution distributi on with parameters 1 and v. v. kk(j) (j) = (N j -1) r ~ 2 (j ) // [[1-r~ 1- r ~ 22((jj)1 )] ~- F ( :-J --j-l)r~2(j) F(( 1l,, v)
((12) 12 )
Choose a threshold Fa from the th e tables tabl e s of F distribution corresponding to the th e prespecipre sp ec ified level of significance a. Then the test t est is as follows: follo ws:
=;'Reject null hypothheesi kk(j)~Fa (j ) ~F a ~ Re j ect tthe he nu 11 hypot si s
= - In lnrI J 1j =-lli t= 1 2 J
is accepted iiff J 1j > J 22'.
k(j)
12
(j
~ j
- I
( 11 ) (11)
We may emphasize that null hypothesis hypothesi s here her e means that rainfall and runoff are not not causally causall y related at lag j. Test 3
In eq. 11, N is the number of residuals w1(t), wj(t), r 12 j2 (j) is the cross correlation coefficient between Wl(t) Wj(t) and W2(t) W2(t ) at lag j with Wl(e) Wj(') leading w2 (e). ('). The cross correlation coefficients r12(j) r j 2(j) and r21(j) r2j(j) are computed by using eqs. 4-7. The hypothesis that that aj == 0
Based on the fact that the set of lagged l agged cross correlations between univariate uni variate residual r es i dual series, which are obtained by fitting fittin g ARMA AR}~ models to two independent series, seri e s, follow f ollow a normal distribution, Haugh (1976) has
224
A. R. Rao and R. L. Kashyap
indepen developed a chi-square test for the independence of the series. series . This is given below as test 3. The test statistic computed in this case is S(j) which is given in eq. eq . 13 where N is the number of observations S(j) S (j)
N
!i
r~ (k) r~2(k)
(13)
2
k=-j
The statistic S(j) is compared to a critical value from the chi-square distribution with (2j+l) degrees of freedom to obtain a significance test for independence of the two series. If Ikl is large compared to N, a different statistic S*(j) is used
S* (j) 5*(j)
= N
2
!
k=-j
(N-(N
IIkl)-lr~2(k) I - 1 r 2I2 (k) k)
(14 ) (14)
Haugh (1976) recommends using S*(j) instead of S(j) for j > N/lO. S*(j) would be compared dis to a critical value from the chi-square distribution with 2j+l degrees of freedom. Further details and properties of S(j) and S*(j) are found in Haugh (1976). Typical cross correlograms of residual sequences are shown in fig. 4. Comparing the correlograms of weekly Waimihia forestforest correlograrns Rangitaiki and Minginui forest-Whirinaki data in fig. 2 with the corresponding cross correlograms of residuals in fig. 3, we see that instead of the first 25 cross correlation coefficients being significant only the first two residual cross correlations are significant. A similar comparison of cross correlaReef ton and Buller tion of monthly data of Reefton with the corresponding residual cross correlation coefficients shows that only the lagzero residual cross correlation coefficient is significant. These conclusions are supported by the results of the three tests mentioned above. By using test 1, the lags at which Uj were not zero were selected. The statistIcs K(j) and S(j) at these lags were estimated by eqs. 12 and 13. Typical K(j) and S(j) values along with the corresponding critical F- and x2-statistics are given in Table 3. The Fand x2-values are given at 5% 5% significance level. The residual cross correlation coefficients and the corresponding two standard error limits are also shown in Table 3. The results presented in Table 3 show that the test 2 is more sensitive than test 3 in the following sense. If the magnitude of a cross correlation coefficient is larger than the corresponding 2-standard error value, the corresponding K(j) statistic is in general greater than the critical F-value. In other 2-test based on the S(j) statiswords, the xX2-test tic is liable to underestimate the significance of cross correlation values in comparison with the F-test.
Another important feature of the results in Table 3 is that the lag-zero cross correlation coefficients are very high compared to others, especially for weekly data. Lag-one weekly cross correlation coefficients are also high for watersheds with large storages. Only the lag-zero and lag-one monthly residual cross correlation coefficients for Whirinaki and Rangitaiki data and lag-zero cross correlation coefficient for Buller data are large. The residual cross correlation coefficients and the results of tests based on them thus indicate that the rainfall and runoff series are instantaneously correlated. Watersheds with large storages have smaller, but significant residual cross correlations at low lags also. The memory in these watersheds is certainly smaller than that indicated by the cross correlograms of rainfall and runoff data. V.
DEVELOPMENT OF TRANSFER FUNCTION MODELS
w1 ( · ) and w2 (·) The residual sequences wI(-) (-) may be used to develop transfer function models between the x and y sequences. In the first stage, models between wI(-) wl( ' ) and w2 (-) (· ) sequences are developed. The initial choice of terms to be included in these models is governed by the results of tests on residuals discussed in the last section. The parameters of these models, because w1l (-) ( · ) and w2(-) W2( ' ) sequences are white, may be estimated by the least squares method. Characteristics of residuals from the models which result are examined for whiteness. The standard errors of parameters are also examined to determine whether the parameters are significant. As a result of these tests models such as those given in Table 4 are derived. It is interesting to note that some of these models require moving average terms as well. By substituting back for w1l (-) (-) ( · ) and w2 (') sequences from the original models, such as those given in Table 3, models relating original variables such as monthly runoff and monthly rainfall may be derived. A monthly model thus derived for the BullerReefton Reef ton data is given below. Models such as those given in eq. 15 are good candidate transfer function models. In deriving eq. 15 some of the terms, the coefficients of y(t+l) = 0.511 + 0.196y(t) - 0.078y(t-l) 0.07Sy(t-l) + 0.165y(t-2) O. 165Y(t - 2) - 0.14y(t-7) O. 14Y(t - 7) + 0.10y(t-14) O. 1Oy (t -14 ) + 0.337x(t+l) 0 . 337x(t+l) + 0.019x(t-l) - 0.309x(t-5) - 0.017x(t-7) 0.017ic(t-7) - 0.081 O.OSl sinAlt sinA1t - 0.069coSA1t O. 069cosA It - 0.1034 sinA 2t - 0.004 0 . 004 COSA COSA 2t + n(t+l)
(15)
Causality in Hydrologic Systems
TABLE 3.
Results from Tests 2 and 3, Cross Correlation Values and the Critical F- and xL-Values xZ-Va1ues (S.E. denotes standard error) WEEKLY DATA 2 S.E. of r (j)
SERIES
LAG j
r (j)
Bu11erReefton Reef ton
0 1 33 51 60 0 1 4 29 51 67 76 82 150 0 1 4 51 SI 67 76 82 112 122 123
0.764 0.172 -0.082 0.101 -0.098 -0 . 098 0.524 0.356 -0.068 -0.066 -0.089 -0.078 0.094 -0.075 0.078 0.430 0.389 0.055 -0.087 -0.057 0.056 -0.061 -0.064 0.074 0.074
0.081
0 6 0 1 0 1 43
0.865 -0.166 0.691 0.187 0.54 0.306 0.153
0.176
WhirinakiMinginui
RangitaikiWaimihia
Bu11erReefton Reef ton WhirinakiMinginui RangitaikiWaimihia
TABLE 4.
225
0.059
0.054
852 18.6 18 . 6 3.86 5.68 5.33 4410 169 5.64 5.00 8.88 6.65 9.62 6.09 6.26 313 246 4.20 10.00 4.32 4.06 4.89 4.38 5.77 5.00
3.84
3.84
3.84
MONTHLY DATA 379 3.41 237 0.124 9.42 0.115 122 30.5 4.93
3.84 3.92 3.84 3.84
2
S (j)
/ (2j+1) X a
N
355.5 18.42 48.12 20.35 49.21 321 148 7 51 47 44 24 17 141 255 210 10 105 42 23 13 20 20 10
3.84 7.42 83 83.54 . 54 59.9 30.1 3.84 7.82 16.91 16 . 91 67.51 61.64 47.37 30.14 22.36 >124 .34 >124.34 3.84 7.82 16.92 118.75 47.37 30.14 22.36 80.23 32.67 7.81
609
96.6 19.4 125.1 43.9 87.1 29.90 84
3.84 12.59 3.84 77.81 .81 3.84 7.82 105.00
Fa (l,N-j) a
K(j)
1168
1381
129 262 300
Models Fitted to Residual Sequences Derived from Runoff-Rainfall Data
STATIONS
MODEL 0. 309w 2(t-5) + n(t+1) (0.003)
Bu11er-Reefton
Wj
wl (t+1)
0.337w2(t+1) 0. 337w 2(t+1) (0.003)
Whirinaki-Minginui
w Wj1 (t+1)
0.103w2 (t) 103w 2(t) 0.360w2(t+1) + 0. (0.001) (0.001)
0.308n(t) + n(t+1) (0.004)
Rangitaiki-Waimihia
Wj
w1 (k+1)
0. 051w 2(t+1) + O. 028w2 (t) (0.0003) (0.0003)
0.OO41w O. 0041 w22 (t-11) (t -11 ) (0.0002)
n(t+1) + nCt+1) W j (.) : w1(e):
Residuals Residua1s from monthly runoff model
e ):: w2 W2C(.)
Residua1s from monthly rainfall model
Numbers in parentheses are the standard errors of parameters.
0.296n(t) 0 . 296n(t) (0.003)
226
A. R. Rao and R. L. Kashyap
which were very small have been ignored. This model may be further refined as desired. Some of the theoretical bases of this type of transfer function model building have been discussed by Haugh and Box (1977). Another important use of models such as those given in Table 4 is in simulation. simulation . As W2(e) w2(-) and nee) sequences in the models given in Table 4 are white noise sequences, they may w1 (-) values be relatively easily simulated. w1(e) thus obtained may be used to generate the sequence of effect variables yet). Thus instead of constructing the usual transfer function models we may use the residual transfer function models to generate sequences of effect variables. Some preliminary work in this direction involving groundwater and precipitation data (Bathala, Rao and Spooner, 1977) has indicated the superiority of this approach over the univariate models. VI.
CONCLUSIONS
On the basis of results presented in this paper we may conclude that the relationships between rainfall and runoff processes may not be as strong as they appear by an analysis of their cross correlation characteristics. The cross correlation characteristics of residuals from validated models fitted to the rainfall and runoff processes and results of tests based on them appear to be better indicators of the causal connections between the processes. These residual sequences may be used to get good candidate transfer function models off the rainfall and runoff processes. ACKNOWLEDGEMENTS Part of A. R. Rao's work reported herein was conducted while he was a Senior Research Fellow of the National Research Advisory Council of New Zealand. I am grateful to the Council for this support.
REFERENCES Bartlett, M. S. (1966). An Introduction to Stochastic Processes, 2nd ed. Cambridge University Press, Cambridge, England. pp. 283-289. Bathala, C. T., A. R. Rao and J. A. Spooner (1977). Applications of Linear Systems Analysis to Groundwater Evaluation Studies. Tech. Rept. No. 91, Water Resources Research Center, Purdue University, W. Lafayette, IN 47907, pp. 69-110. Hannan, E. J. (1970). Multiple Time Series, John Wiley and Sons, Inc., New York, NY. Haugh, L. D. (1976). Checking the Independence of Two Covariance Stationary Time Series: A Uni-Variate Residual Cross-Correlation Approach. J. Am. Stat. Assoc., Vol. 71, 378-385. Haug~ Haug~ D. and G. E. P. Box (1977). Identification of Dynamic Regression (Distributed Lag) Models Connecting Two Time Series. J. Am. Stat. Assoc., Vol. 72, pp. 121-130. Kashyap, R. L. and A. R. Rao (1976). Dynamic Stochastic Models from Empirical Data. Academic Press, New York, pp. 222-226. (1969) . Time Series Analysis Kisiel, C. C. (1969). (Ed.), of Hydrologic Data. In V. T. Chow CEd.), Advances in Hydro Science, Vol. 5, AcaPress . New York, pp. 1-119. demic Press. Rodriguez-Iturbe, I. and V. Yevjevich (1968). The Investigation of Relationship between Hydrologic Time Series and Sunspot Numbers. Colorado State University Hydrology Papers No. 26. Ft. Collins, Colorado.
Causality in Hydrologic Systems
""
227
i
I
I
I
BUtLER AT TEKUHA
I \
~,..,,1.1
(m~JAJCJ;--A_AT_LAN _ _G-.I--I---. .
ANGJTAJKI AT ItNRUPARA
O'}
~
FOIfEST ""'11
WHRINAJ(J
AT GALA TEA
f&10}
FOREST
(BU67f/'
NOMHAN(JO AT MMJ
(31131}
MM~
(Ef&
Fig. 1.
Location of Watersheds
10,.....-------------. 1Or----------------------------,. .
101r---------------------------~ MINGINUI fOREST - WH:Q INAKI [ WJ
WA. IHIHIA FORESTFQRE ST- RANGITAIKI RANGITA IKI[W WAIHIHIA [W
REEF 0 ... RE~FTO
- BULLER !W] BUllE R [loll
~ 06
z ~ ~ ~
~2
0
~
~ -06 ~-06
0 ~
10,00 -80
1~--------------------~
o
~ IM IHI A
--60
10
~O
60
00
10' r---------------------------~ MINGINUI
FORES T- RANGllAl(ilH]
-Ob l AG [ WEEKS)
[ WEEKS]
10)------------------------------
REEF TO'" - BUlLHII'1~
FCREST- WHRINA,l(1
IMl 06 z
g ;:;
01
0 ~
- - --
z
Z 0
et ~
-~
-01
-
--, - - ~
-
-
- - -
--.
~ ~
0
-0 -06 VI
~-O
LAG
MONTHS
u -1 {}-,O:---;';t.O~-3""O--""'20=--""1O,--,O!c---,\1O~2bO,.-,3!"O-I..J,O,......,J"O ~O
Fig. 2.
~
0 '-'
:.;;
! ''''
~
.... S'
.....
-' O_=---..l:-........-,!~--:!_b",....---,\· ~.--+----,L,---="zt~-="3c,....-~ SC
Examples of Cross Correlation Coefficients of Weekly eekly (Top Row) and ~lonthly onthly (Bottom Row) Data
228
A. R. Rao and R. L. Kashyap
10
10r-------.-----..--------, 10,---------------- - - , BULLERIW BULLERI
·Or-------.-----..--------, 1o .-------------------~ RANGITAlKI RAN GI TAlKII[W) WI
r------,----,.-----r----.--_
/
/ /
E
L.::: --- --------.: : -
o -- - ---- - --- -- -
C QI
U
~
C
QI
t"
- - - -- ----.-, - - --0 -\--"7',------
QI
- --- - - --
- --
8'
o
u
'~weeksl
10
"
'() 10
10
/
/
E /
.~
~04
/
'"o
/
~<)~/
"U
~
/
/
~~6
o
Cl
i~
// /
QI
'"
c. 0.
/
/ /
/
/
~~V}/
0·1 02 ~3 03 04 04 Frequency Frequen cy [(PMI I(PMI
Fig. 3.
.
// . / //
/
//
/
...... ~2--003 .....3 --0....... 0 44-....J0.5 05 02 [(py II '~requency requenc y I(PY
01 0·1 0·2 OJ 04 0·4 02 03 Frequency. (CPWJ Frequ ency I(PW}
04 0·4
0-5 OS
'" QI
- - -
--
- /\-
- - - - - -
cC
--
--- - ---- - -----
/
- - - -,:: - - -. 7\ ---
QI ~
~'?-
QI ~
o
- -":' - - - - - ~
/ /
/
/
o
o
u
u LJ
01 0·2 02 0OJ3 Frequency Freq uenc y [( I( PMl PM }
0·4 04
OS 05
-10 . . - - - - - - :fA-----.b:-------.!30 -, 0 L L-----~'hr-O----~20Or-------.1 1O Log [months) Imon th s}
-, 0);0 ------;'-:;-0----~20<;;--------.I -100~--~1'A0---::2~0-----"J30 Log Imonths} (months)
Examples of Correlograms and Cumulative Periodograms of Residuals _..
u: u.: l.L u. w W
o uU
u.i W
Q: 0:: Q: 0:: 0 D-
U
tn Vl tn Vl
oQ:
0::
II
...
OS 05
100 r-------.-----B-U,....LL-E-R-[-H-I---, 1 ,---~--BU~ LL-E-R...,.IM-:I--,
C c
/
/
w
/
// /
/
/
/ / /
/
/
00·66
3" 0·2
/
/
,
,
/
E
/
02
,
// /
/
E
/
d
E
.3
1·0.-------,----..--------, RANGITAIKI IM} rH)
/
~/
~o
/
/
//
/
&
/
/
/
/ //
~\. / //, ~~\9
'"
0
//
o000~---"'0'-1 0 1-
30
BULLER [Ml//
OB 0·8
/
//
~o
/
/
Log [weeks)} Lo g Iweeks
1o.---------,r---R-r- G-I-TA"'T"'-KI-[-MY"T'""-:"/----., AN
0·8
20
/
/
/
/ /
/
/
/
/
LJ W
30
20 20
Log [weeks)
/
02 02
::::l
~
/
//
/
Cj 0, /
/
E E
o
o
u
// / /
/
/ h\o/ ~\o/ / // /
06 04
i:
~
o
LJ
h
//
BULLER (Wl // IW I /// /
/
/
0·8 08 III
10.---------,-----..--..----~10
RA GITAlKI [WI / / RANGI TA IKllwl
D.-
...
.ll
¥
-
... ...
"m
(WEEKS ) LAG (WEEKS)
<~
+
U
-1
u: u. u. l.L
w W
0o
U
l:l
W u.i Q:
0:: Q: 0:: 0D U (J) Vl (J) Vl
o0Q:
+
LAG (WEEKS)
0::: U
I
u: u.: l.L u. W w o u U
W W
--60-------- - -4-------.0.-"A LI.
B_,
(f) Vl (J) Vl
o
Q: 0:: U
6___6_~".
.t._m__-- ~ ~P~~~_ILL
---------------------llllll
""
_ _ A -I>..I:.-~~-- ~~~.~~ ~ __
LAG (MONTHS)
"-A-
--:-;'5
+
Fig. 4. Cross Correlograms of Residuals from Weekly Weekl y and Monthly Monthl y Models: Rangitaiki-Waimihia (weekly, top), Whirinaki-Minginui (weekly, middle), middle ) , and BUller-Reefton Buller-Reefton (monthly, bottom). Positive and negative lags respectively respectiv ely indicate runoff and rainfall as leading variables. vari ables.