Testing and modelling autoregressive conditional heteroskedasticity of streamflow processes

Nonlinear Processes in Geophysics (2005) 12: 55–66 SRef-ID: 1607-7946/npg/2005-12-55 European Geosciences Union © 2005 Author(s). This work is licensed under a Creative Commons License. Nonlinear Processes in Geophysics Testing and modelling autoregressive conditional heteroskedasticity of streamflow processes W. Wang1, 2 , P. H. A. J. M. Van Gelder2 , J. K. Vrijling2 , and J. Ma3 1 Faculty of Water Resources and Environment, Hohai University, Nanjing, 210098, China of Civil Engineering & Geosciences, Section of Hydraulic Engineering, Delft University of Technology, P.O.Box 5048, 2600 GA Delft, The Netherlands 3 Yellow River Conservancy Commission, Hydrology Bureau, Zhengzhou, 450004, China 2 Faculty Received: 24 May 2004 – Revised: 15 December 2004 – Accepted: 5 January 2005 – Published: 21 January 2005 Part of Special Issue “Nonlinear deterministic dynamics in hydrologic systems: present activities and future challenges” Abstract. Conventional streamflow models operate under the assumption of constant variance or season-dependent variances (e.g. ARMA (AutoRegressive Moving Average) models for deseasonalized streamflow series and PARMA (Periodic AutoRegressive Moving Average) models for seasonal streamflow series). However, with McLeod-Li test and Engle’s Lagrange Multiplier test, clear evidences are found for the existence of autoregressive conditional heteroskedasticity (i.e. the ARCH (AutoRegressive Conditional Heteroskedasticity) effect), a nonlinear phenomenon of the variance behaviour, in the residual series from linear models fitted to daily and monthly streamflow processes of the upper Yellow River, China. It is shown that the major cause of the ARCH effect is the seasonal variation in variance of the residual series. However, while the seasonal variation in variance can fully explain the ARCH effect for monthly streamflow, it is only a partial explanation for daily flow. It is also shown that while the periodic autoregressive moving average model is adequate in modelling monthly flows, no model is adequate in modelling daily streamflow processes because none of the conventional time series models takes the seasonal variation in variance, as well as the ARCH effect in the residuals, into account. Therefore, an ARMA-GARCH (Generalized AutoRegressive Conditional Heteroskedasticity) error model is proposed to capture the ARCH effect present in daily streamflow series, as well as to preserve seasonal variation in variance in the residuals. The ARMA-GARCH error model combines an ARMA model for modelling the mean behaviour and a GARCH model for modelling the variance behaviour of the residuals from the ARMA model. Since the GARCH model is not followed widely in statistical hydrology, the work can be a useful adCorrespondence to: W. Wang (w.wang@126.com) dition in terms of statistical modelling of daily streamflow processes for the hydrological community. 1 Introduction to autoregressive conditional heteroskedasticity When modelling hydrologic time series, we usually focus on modelling and predicting the mean behaviour, or the first order moments, and are rarely concerned with the conditional variance, or their second order moments, although unconditional season-dependent variances are usually considered. The increased importance played by risk and uncertainty considerations in water resources management and flood control practice, as well as in modern hydrology theory, however, has necessitated the development of new time series techniques that allow for the modelling of time varying variances. ARCH-type models, which originate from econometrics, give us an appropriate framework for studying this problem. Volatility (i.e. time-varying variance) clustering, in which large changes tend to follow large changes, and small changes tend to follow small changes, has been well recognized in financial time series. This phenomenon is called conditional heteroskedasticity, and can be modeled by ARCH-type models, including the ARCH model proposed by Engle (1982) and the later extension GARCH (generalized ARCH) model proposed by Bollerslev (1986), etc. Accordingly, when a time series exhibits autoregressive conditionally heteroskedasticity, we say it has the ARCH effect or GARCH effect. ARCH-type models have been widely used to model the ARCH effect for economic and financial time series. The ARCH-type model is a nonlinear model that includes past variances in the explanation of future variances. ARCH- 0 Day 10000 15000 W. Wang et al.: Testing and modelling autoregressive conditional heteroskedasticity Figure 1 Daily streamflow (m3/s) of the upper Yellow River at Tangnaihai 1000 2000 3000 4000 3 Discharge (m /S) 5000 56 1600 1400 1200 1000 800 600 400 200 0 1-Jan daily mean standard deviation 2-Mar 1-May 30-Jun Date 29-Aug 28-Oct 27-Dec Fig. 2. Variation in daily mean and standard deviation of the stream- 0 5000 Day 10000 15000 1.0 0 Figure 2 Variation in daily mean and standard deviation of the streamflow at Tangn flow at Tangnaihai. 1.0 Discharge (cms) 5000 0.8 0.8 tributing watershed, including a permanently snow-covered area of 192 km2 . The length of the main channel of this waFigure 1 Daily streamflow (m /s) of the upper Yellow River at Tangnaihai tershed is over 1500 km. Most of the area is 3000∼6000 naihai. meters above sea level. Snowmelt water composes about 5% 1600 of total runoff. Most rain falls in summer. Because the water1400 type models can generate accurate forecasts of future volatilshed is partly permanently snow-covered and sparsely popdaily mean ity,1200 especially over short horizons, therefore providing a betulated, without any large-scale hydraulic works, it is fairly standard deviation ter1000 estimate of the forecast uncertainty which is valuable for pristine. The average annual runoff volume (during 1956– 800resource management and flood control. And they take water 2000) at Tangnaihai gauging station is 20.4 billion cubic meinto600 account excess kurtosis (i.e. fat tail behaviour), which ters, about 35% of the whole Yellow River Basin, and it is the 40 60 80 100 0 20 40 60 400 is common in hydrologic processes. Therefore, 0ARCH-20 major runoff producing area of the Yellow River basin. Daily Lag Lag 200 type models could be very useful for hydrologic time seaverage streamflow at Tangnaihai has been recorded since 1 0 ries modelling. Some authors propose new models to reproJanuary 1956. Monthly series is obtained from daily data by 1-Jan 2-Mar 1-May 30-Jun 29-Aug 28-Oct 27-Dec Figure 3 ACF and PACF of deseasonalized daily flow series duce the asymmetric periodic Date behaviour with large fluctuataking the average of daily discharges in every month. In this tions around large streamflow and small fluctuations around study, data from 1 January 1956 to 31 December 2000 are small streamflow (e.g. Livina et al., 2003), which basically used. The daily streamflow series from 1956 to 2000 is plot2 Variationcan in daily mean with and standard deviationtime of the streamflow at be handled those conventional series modtedTangnaihai in Fig. 1, and variations in the daily mean discharge and els that have taken season-dependent variance into account, daily standard deviation of the streamflow at Tangnaihai are such as PARMA models and deseasonalized ARMA models. shown in Fig. 2. However, little attention has been paid so far by the hydrologic community to test and model the possible presence of 3 Tests for the ARCH effect of streamflow process the ARCH effect with which large fluctuations tend to follow large fluctuations, and small fluctuations tend to follow small The detection of the ARCH effect in a streamflow series is fluctuations in streamflow series. actually a test of serial independence applied to the serially In this paper, we will take the daily and monthly streamuncorrelated fitting error of some model, usually a linear auflow of the upper Yellow River at Tangnaihai in China as toregressive (AR) model. We assume that linear serial depencase study hydrologic time series to test for the existence dence inside the original series is removed with a well-fitted, of the ARCH effect, and propose an ARMA-GARCH error pre-whitening model; any remaining serial dependence must model for daily flow series. The paper is organized as folbe due to some nonlinear generating mechanism which is lows. First, of testing conditional 20 40 60 the method 80 100 0 20 heteroskedas40 60 80 100 not captured by the model. Here, the nonlinear mechanism Lag Lag ticity of streamflow process is described. Then, the causes of we are concerned with is the conditional heteroskedasticthe ARCH effect and the inadequacy of commonly used seaity. We will show that the nonlinear mechanism remaining sonal time series modelsoffordeseasonalized modelling streamflow are disFigure 3 ACF and PACF daily flow series in the pre-whitened streamflow series, namely the residual cussed. Finally, an ARMA-GARCH error model is proposed series, can be well interpreted as autoregressive conditional for capturing the ARCH effect existing in daily streamflow heteroskedasticity. series. -0.2 -0.0 -0.2 -0.0 Partial ACF 0.2 0.4 0.6 0.8 1.0 0.0 3 0.2 Discharge (m /S) ACF 0.4 0.6 Partial ACF 0.2 0.4 0.6 Fig. 1. Daily streamflow3 (m3 /s) of the upper Yellow River at Tang- 3.1 Linear ARMA models fitted to daily and monthly flows 2 Case study area and data set The case study area is the headwaters of the Yellow River, located in the northeastern Tibet Plateau. In this area, the discharge gauging station Tangnaihai has a 133 650 km2 con- Three types of seasonal time series models are commonly used to model hydrologic processes which usually have strong seasonality 18 (Hipel and McLeod, 1994): 1) seasonal autoregressive integrated moving average (SARIMA) mod- 80 0 1-Jan 2-Mar 1-May 30-Jun 29-Aug Date 28-Oct 27-Dec W. Wang et al.: Testing and modelling autoregressive conditional heteroskedasticity 57 0.0 -0.2 -0.0 0.2 ACF 0.4 0.6 Partial ACF 0.2 0.4 0.6 0.8 0.8 1.0 1.0 Figure 2 Variation in daily mean and standard deviation of the streamflow at Tangnaihai 0 20 40 Lag 60 80 0 20 40 Lag 60 80 100 of deseasonalized daily flow series 0.8 0.6 Partial ACF 0.2 0.4 18 0.0 0.0 0.2 ACF 0.4 0.6 1.0 Fig. 3. ACF and PACF of deseasonalized flowand series. Figure 3daily ACF PACF 100 0 10 20 30 Lag 40 50 Fig. 4. ACF and PACF of deseasonalized monthly flow series. of Figure 4 ACF and PACF 60 0 10 20 30 Lag 40 50 60 deseasonalized monthly flow series Residuals -1.5 -1.0 -0.5 0.0 0.5 0.2 0.1 ACF 0.0 0.00 ACF 0.05 0.10 -1.0 -0.5 Residuals 0.0 0.5 1.0 1.0 1.5 els; 2) deseasonalized ARMA models; and 3) periodic Firstly, we inspect the ACF of the residuals. It is wellARMA models. The deseasonalized modelling approach is known that for random and independent series of length n, adopted in this study. The procedure of fitting deseasonalized the lag k autocorrelation coefficient is normally distributed ARMA models to daily and monthly streamflow at Tangwith a mean of zero and a variance of (a) (b) and the 95% √ 1/n, naihai includes two steps. First, logarithmize both flow seconfidence limits are given by ±1.96/ n. The ACF plots ries, and deseasonalize them by subtracting the seasonal (e.g. in Fig. 6 show that there is no significant autocorrelation left daily or monthly) mean values and dividing by the seasonal in the residuals from both ARMA-type models for daily and standard deviations of the logarithmized series. To alleviate monthly flow. the stochastic fluctuations of the daily means and standard Then, more formally, we apply the Ljung-Box test (Ljung deviations, we smooth them with first 8 Fourier harmonics and Box, 1978) to the residual series, which tests whether before using them for standardization. Then, according to the the first L autocorrelations r̂k2 (ε2 ) (k = 1, ..., L) from a proACF (AutoCorrelation Function) and PACF (Periodic Autocess are collectively small in magnitude. Suppose we have 0 structures 200 800well as 1000 40 60 Correlation Function) of400 theDay two600 series, as the 0first L20autocorrelations r̂k80 (ε) (k100 = 1, 120 ..., L) from any Month the model selection criterion AIC, two linear ARMA-type ARMA(p, d, q) process. For a fixed sufficiently large L, models (one ARMA(20,1) and one AR(4)) are fitted to the the usual Ljung-Box Q-statistic is given by Figure Segments of daily the residual series logarithmized and5deseasonalized and monthly flowfrom se- (a) ARMA(20,1) for daily flow and (b) AR(4) L X r̂k2 (ε) ries, respectively, following theTangnaihai. model construction procefor monthly flow at Q = N(N + 2) , (1) dures suggested by Box and Jenkins (1976). Figures 3 and N −k k=1 4 show the ACF and PACF of the deseasonalized daily and monthly series. Figure 5 shows parts of the two residual sewhere N = sample size, L= the number of autocorrelations (a) (b) sample autories obtained from the two models. included in the statistic, and r̂k2 is the squared Before applying ARCH tests to the residual series, to encorrelation of residual series {εt } at lag k. Under the null sure that the null hypothesis of no ARCH effect is not rehypothesis of model adequacy, the test statistic is asympjected due to the failure of the pre-whitening linear models, totically χ 2 (L−p−q) distributed. Thus, we would reject we must check the goodness-of-fit of the linear models. the null hypothesis at level α if the value of Q exceeds the Figure 4 ACF and PACF of deseasonalized monthly flow series -1.0 -0.5 (a) 200 400 -1.0 0 0 200 400 Day Day 600 600 800 1000 800 1000 1.5 (a) Residuals Residuals -1.5 -0.5 -1.0 0.0 -0.5 0.5 0.0 1.00.5 1.51.0 -1.5 -1.0 1.0 W. Wang et al.: Testing and modelling autoregressive conditional heteroskedasticity Residuals Residuals -0.5 0.0 0.5 0.0 0.5 1.0 58 (b) (b) 0 20 40 0 20 40 60 Month 60 Month 80 80 100 120 100 120 Figure 5 Segments of the residual series from (a) ARMA(20,1) for daily flow and (b) AR(4) Fig. 5. Segments of5 the residual fromresidual (a) ARMA(20,1) daily(a) flowARMA(20,1) and (b) AR(4) forfor monthly at Tangnaihai. Figure Segments of the series for from dailyflow flow and (b) AR(4) for monthly flow atseries Tangnaihai. 0.2 0.10 for monthly flow at Tangnaihai. 0.05 0.00 ACF 0.00 ACF 0.05 (a) 0 100 0 100 200 Lag 200 Lag (b) ACF ACF 0.0 0.0 0.1 0.1 0.2 0.10 (a) (b) 300 0 2 4 300 0 2 4 6 Lag 6 Lag 8 8 10 10 12 12 Fig. 6. ACFs of residuals from (a) the ARMA(20,1) model for daily flow and (b) the AR(4) model for monthly flow at Tangnaihai. Figure 6 ACFs of residuals from (a) ARMA(20,1) model for daily flow and (b) AR(4) model Figure 6 ACFs from (a) ARMA(20,1) model for daily flow and (b) AR(4) model for monthly flowofatresiduals Tangnaihai (1−α)-quantile of the χ 2 (L−p−q) distribution. The Ljungfor monthly flow at Tangnaihai Box test results for ARMA(20,1) and AR(4) are shown in Fig. 7. The p-values’ exceedance of 0.05 indicates the acceptance of the null hypothesis of model adequacy at significance level 0.05. However, while the residuals seem statistically uncorrelated according to ACF and PACF shown in Fig. 6, they are not identically distributed from visual inspection of the Fig. 5, that is, the residuals are not independent and identically distributed (i.i.d.) through time. There is a tendency, especially for daily flow, that large (small) absolute values of the residual process are followed by other large (small) values of unpredictable sign, which is a common behaviour of GARCH processes. Granger and Andersen (1978) found that some of the series modelled by Box and Jenkins (1976) exhibit autocorrelated squared residuals even though the residuals themselves do no seem to be correlated over time, and therefore suggested that the ACF of the squared time series could be useful in identifying nonlinear time series. Bollerslev (1986) stated that the ACF and PACF of squared process are useful in identifying and checking GARCH behaviour. Figure 8 shows the ACFs of the squared residual series from the ARMA(20,1) model for daily flow and the AR(4) model for monthly flow at Tangnaihai. It is shown that although the residuals are almost uncorrelated, as shown in Fig. 6, the squared residual series are autocorrelated, and the ACF structures of both squared residual series exhibit strong seasonality. This indicates that the variance of residual series is conditional on its past history, namely, the residual series may exhibit an ARCH effect. There are some formal methods to test for the ARCH 19 effect of a process, such as the McLeod-Li test (McLeod 19 and Li, 1983), the Engle’s Lagrange Multiplier test (Engle, 1982), the BDS test (Brock et. al., 1996), etc. McLeodLi test and Engle’s Lagrange Multiplier test are used here to check the existence of an ARCH effect in the streamflow series. 3.2 McLeod-Li test for the ARCH effect McLeod and Li (1983) proposed a formal test for ARCH effect based on the Ljung-Box test. It looks at the autocorrelation function of the squares of the pre-whitened data, and tests whether the first L autocorrelations for the squared residuals are collectively small in magnitude. Similar to Eq. (1), for fixed sufficiently large L, the LjungBox Q-statistic of Mcleod-Li test is given by Q = N(N + 2) L X r̂ 2 (ε 2 ) k k=1 N −k , (2) 0. 3 (b) 0. 25 0. 2 p-value p-value (a) 1 0. 8 0. 6 0. 15Testing and modelling autoregressive conditional heteroskedasticity W. Wang et al.: 0. 4 0. 1 0. 05 0. 2 0. 3 (b)0 26 0. 15 Lag 31 36 p-value 0. 25 21 0. 2 p-value (a)0 0. 1 59 1 0.4 8 6 8 10 12 Lag 0. 6 14 16 18 20 0. 4 0. 2 Figure 7 0.Ljung-Box lack-of-fit tests for (a) ARMA(20,1) model for daily flow and (b) AR(4) 05 0 0 model for monthly flow. 21 26 Lag 31 4 36 6 8 10 12 Lag 14 16 18 20 Fig. 7. Ljung-Box lack-of-fit tests for (a) the ARMA(20,1) model for daily flow and (b) the AR(4) model for monthly flow. Figure 7 Ljung-Box lack-of-fit tests for (a) ARMA(20,1) model for daily flow and (b) AR(4) 0.2 0.2 model for monthly flow. (b) 0 0.0 0.2 0.1 (b) ACF 0.0 ACF 0.1 0.0 (a) ACF 0.0 ACF 0.2 0.1 0.1 (a) 100 200 Lag 300 0 2 4 6 Lag 8 10 12 Fig. 8. ACFs of the8squared (a) the ARMA(20,1) model(a) for daily flow and (b) the AR(4)for model for monthly flow at(b) Tangnaihai. Figure ACFsresiduals of thefrom squared residuals from ARMA(20,1) model daily flow and 0 100 200 300 AR(4) model for monthly flowLagat Tangnaihai 0 r̂k2 2 4 6 Lag 8 10 12 p-value p-value p-value p-value where N is the sample size, and is the squared sample lags less than 4 are removed by the AR(4) model, when we autocorrelation of squared residual series at lag k. Under autocorrelations at longer lags into consideration, sigFigure 8 ACFs of the squared residuals from take (a) ARMA(20,1) model for daily flow and (b) the null hypothesis of a linear generating mechanism for the nificant autocorrelations remain and the null hypothesis of data, namely, no ARCH effectfor in monthly the data, the testatstatistic is no ARCH effect is rejected. Because it is required for the AR(4) model flow Tangnaihai asymptotically χ 2 (L) distributed. Figure 9 shows the results McLeod-Li test to use sufficiently large L, namely, a suf0. 8 (a) 0. 06test for daily and monthly flow. It illus- (b)ficient of the McLeod-Li number of autocorrelations to calculate the Ljung0. 7 0. 05 trates that the null hypothesis of no ARCH effect is rejected Box statistic (typically around 20), we still consider that the 0. 6 0. 04 0. 5 for both daily and monthly flow series. monthly flow has the ARCH effect. 0. 4 the whole, evidences are clear with the McLeod-Li test 0. 03 On 0. 06 0. 3 0. 8 (a) (b) 0. 02 3.3 Engle’s Lagrange Multiplier test for the ARCH effect 0. 7 and LM test about the existence of conditional het0. 2Engle’s 0. 05 0. 6 0. 01 eroskedasticity in the residual series from linear models fitted 0. 1 0. 04 0. 5 Since the ARCH 0model has the form of an autoregres0.the 0 logarithmized and deseasonalized daily and monthly to 0. 4 0. 03 sive model, Engle (1982) proposed Multiplier 5 10 of the 15 upper 20 Yellow 25 River 30 at Tangnai0 5 10 the15Lagrange 20 25 30 0. 0 3 streamflow processes 0. 02 Lag Lag 0. 2 (LM) test, in order to test for the existence of ARCH behai. 0. 01 0. 1 haviour based on the regression. The test statistic is given 0. 0 0 by TR2 , where R is the sample multiple correlation coef0 5 model 10 for 15 daily 20 flow 25 and30 Figure 9 McLeod-Li residuals 0 5test for 10 the 15 20 25from30(a) ARMA(20,1) 2 ficient computed from the regression of εLag Lagof ARCH effects and inade4 Discussion of the causes t on a constant 2 ,. . . ,ε 2 , and T is the sample size. Under the null and εt−1 quacy of commonly used seasonal time series models t−q (b) AR(4) model for monthly flow hypothesis that there is no ARCH effect, the test statistic Figure 9 McLeod-Li test for the residuals (a) ARMA(20,1) model for daily flow and is asymptotically distributed as chi-square distribution withfrom4.1 Causes of ARCH effects in the residuals from ARMAq degrees of freedom. As Bollerslev (1986) suggested, it type models for daily and monthly flow (b)power AR(4)against modelGARCH for monthly flow should also have alternatives. Figure 10 shows Engle’s LM test results for the residuFrom the above analyses, it is clear that although20the residals from the ARMA(20,1) model for daily flow and from the uals are serially uncorrelated, they are not independent AR(4) model for monthly flow. The results also firmly inthrough time. At the mean time, we notice that seasonal20residual dicate the existence of an ARCH effect in both the residual ity dominates autocorrelation structures of squared series. series for both daily and monthly flow processes (as shown in Fig. 8). This suggests that there are seasonal variations in One point that should be noticed is that although Figs. 8b, the variance of the residual series, and we should standardize 9b and 10b show that for monthly flow, autocorrelations at AR(4) model for monthly flow at Tangnaihai 60 W. Wang et al.: Testing and modelling autoregressive conditional heteroskedasticity (a) (b) 0. 06 0. 04 p-value p-value 0. 05 0. 03 0. 02 0. 01 0 0 5 10 15 Lag 20 25 0. 8 0. 7 0. 6 0. 5 0. 4 0. 3 0. 2 0. 1 0. 0 0 30 5 10 15 Lag 20 25 30 Fig. 9. McLeod-Li test for the residuals from (a) the ARMA(20,1) model for daily flow and (b) the AR(4) model for monthly flow. Figure 9 McLeod-Li test for the residuals from (a) ARMA(20,1) model for daily flow and 0. 06 model for monthly flow (b) (a)AR(4) 8 (b) 0. 0. 7 0. 04 p-value p-value 0. 05 0. 03 0. 02 0. 01 0 0 5 10 15 Lag 20 25 0. 6 0. 5 0. 4 0. 3 0. 2 0. 1 0. 0 20 0 30 5 10 15 Lag 20 25 30 Fig. 10. Engle’s LM test for residuals from (a) the ARMA(20,1) model for daily flow and (b) the AR(4) model for monthly flow. Figure 10 Engle’s LM test for residuals from (a) ARMA(20,1) model for daily flow and SD 0.0 0.0 ACF 0.1 ACF 0.1 0.2 0.2 0.3 SD AR(4) for monthly flow standard From the above analyses, it is clear that the ARCH effect is the residual (b) series from model linear models with seasonal fully caused by seasonal variances for monthly flow, but only deviations of the residuals first, then look at the standardpartly ized series to check seasonal variances can explain 0.4 (a) whether (b) 1for daily flow. Other causes, besides the seasonal variation in variance, of the ARCH effect in daily flow may inARCH effects. 0. 8 0.3 clude the perturbations of the temperature fluctuations which Seasonal standard deviations of the residual series from is an0. 6influential factor for snowmelt, as well as evapotranthe ARMA(20,1) model 0.2 for daily flow and the AR(4) model spiration, and the precipitation variation which is the domifor monthly flow are calculated and shown in Figs. 11a and 0. 4 0.1 nant0.factor for streamflow processes. As reported by Miller 11b. They are used to standardize the residual series from 2 (1979), when modelling a daily average streamflow series, the ARMA(20,1) model and the AR(4) model. Figure 12 0 0 the residuals from a fitted AR(4) model signaled white-noise shows the ACFs of the squared series 0 60 standardized 120 180 residual 240 300 360 0 3 6 9 12 errors, but the squared residuals signaled bilinearity. When of daily and monthly flow. It is illustrated Month Day that, after seaprecipitation covariates were included in the model, Miller sonal standardized autocorrelation, as well as the seasonality found that neither the residuals nor the squared residuals sigin the squared standardized residual series for monthly flow Seasonal standard deviations (SD) of the residuals form (a)While ARMA(20,1) model naled any problems. we agree that the for autocorrelais basically Figure removed11(Fig. 12b), the significant autocorrelation existing in the squared residuals is basically caused by tion still exists in the squared standardized residual for flow daily flow and (b) AR(4) model for series monthly a precipitation process, we want to show that the autocorredaily flow (Fig. 12a), despite the fact that the autocorrela(Note: the smoothed linewith in Figure 11(a) the squared first 8 residuals Fourier harmonics the by an lationby in the can be well of described tions are significantly reduced compared Fig. 8a and the is given ARCH model, which is very close to the bilinear model (Enseasonality in the ACFSD structure is removed. This means that seasonal series.) gle, 1982). the seasonality, as well as the autocorrelation in the squared residuals from the AR model of monthly flow series is basi4.2 Inadequacy of commonly used seasonal time series cally caused by seasonal variances. But seasonal variances models for modelling streamflow processes only explain partly the autocorrelation in the squared residu(a) (b) als of daily flow series. As mentioned in Sect. 3.1, SARIMA models, deseasonalThe residual series of daily flow and monthly flow stanized ARMA models and periodic models are commonly used dardized by seasonal standard deviation are also tested for to model hydrologic processes (Hipel and McLeod, 1994). ARCH effects with the McLeod-Li test and Engle’s LM Given a time series (x ), the general form of SARIMA model, t test. Figure 13 shows that the seasonally standardized denoted by SARIMA(p,d,q)×(P,D,Q) , is S residual series of daily flow still cannot pass the LM test (Fig. 13a), whereas the seasonally standardized residual seφ(B)8(B s )∇ d ∇sD xt = θ(B)2(B s )εt , (3) ries of monthly flow pass the LM test with high p-values (Fig. 13b). The McLeod-Li test gives similar results. where φ(B) and θ(B) of orders p and q represent the ordis nary autoregressive and 0 5 moving average 10 components; 8(B ) 0 100 200 300 Lag Lag Figure 12 ACFs of squared seasonally standardized residuals from (a) ARMA(20,1) model 0. 8 0.3 0. 6 SD SD Figure 10 Engle’s LM test for residuals from (a) ARMA(20,1) model for daily flow and 0.2 0. 4 (b)al.:AR(4) model for monthly flow conditional heteroskedasticity W. Wang et Testing and modelling autoregressive 0.1 0. 2 60 120 180 Day 240 300 (b) 0 1 0 0. 8 360 SD SD (a) 00.4 0 0.3 61 0.2 3 6 Month 9 12 0. 6 0. 4 Figure 11 Seasonal standard deviations (SD) of the residuals form (a) ARMA(20,1) model for 0.1 0. 2 daily flow and (b) AR(4) model for monthly flow 0 0 0 60 120 180 240 300 360 0 6 9 12of the (Note: the smoothed line in Figure 11(a) is given by the first3 8 Fourier harmonics Month Day seasonal SD series.) Fig. 11. Seasonal standard deviations (SD) of the residuals from (a) the ARMA(20,1) model for daily flow and (b) the AR(4) model for Figure Seasonalline standard deviations (SD) of the residuals form (a) SD ARMA(20,1) model for monthly flow (note: 11 the smoothed in (a) is given by the first 8 Fourier harmonics of the seasonal series). 0.2 0.3 daily flow and (b) AR(4) model for monthly flow (b) of the (Note: the smoothed line in Figure (a) 11(a) is given by the first 8 Fourier harmonics 200 Lag 300 0 5 Lag 10 0.0 100 ACF ACF 0.1 0 (b) 0.1 0.0 0.2 (a) 0.0 0.2 ACF 0.1 0.3 ACF 0.1 0.2 seasonal SD series.) 0.0 Fig. 12. Figure ACFs of 12 squared seasonally standardized residuals standardized from (a) the ARMA(20,1) daily flow and (b) themodel AR(4) model for ACFs of squared seasonally residualsmodel fromfor(a) ARMA(20,1) monthly flow. for daily flow and (b) AR(4) model for monthly flow and 2(B s ) of orders P and Q represent the seasonal autore0 100 200 d 300 d gressive and moving average components; Lag ∇ =(1−B) and D s D ∇S =(1−B ) are the ordinary and seasonal difference components. resent a day, week, month or season), we have the following 0 5 10 21 PAR(p) model (Salas, 1993): Lag p X x = µ + φj,s (xv,s−j µs−j ) + εn,s , model (5) n,s s Figure 12 ACFs of squared seasonally standardized residuals from (a) −ARMA(20,1) j =1 The general form of the ARMA(p, q) model fitted to deforseries dailyisflow and (b) AR(4) model for monthly flow seasonalized where εn,s is an uncorrelated normal variable with mean zero and variance σ 2s . For daily streamflow series, to make the model parsimonious, we can cluster the days in21the year φ(B)xt = θ(B)εt . (4) into several groups and fit separate AR models to separate groups (Wang et al., 2004). Periodic models would perFrom the model equations we know that although the seaform better than the SARIMA model and the deseasonalsonal variation in the variance present in the original time ized ARMA model for capturing the ARCH effect, because series is basically dealt with well by the deseasonalized apit takes season-varying variances into account. However, proach, the seasonal variation in variance in the residual seas analyzed in Sect. 4.1, while considering seasonal variries is not considered by either of the two models, because in ances could be sufficient for describing the ARCH effect in both cases the innovation series ε t is assumed to be i.i.d. N(0, monthly flow series because the ARCH effect in monthly σ 2 ). Therefore, both SARIMA models and deseasonalized flow series is fully caused by seasonal variances, it is still models cannot capture the ARCH effect that we observed in insufficient to fully capture the ARCH effect in daily flow the residual series. series. In contrast, the periodic model, which is basically a group In summary, while the PARMA model is adequate for of ARMA models fitted to separate seasons, allows for seamodelling the variance behaviour for monthly flow, none of sonal variances in not only the original series but also the the commonly used seasonal models is efficient enough to residual series. Taking the special case PAR(p) model (pedescribe the ARCH effect for daily flow, although PARMA riodic autoregressive model of order p) as an example of can partly describe it by considering seasonal variances. It a PARMA model, given a hydrological time series xn,s , in is necessary to apply the GARCH model to achieve the purwhich n defines the year and s defines the season (could reppose. 62 15 Lag 20 0. 05 p-value 10 0. 8 0. 6 W. Wang et al.: Testing and modelling autoregressive conditional heteroskedasticity 0. 4 0. 2 0. 06 (a) 1. 0 (b) 0. 0 0. 04 0 250. 03 30 0. 02 5 10 1. 0 0. 8 p-value p-value (b) 15 0. 620 Lag 0. 4 0. 01 0. 2 0 0. 0 25 30 ’s LM test for seasonally standardized residuals from30 (a) ARMA(20,1) 0 5 10model 15 0 5 10 15 20 25 Lag Lag d (b) from AR(4) model for monthly flow 20 25 30 Fig. 13. Engle’s LM test for seasonally standardized residuals from (a) the ARMA(20,1) model for daily flow and (b) from the AR(4) model Figure 13 Engle’s LM test for seasonally standardized residuals from (a) ARMA(20,1) model for monthly flow. GARCH(p, q) model has the form (Bollerslev, 1986)  εt |ψt−1 ∼ N(0, ht )  q p P P , 2 + αi εt−i βi ht−i  ht = α0 + i=1 Partial ACF 0.05 0.10 0.15 Partial ACF 0.05 0.10 0.15 0.20 0.20 for daily flow and (b) from AR(4) model for monthly flow (6) i=1 Conditional standard deviation 1.0 1.5 2.0 2.5 3.0 Conditional standard deviation 1.0 1.5 2.0 2.5 3.0 -4 -2 Residuals 0 2 4 6 0.00 0.00 where, ε t denotes a real-valued discrete-time stochastic process, and ψ t the available information set, p ≥0, q >0, α0 >0, αi ≥0, β i ≥0. When p=0, the GARCH(p,q) model reduces to the ARCH(q) model. Under the GARCH(p, q) model, the conditional variance of εt , ht , depends on the squared residuals in the previous q time steps, and the conditional variance in the previous p time steps. Since GARCH be treated 0 20 40 models 60 can 80 100 as ARMA models for squared residuLag 0 20 40 60 80 100 als, the order of GARCH can be determined with the method Lag for selecting the order of ARMA models, and traditional model selection criteria, such as Akaike information criterion Figure 14 PACF of the squared seasonally starndardized residual series from ARMA(20,1) Fig. 14. PACF of the squared seasonally starndardized residual se(AIC) and Bayesian information criterion (BIC), can also be F of the squared seasonally starndardized residual series from ARMA(20,1) ries from ARMA(20,1) daily flow. for dailyfor flow used for selecting models. The unknown model parameters α i (i = 0, · · · , q) and β j (j = 1, · · · , p) can be estimated using (conditional) maximum likelihood estimation (a) (b) 1/2 (MLE). Estimates of the conditional standard deviation ht 5 Modelling the daily steamflow with ARMA-GARCH are also obtained as a side product with the MLE method. error model (a) When there is(b) obvious seasonality present in the residuals (as in the case of daily streamflow at Tangnaihai), to preserve the seasonal variances in the residuals, instead of fitting the 5.1 Model building ARCH model to the residual series directly, we fit the ARCH model to the seasonally standardized residual series, which is obtained by dividing the residual series by seasonal stanWeiss (1984) proposed ARMA models with ARCH errors. dard deviations (i.e. daily standard deviations for daily flow). This approach is adopted and extended by many researchers the 0 200general 400ARMA-GARCH 600 800 model 1000with seasonal 0 time 200 600 and 800Kunst, 1000 Therefore, for modelling economic series 400 (e.g. Hauser Day Day standard deviations we propose here has the following form 1998; Karanasos, 2001). In the field of geo-sciences, Tol (1996) fitted a GARCH model for the conditional variance  φ(B)xt = θ (B)εt and the conditional standard deviation, in the conjunction with Figure 15 A segment of (a) seasonally standardized from ARMA(20,1) and (b)  ε =residuals σs zt , zt ∼N (0, ht ) t an AR(2) model for the mean, to model daily mean temper, (7) q p P P  corresponding standard deviation sequence estimated ARCH(21) model 0 ARMA-GARCH 200 400 600 800 1000 0 400 ature. 600 800 2with In thisitspaper, we1000 proposeconditional to use erh = α + α z + β h t i t−i i t−i 0 Day Day i=1 i=1 ror (or, for notation convenience, called ARMA-GARCH) model for modelling daily streamflow processes. where σ s is the seasonal standard deviation of εt , s is the season number depending on which season the time t belongs to. The ARMA-GARCH model may be interpreted as a comment of (a) the seasonally standardized residuals from ARMA(20,1) and (b) For daily series, s ranges from 1 to 366. Other notations are bination of an ARMA model which is used to model mean 22 the same as in Eqs. (4) and (6). behaviour, and an ARCH model which is used to model g conditional standard deviation sequence estimated with ARCH(21) model The model building procedure proceeds in the following the ARCH effect in the residual series from the ARMA model. The ARMA model has the form as in Eq. (4). The steps: Figure 14 PACF of the squared seasonally starndardized residual series from ARMA(20,1) W. Wangfor et al.: Testing daily flowand modelling autoregressive conditional heteroskedasticity 6 63 4 Residuals 0 2 -2 -4 0 200 400 Day 600 800 (b) Conditional standard deviation 1.0 1.5 2.0 2.5 3.0 (a) 0 1000 200 400 Day 600 800 1000 Fig. 15. A segment of (a) the seasonally standardized residuals from ARMA(20,1) and (b) its corresponding conditional standard deviation Figure 15 Athe segment of model. (a) the seasonally standardized residuals from ARMA(20,1) and (b) sequence estimated with ARCH(21) 0.2 0.2 its corresponding conditional standard deviation sequence estimated with ARCH(21) model (a) (b) 0.0 0.0 ACF 0.1 ACF 0.1 22 0 20 40 Lag 60 80 0 100 20 40 Lag 60 80 100 Fig. 16.Figure ACFs of the standardized and (b) squared standardized fromstandardized the ARMA(20,1)-ARCH(21) model. The 16(a)ACFs of (a) theresiduals standardized residuals and (b)residuals squared residuals from standardization is accomplished by dividing the seasonally standardized residuals from ARMA(20,1) by the conditional standard deviation estimated with ARCH(21). ARMA(20,1)-ARCH(21) model. The standardization is accomplished by dividing the seasonally standardized residuals from ARMA(20,1) by the conditional standard deviation 0.2 2. Fit an ARMA model to the logarithmized and deseasonalized flow series; ACF 0.1 3. Calculate seasonal standard deviations of the residuals obtained from ARMA model, and seasonally standardize the residuals with the first 8 Fourier harmonics of the seasonal standard deviations; 0.0 4. Fit a GARCH model to the seasonally standardized residual series. the squared seasonally standardized residuals are shown in Fig. 12a and Fig. 14, respectively. According to the AIC, as well as the ACF and PACF structure, a GARCH(0,21) model, i.e. ARCH(21) model, which has the smallest AIC value is selected. Therefore, the prelimilary ARMA-GARCH model fitted to the daily streamflow series at Tangnaihai is composed of an ARMA(20,1) model and an ARCH(21) model. The model is constructed with statistics software S-Plus (Zivot and Wang, 2003). 0.2 estimated with ARCH(21). Partial ACF 0.1 1. Logarithmize and deseasonalize the original flow series; 5.2 Model diagnostic and modification 0.2 0.2 CF 0.1 CF 0.1 0.0 If the ARMA-GARCH model is successful in modelling the For forecasting and simulation, inverse transformation (inserial correlation structure in the conditional mean and concluding logarithmization and deseasonalization) is needed. ditional variance, then there should be no autocorrelation left When forecasting, the ARMA part of the ARMA-GARCH 0 20 values 40of the underlying 60 80 time se100 0 the 20 100 standardized in both residuals40and the60squared80residuals model forecasts future mean Lag Lag by the estimated conditional standard deviation. ries following the traditional approach for ARMA prediction, whereas the GARCH part gives forecasts of future volatility, A segment of the seasonally standardized residual seespeciallyFigure over short ries from the ARMA(20,1) model and its corresponding 17horizons. ACF and PACF of seasonally standardized residuals from ARMA(20,1) model Following the above-mentioned steps, a preliminary conditional standard deviation sequence estimated with the ARMA-GARCH model is fitted to the daily streamflow ARCH(21) model are shown in Figs. 15a and 15b. We (a) of (b)residual series from series at Tangnaihai. The ACF and PACF structure standardize the seasonally standardized ARMA(20,1)-ARCH(21) model. The standardization standardization is is accomplished accomplished by by dividing dividing the the seasonally standardized residuals from ARMA(20,1) ARMA(20,1) by by the the conditional conditional standard standard deviation deviation estimated with ARCH(21). W. Wang et al.: Testing and modelling autoregressive conditional heteroskedasticity 0.0 0.0 0.0 0.0 ACF ACF 0.1 0.1 Partial ACF ACF Partial 0.1 0.1 0.2 0.2 0.2 0.2 64 0 0 20 20 40 40 Lag Lag 60 60 80 80 100 100 00 20 20 40 40 Lag Lag 60 60 80 80 100 100 Fig. 17. ACF and PACF of seasonally standardized residuals from the ARMA(20,1) model. 0.2 0.2 0.2 0.2 Figure 17 Figure 17 ACF ACF and and PACF PACF of of seasonally seasonally standardized standardized residuals residuals from from ARMA(20,1) ARMA(20,1) model model (b) (b) 0.0 0.0 0.0 0.0 ACF ACF 0.1 0.1 ACF ACF 0.1 0.1 (a) (a) 0 0 20 20 40 40 Lag Lag 60 60 80 80 100 100 00 20 20 40 40 Lag Lag 60 60 80 80 100 100 Figure 18 of second-residuals and the squared second-residuals Fig. 18. ACFs of (a) second-residuals and (b) the squared second-residuals ARMA(20,1)-AR(16) model. from Figure 18theACFs ACFs of (a) (a) the the second-residuals and (b) (b) from the the squared second-residuals from the the ARMA(20,1)-AR(16) model ARMA(20,1)-AR(16) model the ARMA(20,1) model by dividing it by the estimated conditional standard deviation sequence. The autocorrelations of the standardized residuals and squared standardized residuals are plotted in Fig. 16. It is shown that although there is no autocorrelation left in the squared standardized residuals, which means that the ARCH effect has been removed (Fig. 16b), however, in the non-squared standardized residuals of daily flow significant autocorrelation remains (Fig. 16a). Because the GARCH model is designed to deal with the conditional variance behavior, rather than mean behavior, the autocorrelation in the non-squared residual series must arise from the seasonally standardized residuals obtained in step 3 of the ARMA-GARCH model building procedure. Therefore we revisit the seasonally standardized residuals. It is found that although the residuals from the ARMA(20,1) model present no obvious autocorrelation as shown in Fig. 6a, weak but significant autocorrelations in the residuals are revealed after the residuals are seasonally standardized, as shown by the ACF and PACF in Fig. 17. We refer to this weak autocorrelation as the hidden weak autocorrelation. The mechanism underlying such weak autocorrelation is not clear yet. Similar phenomena are also found for some other daily streamflow processes (such as the daily stream- flow of the Umpqua River near Elkton and the Wisconsin River near Wisconsin Dells, available on the USGS website 23 http://water.usgs.gov/waterwatch), which have strong 23 seasonality in the ACF structures of their original series, as well as their residual series. To handle the problem of the weak correlations, an additional ARMA model is needed to model the mean behaviour in the seasonally standardized residuals, and a GARCH is then fitted to the residuals from this additional ARMA model. Therefore, we obtain an extended version of the model in Eq. (7) as  φ(B)xt = θ (B)εt    εt = σs yt φ 0 (B)yt = θ 0 (B)zt , zt ∼N (0, ht ) ,  q p  P P  2 + ht = α0 + αi zt−i βi ht−i i=1 (8) i=1 where yt is the seasonally standardized residuals from the first ARMA model, zt is the residuals (for notation convenience, we call it second-residuals) from the second ARMA model fitted to yt . An AR(16) model, whose autoregressive order is chosen according to AIC, is fitted to the seasonally standardized residuals from the ARMA(20,1) model of the daily flow series at Tangnaihai, and we obtain a second-residual se- ACF 0.1 ACF 0.1 0 20 40 0.0 65 Lag 60 80 (b) 0 100 20 40 Lag 60 0.0 0.0 ACF 0.1 0.2 (a) ACF 0.1 0.2 0.0 W. Wang et al.: Testing and modelling autoregressive conditional heteroskedasticity Figure 19 ACFs of the (a) standardized second-residuals and (b) squared stand 0 20 40 60 80 0 100 20 40 60 80 100 Lag from ARMA(20,1)-AR(16)-ARCH(21) Lag residuals model. The second-residua Fig. 19. ACFs of (a) the standardized second-residuals and (b) the squared standardized second-residuals from the ARMA(20,1)-AR(16)from are AR(16) to the residuals form ARMA(20,1 ARCH(21) model. The second-residuals obtained fitted from AR(16) fittedseasonally to the seasonallystandardized standardized residuals from ARMA(20,1). p-value Figure 19 ACFs of theThe (a)autocorrelations standardized second-residuals and (b) squared standardized secondries from this AR(16) model. of the 1 second-residual series andARMA(20,1)-AR(16)-ARCH(21) the squared second-residual series residuals from model. 0. 8 The second-residuals are obtained from the ARMA(20,1)-AR(16) combined model are shown 0. 6 in Fig. 18. From visual inspection, weseasonally find that no autocorrefrom AR(16) fitted to the standardized residuals form ARMA(20,1). lation is left in the second-residual series, but there is strong 0. 4 autocorrelation in the squared second-residual series which 1 0. 2 indicates the existence of an ARCH effect. p-value 0. 8 0. 6 0. 4 0 0 5 10 15 Lag 20 25 30 Because the squared second-residual series has similar 2 ACF and PACF stucture to the seasonally 0.standardized residFig. 20. Engle’s LM test for the standardized second-residuals from 0 uals from the ARMA(21,0) model, the same structure of the ARMA(20,1)-AR(16)-ARCH(21) model. Figuremodel, 20 test for second-residuals from 0 Engle’s 20 the 25 standardized 30 the GARCH model, i.e. an ARCH(21) is 5fitted 10 toLM15 Lag the second-residual series. Therefore, the ultimate ARMAAR(16)-ARCH(21) model 6 Conclusions GARCH model fitted to the daily streamflow at Tangnaihai is ARMA(20,1)-AR(16)-ARCH(21), composed of an Figure 20 Engle’s LM test for the standardized second-residuals from the ARMA(20,1)ARMA(20,1) model fitted to logarithmized and deseasonalThe nonlinear mechanism conditional heteroskedasticity in AR(16)-ARCH(21) modelprocesses has not received much attention in the ized series, an AR(16) model fitted to the seasonally stanhydrologic dardized residuals from the ARMA(20,1) model, and an literature so far. Modelling data with time varying condiARCH(21) model fitted to the second-residuals from the tional variance could be attempted in various ways, includAR(16) model. ing nonparametric and semi-parametric approaches (see Lall, 1995; Sankarasubramanian and Lall, 2003). A parametric approach with ARCH model is proposed in this paper to describe the conditional variance behavior. ARCH-type modWe standardize the second-residual series with the conels which originate from econometrics can provide accurate ditional standard deviation sequence obtained with the forecasts of variances. As a consequence, they can be apARCH(21) model. The autocorrelations of the standardplied to such diverse fields as water management risk analized second-residuals and the squared standardized secondysis, prediction uncertainty analysis and streamflow series residuals are shown in Fig. 19. Compared with Fig. 16, the simulation. autocorrelations are basically removed for both the squared and non-squared series, although the autocorrelation at lag The existence of conditional heteroskedasticity is verified 1 of the standardized second-residuals slightly exceeds the in the residual series from linear models fitted to the daily 5% significance level. The McLeod-Li test and the LMand monthly streamflow processes of the upper Yellow River test (shown in Fig. 20) for standardized second-residuals also with the McLeod-Li test and the Engle’s Lagrange Multiconfirm that the ARCH(21) model fits the second-residual seplier test. It is shown that the ARCH effect is fully caused ries well. The small lag-1 autocorrelation in the standardized by seasonal variation in variance for monthly flow, but seasecond-residual series (shown in Fig. 19) is a hidden autocorsonal variation in variance only partly explains the ARCH relation covered by conditional heteroskedasticity. This aueffect for daily streamflow. Among three types of conventocorrelation can be further modeled with another AR model, tional seasonal time series model (i.e. SARIMA, deseasonal24 but because the autocorrelation is very small, it could be neized ARMA and PARMA), none of them is efficient enough glected. to describe the ARCH effect for daily flow, although the the A 66 W. Wang et al.: Testing and modelling autoregressive conditional heteroskedasticity PARMA model is enough for monthly flow by considering season-dependent variances. Therefore, to fully capture the ARCH effect, as well as the seasonal variances inspected in the residuals from linear ARMA models fitted to the daily flow series, the ARMA-GARCH error model with seasonal standard deviations is proposed. The ARMA-GARCH model is basically a combination of an ARMA model which is used to model mean behaviour, and a GARCH model to model the ARCH effect in the residuals from the ARMA model. To preserve the seasonal variation in variance in the residuals, the ARCH model is not fitted to the residual series directly, but to the seasonally standardized residuals. Therefore, an important feature of the ARMA-GARCH model is that the unconditional seasonal variance of the process is seasonally constant but the conditional variance is not. To resolve the problem of the weak hidden autocorrelation revealed after the residuals are seasonally standarized, the ARMA-GARCH model is extended by applying an additional ARMA model to model the mean behaviour in the seasonally standardized residual series. With such a modified ARMA-GARCH model, the daily streamflow series is well-fitted. Because the ARCH effect in daily streamflow mainly arises from daily variations in temperature and precipitation, and given that we have reasonably good skill in predicting weather two to three days in advance (for example, see http://weather.gov/rivers tab.php), the use in developing an ARMA-GARCH model would be limited. However, because (1) on the one hand, the relationship between runoff and rainfall and temperature is hard to capture precisely by any model so far; (2) on the other hand, usually there are not enough rainfall data available to fully capture the rainfall spatial pattern, especially for remote areas, such as Tibet Plateau, and (3) the accuracy of the weather forecasts for these areas are very limited, the ARCH effect cannot be fully removed even after limited rainfall data and temperature data are included in the model. Therefore, the ARMA-GARCH model would be a very useful addition in terms of statistical modelling of daily streamflow processes for the hydrological community. Acknowledgements. We are very grateful to I. McLeod and an anonymous reviewer. Their comments, especially the detailed comments from the anonymous reviewer, are very helpful to improve the paper considerably. Edited by: B. Sivakumar Reviewed by: I. McLeaod and another referee References Bollerslev, T.: Generalized Autoregressive Conditional Heteroscedasticity, Journal of Econometrics, 31, 307–327, 1986. Box, G. E. P. and Jenkins, J. M.: Time series analysis: Forecasting and control, San Francisco, Holden-Day, 1976. Brock, W. A., Dechert, W. D., Scheinkman, J. A., and LeBaron, B.: A test for independence based on the correlation dimension, Econ. Rev., 15, 3, 197–235, 1996. Engle, R.: Autoregressive conditional heteroscedasticity with estimates of the variance of UK inflation, Econometrica, 50, 987– 1008, 1982. Granger, C. W. J. and Andersen, A. P.: An introduction to bilinear time series models, Göttingen, Vandenhoeck and Ruprecht, 1978. Hipel, K. W. and McLeod, A. I.: Time series modelling of water resources and environmental systems, Elsvier, Amsterdam, 1994. Karanasos, M.: Prediction in ARMA models with GARCH-inmean effects, Journal of Time Series Analysis, 22, 555–78, 2001. Lall, U.: Recent advance in nonparametric function estimation – hydrologic application, Reviews of Geophysics, 33, 1093–1102, 1995. Livina, V., Ashkenazy, Y., Kizner, Z., Strygin, V., Bunde, A., and Havlin, S.: A stochastic model of river discharge fluctuations, Physica A, 330, 283–290, 2003. Ljung, G. M. and Box, G. E. P.: On a measure of lack of fit in time series models, Biometrika, 65, 297–303, 1978. McLeod, A. I. and Li, W. K.: Diagnostic checking ARMA time series models using squared residual autocorrelations, Journal of Time Series Analysis, 4, 269–273, 1983. Miller, R. B.: Book review on “An Introduction to Bilinear Time Series Models” by Granger, C. W. and Andersen, A. P., J. Amer. Statist. Ass., 74, 927, 1979. Hauser, M. A. and Kunst, R. M.: Fractionally Integrated Models With ARCH Errors: With an Application to the Swiss 1-Month Euromarket Interest Rate, Review of Quantitative Finance and Accounting, 10, 95–113, 1998. Salas, J. D.: Analysis and modelling of hydrologic time series. In: Handbook of Hydrology, edited by Maidment, D. R., McGrawHill, 19.1–19.72, 1993. Sankarasubramanian, A. and Lall, U.: Flood quantiles in a changing climate: Seasonal forecasts and causal relations, Water Resources Research, 39, 5, Art. No. 1134, 2003. Tol, R. J. S.: Autoregressive conditional heteroscedasticity in daily temperature measurements, Environmetrics, 7, 67–75, 1996. Wang, W., Van Gelder, P. H. A. J. M., and Vrijling, J. K.: Periodic autoregressive models applied to daily streamflow. In: 6th International Conference on Hydroinformatics, edited by Liong, S. Y., Phoon, K. K. and Babovic, V., Singapore, World Scientific, 1334–1341, 2004. Weiss, A. A.: ARMA models with ARCH errors, Journal of Time Series Analysis, 5, 129–143, 1984. Zivot, E. and Wang, J.: Modelling Financial Time Series with SPLUS, New York, Springer, 2003.